{'fraction_non_alphanumeric': 0.03425186328948541, 'fraction_numerical': 0.02575942037592422, 'mean_word_length': 4.478243188287922, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 0, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
On the Spatially Resolved Star Formation History in M51 I: Hybrid UV+IR Star Formation Laws and IR Emission from Dust Heated by Old Stars

We present Lightning, a new spectral energy distribution (SED) fitting procedure, capable of quickly and reliably recovering star formation history (SFH) and extinction parameters. The SFH is modeled as discrete steps in time. In this work, we assumed lookback times of 0-10 Myr, 10-100 Myr, 0.1-1 Gyr, 1-5 Gyr, and 5-13.6 Gyr. Lightning consists of a fully vectorized inversion algorithm to determine SFH step intensities and combines this with a grid-based approach to determine three extinction parameters. We apply our procedure to the extensive FUV-to-FIR photometric data of M51, convolved to a common spatial resolution and pixel scale, and make the resulting maps publicly available. We recover, for M51a, a peak star formation rate (SFR) between 0.1 and 5 Gyr ago, with much lower star formation activity over the last 100 Myr. For M51b, we find a declining SFR toward the present day. In the outskirt regions of M51a, which includes regions between M51a and M51b, we recover a SFR peak between 0.1 and 1 Gyr ago, which corresponds to the effects of the interaction between M51a and M51b. We utilize our results to (1) illustrate how UV+IR hybrid SFR laws vary across M51, and (2) provide first-order estimates for how the IR luminosity per unit stellar mass varies as a function of the stellar age. From the latter result, we find that IR emission from dust heated by stars is not always associated with young stars, and that the IR emission from M51b is primarily powered by stars older than 5 Gyr.


INTRODUCTION
The UV-to-IR spectral energy distribution (SED) of a galaxy, or galactic region, provides useful information about the star-formation history (SFH), which is defined as the star-formation rate (SFR) as a function of time. As byproducts of the SFH, additional useful quantities like the current SFR and stellar mass (M ⋆ ) can be obtained. Obtaining this valuable information from an SED requires modeling the evolution of the stellar populations and their metallicities throughout a long-term SFH, as well as the absorption and emission from gas and dust surrounding the stars and along the line of sight (see Walcher et al. 2011 for a review and http://sedfitting.org).
Given the many components required to model an SED, several simplifying assumptions are often invoked. For example, single burst or exponentially decaying SFHs (with one or two components typically), as well as single metallicities and rigid attenuation curves, are often assumed. Publicly available SED fitting codes have become numerous, each having certain strengths and weaknesses, including, e.g., the treatment of the evolutionary phases of various stars (e.g., the AGB phase) in the core stellar evolution code, modeling of the IR emission from dust, and flexibility in assumptions (see, e.g., Conroy 2013 for a review of the trade-offs). The variability of available SED fitting codes is a product of the needs of the users who developed them, as well as the continuous advancement in computing capabilities and physical models.
In this paper, we present a new SED fitting procedure that is part of a long-term effort designed to understand (1) how hybrid UV+IR star formation laws vary throughout galaxies due to variations in local SFHs and spatial scales; (2) how stellar populations of various ages contribute to the heating of dust in the interstellar medium; and (3) how X-ray binary populations evolve with time with respect to their parent stellar populations. Addressing these scientific areas requires deriving spatially resolved "maps" of SFH information in nearby galaxies using UV-to-IR data, with supplemental high spatial resolution X-ray data (e.g., from Chandra) to associate X-ray binaries with local stellar populations. To this extent, our SED fitting procedure must be both fast, as to obtain SFHs and associated errors over many thousands of pixels, and it must be sufficiently accurate in terms of obtaining SFH parameters and uncertainties consistent with those generated by more detailed SED fitting codes that take longer to implement due to large computational intensity.
Searching parameter space for the best solutions and assigning uncertainties to all model-derived quantities is usually very computationally intensive, with methods ranging from inversions (Heavens et al. 2000;Cid Fernandes et al. 2005), grid based χ 2 minimization (e.g., Mentuch Cooper et al. 2012;Boquien et al. 2016), and an increasing number of models relying on Markov Chain Monte Carlo (MCMC; e.g., Leja et al. 2017). Recently Iyer & Gawiser (2017) introduced what they called the basis dense approach, a procedure testing many SFH functional forms in order to derive a smooth SFH, after many regularization conditions. Here, we implement SFHs as pre-defined, discretized steps in time and utilize a likelihood maximization based on matrix inversions to optimize computational efficiency, similar to Dye (2008), but simpler and faster. Our approach, implemented in the SED fitting code Lightning, is easily hundreds to thousands of times less computationally intensive than many other methods and leads to consistent results. This is primarily due to the fact that all the parameters describing our SFHs are determined by a series of matrix inversions, which usually requires much less computation than other methods, especially when dealing with large dimensional parameter spaces describing the SFH. This speed improvement is paramount when fitting a large number of SEDs across galaxies on a pixelby-pixel basis, as is required by our long-term scientific goals described above.
In the long term, we will apply our method widely; however, here we initially apply our method to data in the Whirlpool galaxy (M51a or NGC5194) and its companion (M51b or NGC5195). M51 is an ideal first test case due to the proximity of the system (8.58 Mpc, which corresponds to a physical scale of 41.6 pc arcsec −1 ; McQuinn et al. 2016), large angular size on the sky (≈8 ′ ×12 ′ ), the nearly face-on orientation of M51a, the extensive multiwavelength coverage in the public archives, and the availability of many comparison results from the literature to test our procedure. In particular, the relatively well-determined timescale for the interaction between M51a and M51b (≈350-500 Myr ago; e.g., Salo & Laurikainen 2000;Dobbs et al. 2010;Mentuch Cooper et al. 2012) allows us to test whether this epoch is recovered from our derived SFHs. This paper focuses on the UV-to-IR SED fitting procedure, and presents first results related to the variations in the star-formation laws across the M51 system and the heating of dust from stellar populations of various ages. In a companion paper, Lehmer et al. (2017, submitted), we present results related to the evolution of the X-ray binary luminosity function with time, which makes use of the SFH maps derived in this paper. The products derived in this paper (e.g., maps of SFH, SFR, M ⋆ , and extinction parameters), as well as the IDL codes that are used here, are provided publicly at the Astrophysics Source Code Library and at https://lehmer.uark.edu/.
This paper is structured as follows: in § 2, we present the dataset and introduce the regions into which we will subdivide M51; in § 3, we describe our SED fitting model and apply it to M51 on a pixel-by-pixel basis, in § 4, we compare our derived parameters with those presented in the literature. In § 5 provide first results on the spatial variations in the UV + IR star formation law and dust heating from stellar populations of different ages; finally, in § 6 we summarize our results.
2. DATA 2.1. Broadband Photometry A total of 18 bands from the FUV to the FIR were used in this work (see Table 1). Central wavelengths, FWHM spatial resolutions, normalized Galactic extinction, and fractional calibration uncertainties are listed. We have combined archival data from several facilities, including GALEX, SDSS, 2MASS, Spitzer, and Herschel. All these data are publicly available and they were downloaded from their respective archives.

Foreground Stars Recognition and Masking
Foreground stars were masked and replaced with the median local background from the FUV to the Spitzer 4.5 µm band. For longer wavelengths, the contribution from any foreground star is negligible.
First we used the find procedure of the NASA Goddard IDL Library (Landsman 1993) to search for point sources over the whole field of the 3.6µm Spitzer image. Then we used the aper routine also from the Goddard Library to performed aperture photometry in all detected point sources in each Spitzer IRAC band (3.6, 4.5, 5.8, and 8.0µm). We used apertures with diameters of 2 times the FWHM, with a sky background determined between radii of 3 and 5 times the FWHM, with the FWHM varying according to each band. Specifically, we used FWHM of 1.7 ′′ , 1.7 ′′ , 1.9 ′′ , and 2.0 ′′ for 3.6, 4.5, 5.8, and 8.0 µm bands, respectively.
Star-forming regions in M51 and foreground stars exhibit very different Spitzer IRAC SEDs. Star-forming region SEDs contain Polycyclic Aromatic Hydrocarbon (PAH) features, which enhance the fluxes in the 5.8 and 8.0µm IRAC bands, while foreground stars in the Galaxy show declining SEDs for increasing wavelengths, corresponding to the Rayleigh-Jeans tail of the stellar emission (F ν ∝ λ −2 ). We therefore recognize as foreground stars all point sources with ratios F ν (5.8µm)/F ν (3.6µm) < 1 and F ν (8.0µm)/F ν (3.6µm) < 1.
Many of these stars are detected in one band, but not in others. For instance, several of them are not detected in the GALEX FUV image and masking these undetected stars could mask real star-forming regions. We therefore only mask a star if it is detected with at least a signalto-noise of 5, with photometry performed in a circle of 2 times the FWHM diameter and background evaluated between 5 and 6 times the FWHM. The FWHM of each band is shown Table 1. We replaced pixels within a circle of 5 times the FWHM diameter centered on each star with the median background pixel around that star in all bands, with the exception of the GALEX bands, where the masked circles had diameters of 3 times the FWHM. Note the masked areas are larger than the apertures used for detection of the cores of the stars.
We tested masking a different numbers of selected foreground stars, but we ended up selecting the 30 brightest stars in the field as a good trade-off between masking bright stars and not blocking star-forming regions.

Convolution to a Common Spatial Resolution
For the Spitzer and Herschel bands we display the nominal values, most commonly used. b FWHM spatial resolution of each map. c Galactic extinction curve normalized at the d Calibration uncertainty as a fraction of the intensity.
After masking foreground stars and replacing them with the local backgrounds, all maps were convolved to a common 25 ′′ FWHM spatial resolution and registered to a common astrometric frame with 10 ′′ pixels. We transform each native point spread function (PSF) to a final PSF of our choice, using the properties of Fourier transforms of convolutions. The procedure is described in Eufrasio (2015) and Straughn et al. (2015) and it is made available to the public through the Astrophysics Source Code Library.
The worst spatial resolution image of our data set is the Herschel SPIRE 250 µm, which has a FWHM resolution of approximately 18 ′′ . We have chosen 25 ′′ to err on the side of caution and guarantee our results do not depend on the uncertainties of the largest PSFs, mainly the 250 µm and the Herschel PACS 160 µm (12.0 ′′ ). The approximate original FWHM spatial resolutions of all our images are tabulated in Table 1. We could in principle have done this analysis with better spatial resolution (∼18 ′′ FWHM) and smaller pixels (6 ′′ −8 ′′ ), but our spatial resolution of 25 ′′ allows us to include Herschel SPIRE 350 µm in any subsequent study and a pixel size of 10 ′′ considerably decreases the number of SEDs to be fitted, making the process more tractable.
After matching spatial resolutions, we combined the Spitzer MIPS 24 µm, the Herschel PACS 70 µm and 160 µm, and SPIRE 250 µm maps in order to generate an IR map corresponding to the integrated TIR luminosity from dust. Our TIR luminosities correspond to the dust emission from 3 to 1100 µm, based on the calibration from Galametz et al. (2013), more precisely where for each pixel of the image L 24 = νL ν (24 µm) and analogously for all bands. This calibration maps these four FIR bands into Draine & Li (2007) models to recover appropriate TIR luminosities. The 1σ statistical uncertainties for the coefficients in Equation 1 are 0.092, 0.024, 0.024, and 0.063, respectively. Since L 70 or L 160 are the larger of the four luminosities, the final TIR uncertainty is dominated by calibration uncertainties (∼ 10%; see Table 1) and not by these statistical uncertainties. The combined TIR surface density map is displayed in the left panel of Figure 1.

Background Subtraction
The galaxies were masked to include all diffuse emission seen around the galaxies in all bands. Background values outside of the galaxy and their uncertainties were determined following the procedures of Eufrasio et al. (2014). The final uncertainties include background subtraction and calibration uncertainties. Calibration uncertainties used for each band are listed in Table 1.
After background subtraction, we confined our analysis to pixels with signal-to-noise ratio greater than 4 in all bands. We tested different criteria, but this fully enclosed our pixels of interest. This left us with a total of 2043 pixels.
Determining dust emission parameters is beyond the scope of this paper and therefore we do not fit wavelengths longer than 4.5 µm. Our final dataset therefore consists of 2043 generalized SEDs with 13 datapoints each (the first 12 broadband images listed in Table 1 and the TIR luminosity).

Regions
We defined larger regions to simplify some of our later analyses. These regions are shown in the right panel of Figure 1. The largest region, corresponding to M51a, is shown in violet. In red, we display the pixels enclosed by the M51b ellipse. Both ellipses are used here as defined by Mentuch Cooper et al. (2012) and their parameters are given in Table 2. In cyan, we show the pixels outside of both ellipses, which comprise the region we will refer hereafter as the "outskirts". We also assess the properties of a region that includes all 2043 pixels combined.
In Figure 2 we show the SEDs for all these regions, normalized at the K S -band. The filled diamonds correspond to the 13 points of our SEDs and the empty diamonds show the four FIR bands combined to produce the TIR luminosity (displayed in the side panel). These SEDs are qualitatively quite different. M51b has the reddest stellar population and, at least apparently, the hottest dust emission in the FIR (from the highest ratio of 70 and 160 µm luminosities). The outskirts is the bluest region in the UV/optical, indicating it was likely strongly affected by the interaction when star formation was triggered in this area. M51a has an intermediate SED between the previous two and the SED combining all pixels is the sum of the three others.
Step Star Formation History With our code Lightning, the SFH of each pixel is modeled with n steps in time, in the form of where t i and t i+1 represent the respective lower and upper boundaries to the i-th time step. For the current work, we have chosen to model the SFH with five steps (i.e., n = 5). The boundaries of the steps have lookback times of 0, 10 Myr, 100 Myr, 1 Gyr, 5 Gyr, and 13.6 Gyr. Figure 3 shows the spectra corresponding to each of the five SFH steps, for ψ i = 1 M ⊙ yr −1 . We have used the population synthesis code Pégase (Fioc & Rocca-Volmerange 1997, 1999, with a Kroupa (2001) IMF and a solar metallicity (Z = 0.02).
The first SFH step, from 0 to 10 Myr, is chosen to model the youngest stellar population. This population is possibly still subjected to a higher birth-cloud attenuation and the only able to emit a substantial number of hydrogen-ionizing photons that quickly recombines and produces hydrogen recombination continua and lines. Its intensity ψ 1 is directly proportional to the total ionizing photon rate and the intrinsic Hα line intensity.
The second SFH step, from 10 to 100 Myr, models the stellar population that, combined with the first step, produces the bulk of the UV emission. Together with the first step it encompasses the typical timescale by which the current SFR is traced by the UV emission, but it still allows the average SFR averaged over the last 100 Myr to be modeled by the two quantities ψ 1 and ψ 2 .
The other three steps were chosen so we had five total SFH steps, as derived by Dye (2008), and similar bolometric luminosities for each of the three steps. Dye (2008), however, did not use UV data and did not impose energy balance between the absorbed UV/optical and dust-reprocessed IR, which considerably improves SED fits, as we will discuss later. Leja et al. (2017) used six SFH steps to model the integrated SED of the 129 nearby galaxies presented in Brown et al. (2014). They were able to recover all six SFH steps relatively well from broadband SEDs alone. They have also observed remarkable agreement between the derived and observed Hα, even when the line emission is not used in the fits. Table 3 lists the stellar masses, bolometric luminosities, and the rate of hydrogen-ionizing photons Q 0 for each SFH step per 1 M ⊙ yr −1 . The extinction-free high-resolution model spectrum, L ν (ν), is therefore given by the sum of the spectra emitted from all steps of the SFH,L i ν (ν), i.e., whereL i ν (ν) is a specific luminosity per unit SFR, or the template spectrum arising for the i-th SFH step per 1 M ⊙ yr −1 , as depicted in Figure 3.

Dust Extinction Parameters
The high-resolution model spectrum is attenuated by an intrinsic extinction curve and a fixed Galactic extinction curve. For the Galactic extinction, we use  Figure 1. The color code is the same as in Figure 1 Table 2 (i.e., violet, red, cyan, and black respectively correspond to M51a, M51b, outskirts, and all pixels combined). Uncertainties of 10% were used just for plotting purposes, once the combined uncertainty is mostly coming from calibration uncertainties. Filled diamonds correspond to the 12 broadband filters and the TIR luminosity (shown in the side panel). Empty diamonds show the four FIR bands combined to produce the TIR map. Black curves correspond to the response functions of the broadband filters used in the SED fitting. Response function peaks were normalized to 0.05. Figure 3. Five spectral templates corresponding to the 5 steps of the SFHs considered in this paper. High-resolution Pégase spectra are shown in conjunction with the corresponding broadband SEDs. These are the intrinsic, unattenuated spectra, for a constant SFR of 1 M ⊙ yr −1 . Only the first step (ages from 0 to 10 Myr) emit nebular continua and line emission. Table 3 lists some properties associated with each of these spectral templates. Response functions of all filters are shown in black and their peaks were normalized to 10 7 L ⊙ .
the Fitzpatrick (1999) standard curve as coded in the fm unred script of the NASA Goddard IDL Library (Landsman 1993).
For the intrinsic attenuation, we consider a variable extinction curve composed of a diffuse attenuation com-ponent, common to steps of all ages, Calzetti et al. (2000) extinction curve, as modified by Noll et al. (2009) with the addition of a bump and a δ parameter controlling the UV slope. The D(λ) term adds a Drude profile representing the 2200Å bump of the extinction curve, of the form where λ 0 = 2175Å is the central wavelength of the bump and ∆λ = 350Å its FWHM. In this work, we take the bump strength E b to be related to the slope δ as following the findings of Kriek & Conroy (2013). On top of the diffuse extinction, a birth-cloud attenuation is applied to the youngest SFH step, i.e., lookback times up to 10 Myr. This extinction is given by Thus, the attenuated, high resolution spectrum is , where δ i1 is the Kronecker delta, i.e., δ i1 equals unity for i = 1 and 0 otherwise, not to be confused with the δ extinction parameter defined in Equation 4. Our model has therefore three free parameters that control extinction: τ DIFF,V , δ, and τ BC,V .
After extinction is applied, the high resolution attenuated spectrum is convolved with the broad or narrowband filters to produce the model SED that will be compared to the observations. The full SED model is then with R k (ν) being the response curve of the k-th filter, normalized to 1, i.e., ∞ 0 R k (ν) dν = 1. The TIR dust emission provides a strong constraint on the attenuation, since the UV-optical attenuated light will be re-emitted in the IR by the absorbing dust particles. This paper aims at recovering SFHs and not the dust parameters, as done by Mentuch Cooper et al. (2012) and others, so we do not model the dust emission in detail. Instead, we have combined four FIR bands (24, 70, 160, and 250 µm) in a single total luminosity map, following Equation 1.
In order to better constrain the attenuation parameters, we model the TIR luminosity as the total intrinsic (extinction-free) luminosity minus the attenuated luminosity. The model TIR luminosity L i TIR is separately computed for each SFH step, given the three extinction parameters, as and the model TIR luminosity is then The set of observables to be fitted are therefore all the broadband luminosities as well as the TIR luminosity, i.e., so that, for any band k, L i k is therefore the intensity of the i-th step in the kth band, with k varying from 1 to m + 1, with k = m + 1 representing the TIR luminosity of that given step per unit SFR, and a given set of the three extinction parameters.

Likelihood Maximization via Matrix Inversion
For a given choice of extinction parameters, the model SED is a linear combination of the SEDs of the SFH steps, as in Equation 13. This greatly simplifies χ 2 minimization, as we show in the following lines, increasing computational speed for a large number of computations.
We measure the goodness-of-fit as Here L obs k are the observed luminosities and σ k are their corresponding 1σ uncertainties.
To minimize χ 2 , we set its partial derivatives with respect to each SFH step intensity ψ i equal to zero, i.e., which can be rearranged as and simplified to where is calculated from the population synthesis models alone and is calculated from a combination of models and observed luminosities. This procedure leads to a linear system of n equations, that is, with Equations 17 representing as many equations as SFH steps. Hence the whole system of equations can be written in matrix form as where A is a square (n × n) symmetric matrix, with only positive elements, ψ is a column vector with all the SFH coefficients ψ i , and B is column vector of dimension n as well. A is invertible if det(A) = 0 and in that case there is only a single solution for the system of equations. The solution for the SFH is thus given by Since the coefficients ψ i represent the SFRs of the different previously defined SFH steps, they are necessarily non-negative real numbers. However, the numerical inversion in Equation 21 does not require that ψ i ≥ 0, for all steps i. The solution may initially contain one or more negative ψ i coefficients after the first matrix inversion.
To assure all SFH steps are non-negative, we have implemented an iterative algorithm that first sets to zero all coefficients that are found to be negative, next reduces the dimensionality of the problem, and finally solves Equation 21 again for the other coefficients until they are all non-negative.
In order to determine the uncertainties in the SFH solutions, we simulate 400 perturbed SEDs and solve for the SFH each time. For each simulated perturbed SED, the fluxes of each band are drawn independently from Gaussian distributions set by the mean and standard deviation of the measured fluxes at each wavelength.
Lightning, our SED fitting code, is available to the astronomy community through the Astrophysics Source Code Library.

Extinction Parameters
As described above, we solve for the SFH step intensities by employing the matrix inversion procedure described in § 3.3. To determine the extinction parameters and their uncertainties, we create a three-dimensional grid, with 81 steps in each dimension with as in Leja et al. (2017). We adopt flat priors for all three parameters, and also impose the joint prior τ DIFF,V /τ BC,V < 2.0 used by Leja et al. (2017), to reproduce the previous observations suggesting that the total optical depth towards younger stars that ionize the nebular emission lines emitting gas is around twice that of the older stellar component (Calzetti et al. 1994;Price et al. 2014).
For each of the 531,441 (i.e., 81 3 ) discrete points of extinction parameter space grid, we employ the inversion procedure to determine the five ψ i coefficients of the SFH from the observed SED. For each of the 2043 pixels, we calculate the minimum χ 2 for each point of this 3D grid of models.
With the current grid, we have performed both a Bayesian analysis in 3D extinction parameter space to derive marginalized likelihoods and posterior probability densities for each extinction parameter and a frequentist analysis repeating this minimization procedure for each of the 400 simulated SEDs, perturbed around the observed SED. No statistically significant difference was observed between the two approaches regarding the confidence intervals for each of the three extinction parameters.
Henceforth, we quote only frequentist uncertainties and confidence intervals, resulting from the fits of the 400 simulated maps we generated for each of the 2043 observed SEDs. Every quantity derived from the fits is measured 400 times, by selecting the model that maximizes the posterior probability distribution and not necessarily the model with lowest χ 2 statistic. No significant difference is observed between the minimum χ 2 and the χ 2 for the model that maximizes the posterior probability density. For any derived quantity, we therefore quote the median value of the 400 simulations as the solution, with the lower and upper error bars corresponding to the differences to the 16th and 84th percentiles, respectively. . Goodness-of-fit (χ 2 ) map, corresponding to the models that maximize the posterior probability for each pixel. Ellipses are the same as those in Figure 1

Derived Quantities
We have applied our SED fitting procedure to the 13 images of M51, which were registered to a common grid with 2043 square pixels of 10 ′′ (see § 2.3 for details). Figure 4 displays the χ 2 map for the model that maximizes the posterior probability for each pixel. Median χ 2 considering all pixels is 6.8, which is very reasonable, considering we have 13 observables, 3 dust parameters, and 5 SFH steps (5 degrees of freedom). Figure 5 presents the SED fitting of the pixel with median χ 2 , with a highest posterior probability model having χ 2 = 6.8. This exemplifies a typical fit. The observed SED is shown as the empty black diamonds, with the last point in the side panel corresponding to the TIR luminosity. All error bars on the black points are 1σ uncertainties and include background subtraction and calibration uncertainties. The attenuated SED model is shown in red, and in blue we show the intrinsic, unattenuated SED model. The red point fitting the TIR luminosity is the difference in bolometric luminosity between the (highresolution) intrinsic and attenuated model spectra, as in Equation 10.
The results for the four regions defined in Figure 1, right are summarized in Table 4. The table lists all five SFHs step intensities, the average SFR averaged over the last 100 Myr (SFR100 = 0.1ψ 1 + 0.9ψ 2 ), the total stellar mass M ⋆ enclosed in each region, and the specific star formation rate (sSFR), sSFR100 ≡ SFR100/M ⋆ . Figure 6 shows the SFHs of these four regions, using the same colors as the right panel of Figure 1. A very interesting feature is evident from the outskirts region, implying that the peak of star formation activity happened at the third SFH step, i.e., corresponding to stellar ages between 100 Myr and 1 Gyr. M51b has a predominantly old stellar population, with most of its mass com-  Figure 5. SED of the pixel with the median χ 2 (χ 2 med = 6.8) of all of our pixels. The observed luminosities are shown as the black empty diamonds. The observed TIR luminosity, combining the 24, 70, 160, and 250 µm data is displayed in the side panel. The error bars correspond to the 1σ uncertainties, which include background subtraction and calibration uncertainties. In red we display the model with maximum posterior probability and in blue we show its intrinsic SED (unattenuated). The filter response functions are shown as well. The attenuated power corresponds to the red model that fits the observed TIR luminosity.
ing from the oldest SFH step and considerably little SF activity over the last 100 Myr. These results are consistent with the findings of Mentuch Cooper et al. (2012), who modeled the SFH across the M51 system with two exponential decays and recovered an older stellar component with ages between 7 and 10 Gyr and a younger component with ages between 300 and 600 Myr, which they attributed to the interaction between the two galaxies. Figure 7 displays the maps of the solutions for the three extinction parameters. The left panel shows the diffuse optical depth at the V -band, τ DIFF,V , which was applied to SFH steps of all ages. In general, attenuation is less than unity and it increases slightly where star formation is more intense in the last ∼100 Myr. The middle panel displays the parameter δ that controls the slope of the extinction curve. The solutions favor curves around Calzetti extinction curve (δ = 0) and strongly disfavor much steeper solutions (δ < −0.4). The right panel Figure 6. Five-step SFHs of all pixels combined (black), the M51a ellipse (violet), the M51b ellipse (red), and pixels outside the two ellipses (cyan). The solid lines denote the median values from the simulations and dotted lines show the 16th and 84th percentiles, enclosing the hashed areas. The first SFH step is actually from 0 to 10 Myr, but it is displayed here from 1 to 10 Myr for simplicity. The pixels outside of the two ellipses (blue) correspond to the outskirts of the spiral galaxy M51a and the interaction region between M51a and M51b. Their combined SFH points to an enhancement of SFR not associated with the recent SFR (last 100 Myr), corresponding to ages between 100 Myr and 5 Gyr (SFH steps 3 and 4). M51b seems to have the bulk of the oldest stellar population (fifth SFH step; ages between 5 and 13.6 Gyr) and a declining SFH. Also, the M51b ellipse does not contain considerable amounts of stars produced by the interaction. The spiral galaxy (M51a), however, contains mainly stars produced by the third and fourth SFH steps, likely strongly affected by the interaction.
shows the extra V -band optical depth for birth clouds, τ BC,V , that applies only to the most recent SFH step, affecting the stellar population younger than 10 Myr. This component of the attenuation is also small, usually less than unity. It, however, rises to high values in the interarm regions, which does not mean the full attenuation in the inter-arm regions is high. Figure 8 shows the final recovered attenuations at the SDSS g and GALEX FUV bands. The combined g-band attenuation traces the arms and star forming regions very well, but in the FUV the solution is much messier. This is likely due to larger propagated uncertainties in the FUV attenuation, Figure 7. Maps of the three extinction parameters. The left panel displays the optical depth at the V -band, τ DIFF,V , which is applied to SFH steps of all ages. In general attenuation is less than unity and it increases slightly where star formation is more intense in the last ∼100 Myr. The middle panel displays the parameter δ that controls the slope of the extinction curve. The solutions favor curves around Calzetti extinction curve (δ = 0) and strongly disfavor much steeper solutions (δ < −0.4). Right panel shows the extra V -band optical depth for birth clouds, τ BC,V , that applies only to the most recent SFH step, affecting stellar population younger than 10 Myr. This component of the attenuation is also small, usually less than unity. It, however, rises to high values in the inter-arm regions, which does not mean the full attenuation in the inter-arm regions is high. Figure 8. Maps of the combined attenuation at the SDSS g and GALEX FUV bands. The recovered g-band attenuation traces arms and star-forming regions quite well. The FUV attenuation map is much messier. This can happen because of the higher uncertainties on the recovered FUV attenuation or because it indeed has a more complicated morphological structure.
but also because the FUV attenuation does have more complicated morphological structure.
Maps of the SFR densities Σ(ψ i ) at different epochs are shown in Figure 9, in units of M ⊙ yr −1 kpc −2 , as well as the average SFR density in the last 100 Myr, Σ(SFR100). The figure shows the stellar mass of M51b comes mainly from the 5th SFH step (ages from 5 to 13.6 Gyr), as expected from Figure 6, but M51a also has an old component mainly in its center. M51a experienced a lot of SF from between 100 Myr and 5 Gyr ago and the stars produced during this period are well mixed throughout the disk.
Interestingly, in the last 100 Myr, the timescale usually adopted to infer the current SFR from the UV emission, the SFR of M51 has been relatively weak. This implies a considerable fraction of the light in all photometric bands (including the UV) is arising from older stars, not associated with this "recent" SFR.
Our model naturally separates the contributions of each SFH step to any quantity derived by our model. We are able to, for example, quantify that the third SFH step (ages between 100 Myr and 1 Gyr) contributes 22.9 +8.3 −7.5 % of the total intrinsic FUV. The second step (10-100 Myr), which would naively be expected to have the most impact on the total FUV luminosity, contributes 24.1 +13.6 −8.1 %. We also estimate that 11.8 +4.8 −3.3 % of the stellar mass of the whole map arises from the third SFH step, very consistent with the burst stellar mass fraction of between 5-15% derived by Mentuch Cooper et al. (2012). This shows the effect of the interaction on the system as a whole. Our model indicates 36.3 +9.6 −14.3 % of the current stellar mass coming from the fourth, and 51.5 +13.9 −10.2 % from the fifth, oldest step. For completeness, only 0.049 +0.019 −0.017 % and 0.34 +0.28 −0.10 % of the total current stellar mass were produced by the first and second SFH steps, respectively.

Lightning Parameter Recovery
We performed simulations of parameter recovery for a variety of theoretical SFHs, including scenarios where SFR increases or decreases with lookback time, as well as step functions where only one SFH step is non-zero, and cases with a sharp drop. As an example of parameter recovery and of the degeneracies involved, we use the best-fit SFH and extinction parameters derived for the large region M51a as simulation inputs. This large region (with 862 pixels) is detected with high confidence over the background noise in all images and its SED uncertainties arise mainly from calibration uncertainties on the photometry (see Table 1).
The observed SED was fitted once, then we simulated 400 SEDs around the best-fit model SED. The simulated SEDs were generated as described in § 3.3. Figure 10 is one of the diagnostic outputs of Lightning. It shows the results of this simulation. First, histograms of all 8 parameters were plotted. Then the maximum value of Figure 9. Maps of the SFR surface density at different epochs. All maps in units of M ⊙ yr −1 kpc −2 . Maps correspond to the highest likelihood models, but even for pixels where the best solution is zero, we can still provide upper limits and confidence intervals from our simulations. each x-axis was used to discretize the full 8D parameter space. For each of the five SFH intensities, we created 1D grids with with 15 values, linearly increasing from 0 to the respective maximum value. The 3D extinction grid has respectively 21, 14, and 21 values for τ DIFF,V , δ, and τ BC,V . The whole 8D grid has therefore ≈4.7 billion models (15 5 × 21 × 14 × 21).
In Figure 10 we also present histograms of the 400 results for the best-fit χ 2 values, SFR100, stellar mass, bolometric luminosity, and sSFR100 are also shown, with the median results as blue dashed vertical lines and the best fits to the data as red vertical lines. Figure 10 also includes a set of scatter plots with the results of best-fits to the 400 simulations of the SED (clouds of points), the median values recovered from the simulation (empty diamonds), the best-fit to the data (solid, red stars), and the various Bayesian marginalized likelihoods for eight parameters (five SFH step intensities and three extinction parameters, shown as the red contours and curves). The three contour levels enclose 68%, 95%, and 99% of the 2D marginalized likelihoods.
Various degeneracies are at play. Figure 10 shows how well each parameter and derived quantity is recovered both from the Bayesian and from the frequentist analysis. The results from both methods agree, with the best-fit to the data (shown as red stars or red vertical lines) matching, within uncertainties, the median values from the 400 simulations (shown as blue diamonds or dashed vertical lines). Many quantities, however, are derived effectively as upper limits, in agreement with the findings of Leja et al. (2017), who used a Markov Chain Monte Carlo (MCMC) method. For this particular SED, ψ 3 and τ DIFF,V are well recovered and both distributions agree for these parameters. The final likelihood is almost independent of τ BC,V , with lower values around 0.5 being weakly favored, and the favored extinction curve is flatter than the Calzetti curve (δ ∼ 0.2). All other parameters (ψ 1 , ψ 2 , ψ 4 , and ψ 5 ) receive proper upper limits and the best-fit to the data agree with the median from the simulations, within uncertainties. The χ 2 distribution shows the recovered goodness-of-fit for the data is reasonable, with the model not overfitting (χ 2 too low) nor being inadequate (χ 2 too high).
We also show the histograms for the derived SFR100, stellar masses, bolometric luminosities, and sSFR100, as derived from the 400 simulations. We do not show these marginalized Bayesian distributions, as computing them is too computationally intensive. In fact, one of the main advantages of Lightning is generating a vector of derived quantities with the same size of the number of simulations, not having a burning phase as it would happen to MCMC methods.
In this work, Lightning was extensively tested with real data, considering the large variety of SEDs fitted here (2043 pixels and 4 large regions). In the next section, we compare the results of Lightning with a widely tested SED fitting code.

Comparison to CIGALE
We have compared our results with the well-tested, widely used Python implementation (Roehlly et al. 2014;Boquien et al. 2014) of the CIGALE SED modeling code (Burgarella et al. 2005;Noll et al. 2009), PCIGALE. We ran the latest version (0.11.0) of the code on the SEDs of all 2043 pixels and of the four separate large regions.
Initially, we ran CIGALE with the same grid as Boquien et al. (2014). Bruzual & Charlot (2003) stellar population models were used, assuming a Chabrier (2003) and a solar metallicity (Z = 0.02). The SFH consists of two exponentially decaying populations, an old stellar population with formation age of 13 Gyr and a young stellar population with eight possible formation ages (5,10,25,50,100,200,350,or 500 Myr). For a complete list of the parameter values used in CIGALE, see Table 2 of Boquien et al. (2016). The FIR SEDs are fitted with Dale et al. (2014) templates. Excellent agreement is found between the TIR luminosities derived from Galametz et al. (2013) calibration based on Draine & Li (2007) models (Equation 1) and the ones derived with CIGALE, with the CIGALE values being only 5% larger (less than the calibration uncertainties) on average, considering all 2043 pixels. With this initial model, however, we observe a significant discrepancy between the stellar masses recovered from CIGALE and Lightning, with the best-fit masses from CIGALE being 47% higher on average, considering all 2043 pixels. We suspected this disagreement was most likely due to the fixed 13 Gyr age for the older stellar population used by Boquien et al. (2016), which has considerably lower mass-to-light ratio. We added more freedom for the possible ages of the older stellar population. CIGALE accepts at most three ages for the old stellar population. We then chose 5, 8, and 13 Gyr as the new age grid for the old stellar population. We also added 1 and 2 Gyr as possible ages for the young stellar population. After this simple modification, reasonable agreement is reached between the stellar masses derived from Lightning and CIGALE, with the pixel-by-pixel stellar masses from CIGALE being only 2.2% higher on average.
To compare the SFHs derived from CIGALE with those derived from Lightning, we averaged individually the high-resolution SFHs of all 2043 pixels over the five time bins considered here. Figure 11 shows the results of this comparison adding the SFHs over the four large regions considered in this work. SFHs for M51a, M51b, outskirts, and all pixels combined are shown, respectively, in violet, red, cyan, and black. High-resolution models from CIGALE are shown as the solid lines and the thicker solid lines show the high-resolution CIGALE models averaged over our five time bins. Note many exponential decays are seen in a single region's SFH model. This happens because there are hundreds of pixels in each region, with the SFH of each pixel containing two exponential decays of different decay timescales. The number of pixels in each region is actually listed in Table 2. 1-σ confidence intervals (between the 16th and 84th percentiles) from Lightning are shown as the hashed areas and the median values are shown as the dotted lines.
Note these are the best-fit SFHs from CIGALE, and no uncertainties were derived on these high-resolution SFHs, for the simple reason that effort would be too computationally intensive. Qualitatively, however, similar trends are recovered from the two models. Considering CIGALE's uncertainties are of the same magnitude as the ones derived from Lightning, there is virtually no significant discrepancy. For example, one apparently large disagreement is the 0-10 Myr SFR in M51b, which is lower by a factor of 3 or 4 in Figure 11. However, with the simpler, initial CIGALE grid the Bayesian SFR10 (averaged over the last 10 Myr) is 0.114±0.018 M ⊙ yr −1 , while the best-fit SFR is 0.16 M ⊙ yr −1 , and therefore in excellent agreement with our estimate. We have good reason to believe that a more detailed CIGALE model accepting more mass ratios between the younger and the older would bring CIGALE's answer closer to Lightning's. This increase in complexity to fix such a small detail is beyond the scope of this work. Since the SFH steps were chosen here to have similar bolometric luminosity per unit SFR, the large stellar mass in the older SFH bin for M51b outshines the younger SFH bin by at least a factor of 10. Of course, in the FUV the youngest first SFH step still outshines the oldest, fifth step by a factor of 40 and the SFR can be recovered.
All in all, remarkable agreement is found between the results of CIGALE and Lightning, especially considering the two codes use different evolutionary population synthesis models, (slightly) different IMFs, dust extinction treatment, and fitting techniques. Regarding running time, however, for an extinction grid with as many models as CIGALE, Lightning performs ∼2000 times faster. We point out CIGALE models the IR SED and Lightning is using the TIR luminosities from Galametz et al. (2013), which is responsible for part of the speedup. Lightning derives fast and reliable SFHs, which meets our current goal, but we plan to add IR SED fitting in a later version of Lightning.

Hybrid Star Formation Rate Across M51
Here we set out to use our SED fitting results to derive hybrid SFR laws combining UV and IR indicators. As an example, we show the GALEX FUV + Spitzer 24 µm, but the reasoning used here can be applied to all other hybrid SFR laws.
The UV has the convenient property of tracing almost exclusively young stellar populations. If a galaxy, or a region of a galaxy, forms stars at a constant rate, its UV luminosity reaches a plateau and that plateau is commonly used as a SFR calibration. The FUV reaches a steady state (i.e., constant L UV /SFR) on timescales of 100 Myr. The same can be said for tracers of ionizing photon rate, like Hα, Hβ, and Paschen α line emission intensities, hydrogen recombination continuum intensities, or free-free intensities, but that plateau timescale is much shorter, usually less than 10 Myr. Qualitatively, these timescales correspond to lifetimes of the least massive stars that produce significant UV or ionizing photons. Eufrasio (2015) tabulated these timescales for different bands, for various Figure 10. Diagnostic plots for the fitting of the SED of the large region M51a. Scatter plots of the values recovered from the 400 simulated SEDs are shown for the two-by-two combinations os parameters and histograms are shown for all eight parameters, SFR, stellar mass, bolometric luminosity, sSFR, and best-fit chi-squares. The best-fit values to the observed data in the histograms as the red, solid, vertical lines (or red stars) and the median values from the 400 simulations are shown as the blue, dashed, vertical lines (or blue diamonds). On top of the histograms for all eight parameters we display Bayesian 1D marginalized likelihoods. All two-by-two combinations of 2D scatter plots are shown with the 2D marginalized likelihoods overlaid as the red contours (enclosing 68%, 95%, and 99% of the likelihood).
IMFs, and metallicities. One main issue that arises when modeling SFHs is that an inappropriate choice of models might not be able to properly recover independently the UV and Hα fluxes. Interestingly, Leja et al. (2017) shows that one is able to recover the proper Hα fluxes by solely modeling broadband SEDs with step SFHs, very similar to the ones used here.
Based on the above, the UV emission is a strong tracer of the recent SFR, and the conversion from a given UV band intensity to SFR can easily be derived from stellar population synthesis. These conversions depend on IMF and metallicity, as well as the timescale considered for the SFR scenario. From the population synthesis code Pégase, assuming a Kroupa IMF, a stellar population with solar metallicity (Z = 0.02) that formed stars uniformly during the last 100 Myr (and nothing before) will More generally, one can write where k X is the conversion from the intrinsic luminosity of band X to the average SFR over the last 100 Myr (Kennicutt 1998;Kennicutt & Evans 2012;Murphy et al. 2011). Of course, these conversions rely on the assumption of constant SFR over the last 100 Myr and nothing before it. Looking at Figure 6, however, one might expect that this is not the case for any of our regions. Actually for the FUV calibration conversion k FUV (Equation 25), instead of the factor 1.6 for the simple, idealized case, the SFHs derived here for M51a, M51b, outskirts, and all regions lead to conversion factors of 1.29 +0.89 −0.39 , 0.81 +0.39 −0.17 , 1.33 +0.66 −0.30 , and 1.28 +0.83 −0.37 , respectively. For the NUV band, instead of 2.7, the conversion factors are 1.72 +1.22 −0.52 , 1.06 +0.51 −0.21 , 1.77 +0.86 −0.39 , and 1.71 +1.17 −0.50 . Thus the simplified assumption yields a k X higher by 40%-80% for M51, which alone would overestimate the SFR. Note that the conversion being different is not a surprising result of our fits, but is expected from a theoretical standpoint, if the SFH is anything other than a constant SFR for the last 100 Myr preceded by no SFH. It just happens that the SFHs across M51 are quite different from the simplifying assumption. We therefore emphasize the importance of properly modeling the SFH and attenuation parameters in order to derive reliable SFR calibrations.
One difficulty of using the UV luminosities as SFR indicators is that the UV fluxes are highly reprocessed by dust and gas, mainly due to dust attenuation and gas absorption and scattering. Some modeling is necessary to obtain intrinsic UV luminosities. SED fitting, as employed here, with as many observables as possible, is most likely one of most accurate ways of determining dust attenuation and emission.
Standard prescriptions to correct the UV luminosities for dust extinction are widely available, empirically calibrated for large samples of galaxies. For instance, some commonly used include: and L corr NUV = L obs NUV + 2.26 × L obs 24µm , as derived by Hao et al. (2011). After applying these prescriptions, the so-called "corrected" luminosities on the left-hand sides of Equations 28-31 are treated as the intrinsic. These luminosities can then be converted into SFRs following the appropriate relations. Generalizing the equations above, we can define where a corr (X, Y ) is the correction factor multiplying observable luminosity L Y in order to account for the extinction in observable luminosity L X . In fact, a corr (X, Y ) is the ratio of the attenuated light in the X band over the intensity of the Y band, i.e., the efficiency of attenuation in X compared to emission of Y .
In Eufrasio et al. (2014), we tested Equation 32 combining the observed GALEX FUV and WISE 22 µm emission (which is almost numerically identical to Equation 28). The correction factor a corr (FUV, 22 µm) was observed to vary significantly across the giant interacting Condor galaxy, NGC 6872, for regions 10 kpc in diameter. a corr (FUV, 22 µm) were also observed to increase with increasing sSFR, meaning younger, bluer regions were observed to have higher attenuation at FUV wavelengths compared to emission at 22 µm. However, FIR maps were not available for the Condor galaxy and energy balance was not imposed on those SED fittings. This made it more difficult to properly model the intrinsic FUV. Boquien et al. (2016) modeled the spatially resolved SFHs of eight nearby face-on spiral galaxies with CIGALE and derived similar relations, where a variety of these correction factors, a corr (X, Y ), are clearly increasing with increasing sSFR, especially a corr (FUV, TIR). However, they seemed to also be correlated with the attenuations, stellar masses, and SFRs, albeit more weakly, and it was extremely difficult to determine if the trend was driven by the underlying stellar population or by geometry and differential attenuation effects.
These hybrid SFR tracers are widely used, since they are often the only way to determine SFRs for galaxies in the distant universe. Hence it is important to properly determine how the correction factor depends on attenuation and on the underlying stellar population properties.
Reliably measuring these correction factors requires good estimates of the intrinsic UV luminosity, which makes our model very suitable to tackle this problem. Taking advantage of SFRs derived from our previously described SED modeling, we produced maps of various correction factors, a corr (X, Y ), combining the GALEX bands with the four FIR bands available for M51 (24, 70, 160, 250 µm) and the TIR map. Figure 12 displays a corr (FUV, 24 µm) versus sSFR100 = SFR100/M ⋆ = Σ(SFR100)/Σ(M ⋆ ) for all our pixels. Empty violet diamonds correspond to pixels inside the M51a ellipse, empty red triangles are inside the M51b ellipse, and cyan crosses are pixels in the "outskirts" region. The dashed and dotted lines correspond, respectively, to the Hao et al. (2011) andEufrasio et al. (2014) relations. Treating all pixels in a region as a single location, we can also plot global quantities for our regions. The larger, filled symbols in the figure were added for the global quantity derived for each region, where the filled diamond, triangle, cross, and circle correspond to M51a, M51b, outskirts, and all regions combined, respectively. The eight galaxies from Boquien et al. (2016) are shown as filled yellow stars.
The problem of estimating recent SF activity from a combination of a UV and a FIR band is very complicated, with a corr (FUV, 24 µm) increasing on average with increasing FUV attenuation or SFR100, and increasing with decreasing stellar mass density. The fact that a corr (FUV, 24 µm) changes at all suggests that stars older than 100 Myr (not associated with the recent SFR traced by the UV) are capable of heating dust and significantly contribute to all FIR bands and to the TIR luminosity.
In fact, Kennicutt et al. (2009) arrived to the conclusion that a considerable fraction of the TIR radiation in normal star-forming galaxies is heated by stars older than 100 Myr. Their comparisons with evolutionary synthesis models suggested that up to 50% of the TIR emission could be from dust heated from this evolved stellar population. They issued a stern warning that hybrid SFR indicators calibrated by them, Calzetti et al. (2007), and Zhu et al. (2008) should not be applied to map the spatially resolved SFR in galaxies without the risk of introducing significant and possibly large systematic errors in the resulting SFR maps. These systematic uncertainties are exactly what we characterize across M51 in Figure 12.
Our model is perfectly suited to quantitatively check the hypothesis that old stars contribute to the FIR emission, since by construction we have five steps of the SFH probing different epochs. In the next section, we decompose the FIR emission into contributions from the five SFH steps.

Decomposition of the IR Emission Into
Contributions From The SFH A significant body of work is available in the literature on the contribution of old stars to the IR emission of nearby galaxies. Various Herschel studies show warm dust to be associated with the average SFR of the Figure 12. The correction factor acorr(FUV, 24 µm) versus sSFR100 for all pixels of our map. Violet diamonds correspond to pixels inside M51a elipse, red triangles are inside M51b elipse, and cyan crosses are pixels in the "outskirts" region. Global values for the three regions are displayed in larger, filled symbols. Dashed horizontal line shows Hao et al. (2011) correction factor and dotted line represents Eufrasio et al. (2014) relation. The eight galaxies from Boquien et al. (2016) are shown as filled stars. The fact that this correction factor changes with sSFR100 indicates the old stellar population, not associated with the recent SFR100, significantly contributes to the IR emission.
last 100 Myr and colder dust to be associated to older stars (Popescu et al. 2000;Bendo et al. 2010Bendo et al. , 2012Bendo et al. , 2015Boquien et al. 2011;De Looze et al. 2012a,b;Lu et al. 2014;Natale et al. 2015). From hydrodynamical simulations, even in actively star-forming galaxies, up to a third of the TIR luminosity is associated with stars older than 100 Myr (Boquien et al. 2014(Boquien et al. , 2016. There is also strong evidence that the FUV−NUV color is strongly affected by the presence of old stars Boselli et al. (2005); Cortese et al. (2008).
Disentangling IR emission from all ages is crucial to determine what fraction is associated with the recent SFR. Crocker et al. (2013) investigated the SFH effect on the TIR fraction heated by stars older than 100 Myr. They found the result strongly depended on the SFH and extinction prescription used. For an exponentially decaying SFH of decay timescale of 5 Gyr about half the TIR would be heated by old stars, while for a constant SFH, that number would be about a third.
From our model, we can directly determine the fraction of FIR emission coming from each SFH step of each pixel, since the TIR emission was one of the points of our generalized SEDs and it was fit using the linear combination of the five SFH components. However, to better understand the contribution of the old stellar population to the IR emission, we can directly decompose each of the observed FIR images, as well as the TIR map, into contributions from our five SFH steps.
We implemented our procedure using each of the large regions (i.e., M51a, M51b, outskirt, and whole map). We performed a χ 2 minimization using the same inversion method used in our SED fittings (see §3.3). This time, however, we assumed that any of the observed IR band luminosities L obs Y were generated by a linear combination of the five derived SFH maps.
For instance, the model intensity of the j-th pixel in any region S of the 24 µm map is given by where this time M ⋆,i,j is fixed and the five k i constants are unknown. M ⋆,i,j is the current stellar mass of the i-th SFH step of the j-th pixel, which is proportional to ψ i,j (proportionality constants listed in Table 3, for unit SFRs). The k i constants are therefore the mass-to-light ratios.
We then sum Equation 35 over a given region S and rewrite it as where f i,24,S ≡ k i M ⋆,i,S /L mod 24,S is the fraction of the 24 µm luminosity from region S arising from the i-th SFH step. M ⋆,i,S = j∈S M ⋆,i,j and L mod 24,S = j∈S L mod 24,j are, respectively, the total stellar mass and the total 24 µm luminosity from region S, associated with the i-th SFH step.
The goodness-of-fit χ 2 can be written as Note the procedure can be applied to any IR broadband image of the galaxy replacing the 24 µm. Also, Equations 35 and 37 are very similar to Equations 13 and 14, but here we have as many observables as pixels in region S and still five unknown fractions. The number of pixels in each region is listed in Table 2. We therefore apply our previously discussed fitting procedure in order to solve for the five fractions, for each region and each IR band. Table 5 lists the derived fractions f i,λ,S for all IR bands and regions. The largest contributor is in bold to guide the eye. The result is very surprising. For M51a the majority of the IR emission is coming from the stellar populations associated with the third SFH step, with ages between 100 Myr and 1 Gyr. For M51b, most of IR emission is mainly coming from the oldest stellar populations, with older than 5 Gyr. The outskirts region has most of the IR coming from the populations younger than 10 Myr, but comparable fraction arises from populations of ages between 100 Myr and 1 Gyr.
Thanks to the diversity of physical properties around M51, many of them influenced by the interaction between M51a and M51b, we observe different mixes of young and old stars, subjected to different attenuations. This diversity of IR fractions presented in Table 5 shows that even though the young stars do heat dust much more efficiently per stellar mass, a much larger population of old stars may compete (as in M51a) and even dominate the IR emission (as in M51b).
In reality, we are able to derive the age-dependent light-to-mass ratios through this procedure for each SFH component. The last columns of Table 5 shows, for each region, the derived IR luminosity at any given band over the stellar mass, separately for each SFH step component. In Figure 13, we display the IR luminosities over current stellar mass associated to each of the five SFH steps. The figure shows the result of decomposing the IR bands using all pixels (as in the last rows of Table 5). This way, we can actually see that dust in different environments, heated by stars of the same age range, generally have the same luminosity per stellar mass at any region. Some of these k i values, however, are changing from region to region and sometimes are significantly different in the outskirt region. We plan to investigate this variation in more detail in a future paper, with a larger sample of spatially-resolved galaxies.
We find old stars do heat dust much less efficiently than young stars (per stellar mass), as shown in Table 5 and Figure 13. Even though the ratio of IR luminosity over stellar mass does drop consistently with age, old stars may outshine younger stars at any of the IR wavelengths studied here, including the TIR luminosity. In fact they do in M51b.
6. SUMMARY We presented a new broadband SED fitting procedure to efficiently derive star-formation histories for a large number of SEDs. We used the Whirlpool system, M51 (composed of the galaxies M51a and M51b) as a initial testbed for our code, Lightning, taking advantage of the fact that extensive data and previous studies are available for M51. For this work, we made use of FUVto-FIR data, including 18 broadband images of M51 from GALEX, SDSS, 2MASS, Spitzer, and Herschel.
We developed a method based solely on Spitzer bands to detect foreground stars in the Galaxy, masked the brightest foreground stars, and replaced these pixels with the local background in all images from the FUV to Spitzer 4.5 µm. We then convolved all images to a common spatial resolution of 25 ′′ FWHM and registered them to a common astrometric frame with 10 ′′ pixels. At the distance of 8.58 Mpc, the pixel size corresponds to 416 pc. The FIR images (24, 70, 160, and 250 µm) were combined to produce a TIR luminosity map, based on Galametz et al. (2013) calibration of Draine & Li (2007) models.
We employed Lightning to individually fit 2043 pixels, each with 12 photometric bands from GALEX FUV to Spitzer 4.5 um and a TIR luminosity. We model the SFH of each pixel as a function of time with five steps: 0-10 Myr, 10-100 Myr, 0.1,1 Gyr, 1-5 Gyr, and 5-13.6 Gyr. The total attenuated power is constrained by the observed TIR and we model the extinction curve with three parameters. The code has a 3D parameter grid for the extinction parameters and maximizes likelihood in each position of parameter space by an efficient inversion procedure that uniquely determines all five SFH step intensities at once. We have derived 400 extra simulated maps to properly derive uncertainties for all quantities. We do not employ any dust emission modeling in this version of the code, but we plan to add this feature in a later version. This way we have produced reliable SFH and extinction maps for M51, which we make publicly available to the astronomy community at https://lehmer.uark.edu/. We have compared the TIR luminosities, stellar masses, and overall SFHs on a pixel-by-pixel basis with the publicly available code CIGALE and found remarkable agreement.
We also tested the hybrid SFR tracer combining UV and IR wavelengths and observed similar trends to the ones previously found by Eufrasio et al. (2014) and M51a 70 33.9 +12.3 convert a given IR band intensity into the attenuated UV luminosity was found to increase for increasing sSFR100, i.e., younger regions show more efficient UV absorption per unit IR emission. This seems to imply the older stellar population, not associated with the recent SFR of the last 100 Myr, is still able to considerably heat dust, many times outshining younger regions in the IR.
With this in mind, we employed a procedure to determine the fraction of luminosity associated with each SFH step at each FIR wavelength. Interestingly, we found significant emission from all five SFH steps, at all FIR wavelengths. In fact, dust heated by the oldest stellar population (5-13.6 Gyr) outshines each of the other four in M51b. In M51a, the IR emission mainly arises from dust heated by the population most likely produced with the interaction (third SFH step, 0.1-1 Gyr). And the outskirts of the interacting system emits almost equal amounts of IR emission associated with the youngest (0-10 Myr) SFH step and well as the third (0.1-1 Gyr). From this decomposition, we were also able to derive the ratio of IR luminosity at a given band and stellar mass associated with each SFH step, for the four large regions investigated here. The general trend is a decline of dust luminosity per stellar mass for all the wavelengths consid-ered here (24, 70, 160, and 250 µm and the TIR). In other words, we find stars of all ages heat dust. Also, younger the stellar population, more efficient this heating. But surprisingly, older stars may outshine the youngest stellar population even in a spiral galaxy, depending on the SFH of the galaxy.
Our work shows how a simple SFH model as steps of time can be useful and opens a promising possibility of decomposing any given map of a galaxy (be it X-rays, IR, or radio) into various age bins to determine its timeevolution. In fact, we have already successfully decomposed the X-ray binary luminosity functions across M51 into the five steps used here to empirically determine its evolution. This will be published under a separately paper (Lehmer et al. 2017, submitted).
We thank the anonymous referee for comments that significantly improved the paper. We gratefully acknowledge support from NASA/ADAP grant NNX13AI48G (B.D.L., R.T.E., A.Z.). R.T.E. and E.D. acknowledge NASA ADAP proposal NNH11ZDA001N. A.Z. acknowledges funding from the European Union's Seventh Framework Programme (FP/2007(FP/ -2013/ERC Grant Agreement n. 617001. Based on observations made with the NASA Galaxy Evolution Explorer. GALEX is oper- Figure 13. Age-dependent light-to-mass for all FIR bands and the TIR associated to all our five SFH steps, as estimated from the whole map of M51. The first SFH step is actually from 0 to 10 Myr, but it is displayed here from 1 to 10 Myr for simplicity. Violet, cyan, green, red, and black respectively represent the 24, 70, 160, 250, and TIR luminosities divided by the stellar mass of each SFH step. Dotted lines show the 16th and 84th percentiles for each quantity. The ratio of IR luminosity over stellar mass drops consistently with age, meaning the old stars do heat dust much less efficiently than young stars (per stellar mass). Nevertheless, old stars may, and in fact they do in M51b, outshine younger stars at any of the IR wavelengths studied here, including the TIR luminosity (see Table 5). ated for NASA by the California Institute of Technology under NASA contract NAS5-98034. This publication makes use of data products from the Two Micron All Sky Survey, which is a joint project of the University of Massachusetts and the Infrared Processing and Analysis Center/California Institute of Technology, funded by the National Aeronautics and Space Administration and the National Science Foundation. This work is based on observations made with the Spitzer Space Telescope, obtained from the NASA/IPAC Infrared Science Archive, both of which are operated by the Jet Propulsion Laboratory, California Institute of Technology under a contract with the National Aeronautics and Space Administration.