Hierarchical Bayesian approach for estimating physical properties in spiral galaxies: Age Maps for M74

One of the fundamental goals of modern Astronomy is to estimate the physical parameters of galaxies from images in different spectral bands. We present a hierarchical Bayesian model for obtaining age maps from images in the \Ha\ line (taken with Taurus Tunable Filter (TTF)), ultraviolet band (far UV or FUV, from GALEX) and infrared bands (24, 70 and 160 microns ($\mu$m), from Spitzer). As shown in S\'anchez-Gil et al. (2011), we present the burst ages for young stellar populations in the nearby and nearly face on galaxy M74. As it is shown in the previous work, the \Ha\ to FUV flux ratio gives a good relative indicator of very recent star formation history (SFH). As a nascent star-forming region evolves, the \Ha\ line emission declines earlier than the UV continuum, leading to a decrease in the \Ha\/FUV ratio. Through a specific star-forming galaxy model (Starburst 99, SB99), we can obtain the corresponding theoretical ratio \Ha\ / FUV to compare with our observed flux ratios, and thus to estimate the ages of the observed regions. Due to the nature of the problem, it is necessary to propose a model of high complexity to take into account the mean uncertainties, and the interrelationship between parameters when the \Ha\ / FUV flux ratio mentioned above is obtained. To address the complexity of the model, we propose a Bayesian hierarchical model, where a joint probability distribution is defined to determine the parameters (age, metallicity, IMF), from the observed data, in this case the observed flux ratios \Ha\ / FUV. The joint distribution of the parameters is described through an i.i.d. (independent and identically distributed random variables), generated through MCMC (Markov Chain Monte Carlo) techniques.


Introduction
The study of the star formation history (SFH) and star formation rate (SFR) in galaxies provide vital information on the evolutionary properties of galaxies and the physical processes which drive that evolution. The variation in star formation (SF) activity across galaxies of different types is influenced by many factors, including gas content, mass and dynamical environment ([Kennicutt (1998)]). But even galaxies with the same morphological type show a diversity in SF and enrichment histories (e.g. [Grebel (2000)]). Some SFR measurements using Hα emission suggest episodic bursts of SF, rather than continuous SF ( [Glazebrook et al. (1999)]). Trends in SFR and SFH across galaxy types allow predictions of galactic evolution with cosmic lookback time. In those galaxies with resolved stars it is possible to determine the SFH from analysis of the color-magnitude diagrams in regions with enough stars to reach sensible results. However individual stars are only resolved in the closest galaxies, so global SF properties of galaxies are mostly obtained by integrated light measurements. The most common diagnostic methods use measurements in the ultra-violet (UV), far-infrared (FIR), or nebular recombination lines ([Kennicutt (1998)]). Of particular interest for this work are measurements in the UV and the Hα recombination line.
One way to deal with the study of the SFH in galaxies is to obtain the age map of the galaxy. In this paper we propose to study the spatial distribution of SFH in galaxies on a pixel by pixel basis, without any a priori hypothesis about the extent and spatial delineation of the star forming regions. Ages are obtained from comparison between Hα and UV emission, so this method is only sensitive to exploring the age map of the youngest stellar population, less than about 10 Myr. This age map should provide a global vision of the current star formation processes taking place in the galactic disk, the maximum scale of coherent star formation, and its relation with other large scale processes such as density waves.
This work is part of a larger study of a sample of nearby galaxies, where star forming regions are spatially resolved, in order to place the relationship between star formation, ultraviolet and Hα emission on a stronger empirical foundation ( [Sánchez-Gil et al. (2011)]; hereinafter Paper I).
The comparison of UV and Hα fluxes is used as a tracer of the recent star formation history. While the infrared images are used to correct for extinction via the TIR to FUV ratio. This flux ratio is defined as a robust and universal tracer of dust extinction ( [Buat & Xu (1996)]; [Buat et al. (1999)]; [Gordon et al. (2000)]). The measured Hα/FUV ratios are compared with the stellar population synthesis (SPS) models to study the dynamics of stellar populations of different ages within the galaxies; we use Starburst99, SB99 ( [Leitherer et al. (1999)]; [Vázquez & Leitherer (2005)]) 1 . Model flux ratios are compared to the observed values for each pixel in the images, resulting in two dimensional spatial age maps.
We address the problem of deriving the galaxy age map from Hα/FUV flux ratio images by establishing a probabilistic framework which explain the relationship of the random variables involved in the problem. This relationship will be formulated in terms of a joint probability distribution given the observations and their uncertainties. We define this joint probability distribution in terms of a hierarchical Bayesian model (HBM) and use Markov Change Monte Carlo (MCMC) to describe it (Nested Sampling) In order to set up the HBM we need to generate a model flux ratio 2 from SB99 model and compare it with the observations by using a likelihood. Finally, we explore the parameters of interest (age) marginalizing nuisance parameters in the posterior distribution defined by HBM model.

Observations
The Hα image was taken with the Taurus Tunable Filter (TTF; [Bland-Hawthorn & Jones (1997)]) on the William Herschel Telescope (WHT) on 1999 March 4−6. Conditions were photometric with stable seeing of 1.0 arcsec. TTF was tuned to a bandpass of width ∆λ = 20 was used to remove transmissions from all but a single interference order. The pixel scale is 0.56 arcsec. More details on the data and data reduction can be found in Paper I.
The UV image comes from the Nearby Galaxies Survey of the Galaxy Evolution Explorer mission (NGS survey, GALEX, [Martin et al. (2005)]). This survey contains well-resolved imaging (1.5 arcsec pix) of 296 and 433 nearby galaxies for GR2/GR3 and GR4 releases, respectively, in two passbands: a narrower far-ultraviolet band (FUV; λ eff /∆λ = 1516/268 • A), 1 http://www.stsci.edu/science/starburst99/ 2 We prefer use the term model instead than true because its value depends on the assumption that SB99 is an appropriate model Figure 1: Examples of processed frames in Hα (left), far ultraviolet (FUV; centre), and total infrared (TIR; right) for the galaxy M74. The images have been resampled to have identical size, orientation, and pixel scale (1.5"/pix). and a broader near-ultraviolet band (NUV; λ eff = 2267/732 • A). Archival Spitzer 3 images, at 24, 70, and 160µ, were used to provide additional estimates of extinction. They were resampled to a common 1.5 arcsec/px scale and combined into an image of TIR flux, The extinction within our galaxy is corrected using the [Schlegel et al. (1998)] 4 dust maps, which measure the Galactic extinction in all directions. The internal reddening is corrected using a straight relation between the A F UV extinction factor and the F T IR /F F UV flux ratio: Buat et al. (2005)]). The F T IR /F F UV ratio is a robust and universal tracer of dust extinction, almost independent of the dust/stars geometry and of the dust properties, provided that the galaxies are actively forming stars ( [Buat & Xu (1996)]; [Buat et al. (1999)]; [Gordon et al. (2000)]). The A(Hα) extinction factor was derived through the relation A F UV = 1.4A(Hα) ( [Boissier et al. (2005)]).
The resampled and aligned Hα, FUV and TIR images are shown in Fig.1, center panel.

Methodology
In this section, we address the problem of deriving the galaxy age map (GAM) from F Hα /f F UV flux ratio images by establishing a probabilistic framework which explain the relationship between random variables involved in the problem. This relationship will be formulated in terms of a joint probability distribution given the observations and their uncertainties using a hierarchical Bayesian model (HBM) [Gelman et al. (2003)]. Specifically, we want to describe the probability distribution p(age|r obs ), where r obs is the observed Hα/FUV flux ratio. Because the use of SB99 we extend the previous distribution to where the first component of θ is related to the age, the second is related to the initial mass function (IMF) and third component is related to metallicity. Furthermore, by marginalizing we can derive the age distribution according to p(age|r obs ) = p(θ 1 |r obs ) = p(θ|r obs ) dθ 2 dθ 3 .
Figure 2: Hierarchical model using plate notation. Hyperparameters and parameters are written with φ k and θ k respectively. Observations and model ratios are inside the plate. The letter N notes the sample size per pixel, in our case, we only observed one pixel at time, then N = 1. Squares are fixed quantities and circles random variables. If the node is shaded then the variable is known (observed).
The joint probability distribution p(θ|r obs ) can be rewritten by using the Bayes theorem according to where p(θ) is the prior distribution and p(r obs |θ) is the likelihood function and p(r obs ) is a normalization constant according to p(r obs ) = p(r obs |θ)p(θ) dθ (4)

The prior distribution
The HBM is plot in figure 2 using plate notation. Prior distributions are represented with solid line linking hyperparameters φ i and θ i with i ∈ {1, 2, 3}. These hypermaramters are konwn and fixed in order to draw the parameters θ according to different prior distributions. Priors for θ were set as (5) where Cat(n) is a categorical distribution with n categories. For the IMF, the categories are: first, a Salpeter law with α = 2.35 and M up = 100M ⊙ ; second, a truncated Salpeter law with α = 2.35 and M up = 30M ⊙ ; third, and a Miller-Scalo law with α = 3.3 and M up = 100M ⊙ . For metallicity, categories are Z = 0.04, 0.02 (solar, Z ⊙ ), 0.008, 0.004 and 0.001. Finally the prior probability distribution for age was set to the uniform distribution in the interval 0-10 Myr.

The likelihood
The likelihood function is represented by lines linking the variables θ i with r obs in Fig. 2. In order to obtain model flux ratio 5 r mod = F Hα /f F UV from SB99 model, we need to set up SB99 with parameter θ = (θ 1 , θ 2 , θ 3 ), but let us notice that the SB99 model only permits a finite set of values for θ 2 and θ 3 , see Fig. 3. This situation means we cannot work with a continuous parameter space in three dimensions but a grid of parameters {θ} grid . Because this finite combination of parameters in SB99, we can only derive a finite number of model flux ratios from them. Therefore, in order to obtain more model flux ratios different from those on the grid, we use an Artificial Neural Network (ANN) which will interpolate the grid of parameters. A deep information about ANN can be found in [Haykin, S. (1999)]. Specifically the ANN is a multilayer perceptron with four layers: two hidden layers with 20 and 50 nodes, the input layer for the parameters θ and the output layer for the Hα and FUV flux. The ANN was trained with 70 % and validate with 30 % of the grid points. The uncertainties due this interpolation were considered negligible. Dashed lines in Fig. 2 represent the ANN.

Independence between pixels
For explanation purposes this section assume independence between the fluxes of different pixels. Also the results are established for an unique pixel with r obs > 0. Once the SB99 have generated the model fluxes we assume that the observed flux from each pixel has a 5 We prefer use the term model instead than true because its value depends on the assumption that SB99 is an appropriate model normal distribution where F Hα,obs , f F UV,obs , F Hα,mod and f F UV,mod are the observed and model Hα and FUV fluxes for the selected pixel respectively. The quantitiesσ Hα andσ Hα,obs are the flux uncertainties due to the instrumental, reduction and correction processes. Therefore, r obs = F Hα,obs /f F UV,obs is the ratio of two correlated normal random variables, which exact distribution is given by ( [Hinkley, D. V. (1969) where, where ρ is the correlation between Hα and FUV, and Φ is the cumulative density distribution of the standard normal, The F Hα and f F UV are correlated, in fact in normal star forming galaxies it is expected a constant flux ratio FHα fF U V ≃ 11 • A ( [Meurer et al. (1999)], [Iglesias-Páramo et al. (2004)]). In order to compare model and observed flux ratio, pixel-by-pixel, we consider the likelihood according to p(r obs |θ) = ψ(r obs |F Hα,mod , f F UV,mod ,σ Hα ,σ F UV , ρ) where F Hα,mod and f F UV,mod are obtained from θ through the ANN. The theoretical correlation, ρ, between F Hα and f F UV is estimated from all the pixels of the images .

The posterior distribution
We characterize the posterior distribution in equation 1 by extracting independent samples using the Nested Sampling Algorithm (NSA) [Skilling, J. (2006)]. The aim of the NSA is the estimation of p(r obs ) defined in 4, and as by-product an independent sample of the posterior distribution is obtained. The NSA is based on the relationship between the likelihood p(r obs |θ) and the prior volume X(λ) defined by Figure 4: Sample of the marginal probability distribution for the age parameter, for the pixel #34872 of the image. The blue and red dashed lines represent the median and the mean respectively. The solid black line shows the age with the maximum likelihood.
which means the bulk of prior distribution contained within an iso-contour of the likelihood.
Finally we obtain the posterior distribution of our parameter of interest (age) by marginalizing nuisance parameters in the posterior distribution as we set in equation 2. Figure 4 shows the sample of the marginalized the age distribution for a certain pixel of the image, as an example of the resulting posterior distribution.

Results
Figures 5 and 6 show some age maps for M74. In the Figure 5, we can compare between the age maps for M74 obtained in Paper I, by applying a robust, only based on astrophysical models. Whereas on the left, it is shown the age map when applying the Bayesian approach. The main difference, and more remarkable is that the first one is a discrete age map mode. However, with the Bayesian inference we can determine a single age to each pixel. Resulting therefore in more rich age patterns ans structures, which could not been observed with the latter.
To have an idea of the uncertainty in the age estimation, it is included the Figure 6, where it is represented the 50% central distribution of the age. On the left it is shown the Percentile 25, P 25 , of the age distribution. The central age map corresponds to the median ages, P 50 , and on the left the P 75 . With this comparison we get as a confidence interval, checking out the robust of this methodology for age estimation.
Finally, with respect to the physical conclusions we can observe that these are basically the same that were obtained in Paper I. The age map shows an age gradient from the inner to the outer parts of the galaxy, from very recent to less recent episodes of star formation, in agreement with previous authors ( [ Cornett et al. (1994)], [Lelièvre & Roy (2000)]).
Specifically, we find that the Hα luminosity decrease in radius is more pronounced in the inner 5 to 6 kpc, while the UV luminosity shows a shallower rate of change. Consequently, the Hα/FUV ratio decreases with radius indicating an age increases in the outward direction.
On more localised scales, the short arm that opens S-SW at 4 kpc shows a clear age gradient across it. However, the outer longer arm, that runs SE-S at ∼ 5 to 10 kpc, showed a less marked age gradient across its width in Paper I. But improving the "resolution" of our age maps with the bayesian technique, it is more remarkable the age gradient across the outer longer arm as well. If the age gradients across spiral arms are a direct product of the spiral density wave, then the dilution of the gradient in this southern arm may be related to the weakness of the density wave or the approach to corotation. As discussed by [Efremov (2009)], the presence of a shock produced by the spiral density wave, (made visible by dust lanes along the spiral arms), is incompatible with the creation of star forming complexes, because of the absence of visible dust lanes in this arm, despite a chain of complexes along its length. The thickness of the longer arm is basically dominated by a single age range, and the youngest population located in the inner edge of the arm maps the location of the chain of complexes observed in this arm.
Upcoming work is applying this method to the rest of the sample of 6 galaxies from the original work, Paper I. To analyze with more detail now the age patterns. Next, we are studying also the more complex case and model of dependence between pixels of the image. To develop a parallelization of the algorithm, high perfomance computing to emcee, and / or HMC (programming in Python and R) to get better computing times. And finally, to get a wider, new sample of galaxies, to apply the method and determine not only their age maps, but also to study possible patterns/gradients in the metallicity maps.

A Empirical Model: Starbust 99
Comparison of UV and Hα fluxes provides a good indication of the recent star formation history. As a young star forming region evolves, the Hα emission line declines before the UV emission does, leading to a decrease in F Hα /F F UV , which is very sensitive to the age of the youngest stellar population. Images at these wavelengths can be examined to track the recent star formation history and distribution, and to analyze how the SFH correlate with the galactic structure and the galactic environment.
Model Hα/FUV flux ratios were generated from Starburst99, as a function of age, to be compared with the measured flux ratios at each pixel in the image. The full models cover the age range 10 6 to 10 9 yr in steps of 1 Myr with spectral energy distributions (SEDs) spanning 100 • A to 1 µm in wavelength. However, we restricted the ages from 1 to 10 Myr as we are only interested in the youngest stellar populations responsible for the Hα and UV output.
The code can be run with two different star formation modes: an instantaneous burst or continuous. It has three alternatives in the stellar initial mass function, (IMF): a Salpeter law with α = 2.35 and M up = 100M ⊙ ; a truncated Salpeter law with α = 2.35 and M up = 30M ⊙ ; and a Miller-Scalo law with α = 3.3 and M up = 100M ⊙ . And five metallicities are available for each IMF: Z = 0.04, 0.02 (solar, Z ⊙ ), 0.008, 0.004 and 0.001.
We take as our reference model a Salpeter IMF, the interpolated metallicity value Z=0.02 (according to [Miller & Hodge (1996)], M74 has 12+log(O/H)∼8.15), and an instantaneous star formation law, since it is more sensitive to the age variation for younger regions (see paper I).
The age map is contoured in four age bins from the set 0 -4 Myr, 4 -6 Myr, 6 -9 Myr and older than 9 Myr, (this binning corresponds to the internal precision of the flux ratio, Paper I) see right panel of Fig.1 . We find in the resulting age map that the youngest star forming population (0 -4 Myr) surrounds an older one (>9 Myr) located at the centre of the complex, and a clear center to outer rim age gradient.
The F Hα /F F UV depends on the model parameters, such as the IMF, the metallicity, the SFH, etc. To assess the suitability and robustness of our reference model we calculated the effect of changing various model inputs on the SB99 F Hα /F F UV ratios (cf. Fig.2 of Paper I). The photometric uncertainty on the F Hα /F F UV ratio also affects to the age calibration. Aside the model and photometric uncertainties, the age-dating technique is subject to a number of potential sources of systematic error as: (i) the lowest limit on cluster mass allowable for our assumptions on ionizing flux, (ii) the effect of changing the spatial bin size, and (iii) the effect of changes to the metallicity and IMF assumptions in the model. We analyze the reliability and robustness of our age-dating technique taking into account all these uncertainties. See Paper I for further details.
The total uncertainty of the extinction corrected Hα/FUV flux ratio is estimated, by  [Buat et al. (2005)]). These uncertainties are computed pixel by pixel in the image. An average of the relative errors of the respective images are E Hα ≃ 5%, E F UV ≃ 25%, and E T IR ≃ 6%, resulting in an overall uncertainty in the F Hα /F F UV flux ratio of ∼28%. Now, we calculate the underlying stellar mass per pixel and compare with the values of the lower mass limit M min , given by [Cerviño et al. (2003)]. The pixel stellar mass is a lower limit derived by comparing the linearly scaled extinction-corrected observed L F UV with the highest expected value from SB99, for a stellar population mass of 10 6 M ⊙ at the youngest cluster age (1 Myr). In our map, these lowest mass values range logM min = 3.2 − 5.2, with an average of 3.84 ± 0.33, around the M min estimated by [Cerviño et al. (2003)]. Binning the IC 2574 image by 3 × 3 puts the spatial sampling of IC 2574 well above the threshold of M min (cf. Fig.5 of Paper I ), to remove any cause of bias by IMF fluctuations.