. : http://www.eso.org/~rsiebenm/agn_models/Johnson13.pdf
: Fri Jul 17 12:04:52 2015
: Sun Apr 10 01:39:22 2016
MNRAS 436, 25352549 (2013)


Advance Access publication 2013 October 9


Spectral energy distribution Analysis Through Markov Chains

S. P. Johnson,1 < G. W. Wilson,1 Y. Tang1 and K. S. Scott2
1 2

Department of Astronomy, University of Massachusetts, Amherst, MA 01003, USA North American ALMA Science Center, National Radio Astronomy Observatory, Charlottesville, VA 22903, USA

Accepted 2013 September 14. Received 2013 September 12; in original form 2013 August 5


We present the general purpose spectral energy distribution (SED) fitting tool SED Analysis Through Markov Chains (SATMC). Utilizing Monte Carlo Markov Chain (MCMC) algorithms, SATMC fits an observed SED to SED templates or models of the user's choice to infer intrinsic parameters, generate confidence levels and produce the posterior parameter distribution. Here, we describe the key features of SATMC from the underlying MCMC engine to specific features for handling SED fitting. We detail several test cases of SATMC, comparing results obtained from traditional least-squares methods, which highlight its accuracy, robustness and wide range of possible applications. We also present a sample of submillimetre galaxies (SMGs) that have been fitted using the SED synthesis routine GRASIL as input. In general, these SMGs are shown to occupy a large volume of parameter space, particularly in regards to their star formation rates which range from 30 to 3000 M yr-1 and stellar masses which range from 1010 to 1012 M . Taking advantage of the Bayesian formalism inherent to SATMC, we also show how the fitting results may change under different parametrizations (i.e. different initial mass functions) and through additional or improved photometry, the latter being crucial to the study of high-redshift galaxies. Key words: methods: statistical techniques: photometric galaxies: fundamental parameters galaxies: high-redshift submillimetre: galaxies.

Downloaded from http://mnras.oxfordjournals.org/ at European Southern Observatory on July 17, 2015

1 I NTR O DUCTION The complete panchromatic spectral energy distribution (SED) of a source encodes a wealth of information concerning its age, mass, metallicity, dust/gas content, star formation rate (SFR), star formation history (SFH) and more. Various emission mechanisms account for the apparent shape and typical features found in SEDs. For example, dust reprocessing results in attenuation of optical/ultraviolet (UV) photons produced from stellar populations whose energy is then re-emitted at infrared (IR) wavelengths while specific species of ions and molecules produce emission/absorption features across the electromagnetic spectrum. Unfortunately, our ability to extract information from observations is hampered both on the observational and theoretical sides. As observers, we try to make due with either coarse, broad-band sampling of a source's SED or with highresolution sampling of a small portion of the SED through spectroscopy. Theorists, on the other hand, must make difficult decisions regarding which physics, spatial scales and evolutionary histories to include in the creation of their SED libraries or synthesis models. Within the literature, SED models can be subcategorized into two main types. Empirical models are derived for particular classifications based on a subset of similar sources (e.g. Arp 220 and
E-mail: spjohnso@astro.umass.edu

M82 for starburst galaxies, the quasar mean template, etc.). These models offer the simplest approximation of a source's SED and are generally preferred when only sparse photometry is available or for coarse estimates of basic properties (e.g. luminosity, SFR, colours). The underlying assumption behind empirical models, namely that all sources of that `type' have the same SED, is difficult to verify and so their use to probe all but the grossest properties is limited. Theoretical models are constructed from sets of physical and radiative processes believed to be the dominant contributors to the emission of a source. Functionally, these again are divided into two classes. Pre-computed template libraries for particular source types are the most widely used (e.g. Calzetti et al. 2000; Efstathiou, Rowan-Robinson & Siebenmorgen 2000; Siebenmor gen et al. 2004; Siebenmorgen & Krugel 2007; Gawiser 2009; Michalowski, Watson & Hjorth 2010 and references therein). Since they are pre-computed, these libraries allow the rapid exploration of a pre-defined parameter space at the expense of being limited to the resolution and scope of the parameter space provided by the authors. For more generalized applications, SED synthesis packages are now available that offer a wider set of input parameters and the exploration of a continuous parameter space [e.g. stellar population synthesis codes such as GALAXEV (Bruzual & Charlot 2003) and panchromatic galaxy synthesis codes like GRASIL (Silva et al. 1998) and CIGALE (Noll et al. 2009)]. Of course, one may use these

C 2013 The Authors Published by Oxford University Press on behalf of the Royal Astronomical Society


S. P. Johnson et al.
where the product runs over the M observed data points D whose individual variances are given by i2 and the function f ( d | x ) represents the model observations d for the set of x parameters. Note, however, that this prescription of L assumes that the observations have Gaussian distributed errors, an assumption we will return to later in Section 3.1. In the following sections, we detail specific MCMC features that define the basic operation of SATMC.

packages to construct one's own template libraries (e.g. Michalowski et al. 2010) for specific applications as well. In both cases, as the generality of the underlying physics increases to provide relevance to a wider class of sources, so does the number of free parameters in the models. The unavoidable existence of correlations in these parameters insists that simply finding a best-fitting parametrization is no longer sufficient. Rather, we now require tools that properly reveal the complexities and correlations in the adopted model parameter space. With this in mind, here we present the general purpose MCMCbased SED fitting code SED Analysis Through Markov Chains or 1 SATMC. Monte Carlo Markov Chain (MCMC) techniques are a set of methods based in the Bayesian formalism which enable efficient sampling of multidimensional parameter spaces in order to construct a distribution proportional to the probability density distribution of the input parameters, known as the posterior parameter distribution or simply the posterior. The posterior identifies the nuances in parameter space, including any possible correlations, making MCMCs particularly useful for exploring SED models and libraries with high dimensionality. MCMC-based SED fitting codes are relatively new but have been growing in popularity (e.g. Sajina et al. 2006; Acquaviva, Gawiser & Guaita 2011; Serra et al. 2011; Pirzkal et al. 2012). In addition to deriving the posterior for bestfitting parameter and confidence level estimation, SATMC includes many features to aid in improving performance and allows users to easily incorporate additional knowledge and constraints on parameter space in the form of priors. Additionally, SATMC versions exist in both IDL and PYTHON and both are modular and straightforward to use in any wavelength regime and for any class of sources. This paper is organized as follows. We start by detailing the MCMC algorithm and the basic process for MCMC-based SED fitting. We then provide case examples displaying the versatility and accuracy of SATMC compared to standard least-squares methods. As part of this demonstration, we also present a set of SEDs derived from SATMC when used in conjunction with the SED synthesis routine GRASIL to a sample of X-ray selected starburst galaxies. These fits highlight the key advantages obtained from the MCMC-based approach and their agreement with similar classes of composite SEDs.

2.1 MCMC acceptance and convergence Within the literature, there are a variety of sampling algorithms one may use to construct an MCMC (e.g. MetropolisHastings or Gibbs sampling; Metroplois et al. 1953; Hastings 1970; Geman & Geman 1984; Geyer 1992; Chib & Jeliazkov 2001; Verde et al. 2003; Mackay 2003). For SATMC, we employ the MetropolisHastings algorithm which works as follows: (i) Generate a proposal distribution q ( x i |xi -1 ) from which the candidate steps x i will be drawn. (ii) Calculate the acceptance probability according to the likeli(i )(i - hood ratio ( = min(1, PPxx-|1DDqqxx|ixi1 |1x)i ) )). (i | )( - (iii) Draw a uniformly distributed random number u from 0 to 1 and accept the step if u < , reject otherwise. (iv) Repeat for the next step. Though the exact choice of q ( x i |xi -1 ) is arbitrary, it is common to adopt an nD -dimensional multivariate-normal distribution 2 N (, ) with = x I where I is the identity matrix, x is the variance of each parameter in x and N (, ) | |- 2 e-
1 1 2

Downloaded from http://mnras.oxfordjournals.org/ at European Southern Observatory on July 17, 2015

( x -)


( x -)





The primary motivation for SATMC is to provide a means for efficient sampling of a parameter space with nD free parameters in order to derive parameter estimates and their associated confidence intervals. This is accomplished by sampling an nD -dimensional surface proportional to the probability density function of the parameters given the data P ( x | D), also referred to as the posterior parameter distribution. We determine the posterior through Bayes theorem P ( x | D ) P ( D | x )P ( x ), (1)

for generating proposed steps (e.g. Gelman, Roberts & Gilks 1995; Roberts, Gelman & Gilks 1997; Roberts & Rosenthal 2001; Atchade & Rosenthal 2005). Following these works, it is found that the distribution N (, ) should be tuned for optimal performance between the time necessary to reach a stable solution (referred to as the burn-in period) and sampling of the posterior, which may be achieved by adjusting N (, ) until the resulting chains have acceptance rates of 23 per cent in the limit of large dimensionality. We therefore base our proposal distribution around an adaptive covariance matrix (similar to Acquaviva et al. 2011) such that samples are drawn from the multivariate-normal distribution now given by N ( x i -1 , ). To obtain an acceptance rate of 23 per cent, 2 we start by initializing the covariance matrix for each chain as x I, where x is set proportional to the input parameter ranges. After a period of steps, we then calculate directly from the chain over the previous interval, i.e. = E [( x - E [ x ])( x - E ( x ))T ], (4)

where P ( x ) represents our current knowledge of the parameters or priors and P ( D| x ) is the probability of the data given the model parameters, often referred to as the likelihood L. We shall define the general form of the likelihood according to

P ( D| x ) = L =
i -1


(Di - f (di | x )) 2i2






where E [ x ] is the expectation value or weighted average of x . Following each period of steps, we compute the acceptance rate and scale if the acceptance rate is too high (>26 per cent) or too low (<20 per cent). The covariance matrix is continuously updated until the target 23 per cent acceptance is reached, which allows to take on the shape of the underlying posterior to readily identify and account for possible correlations in the parameters. Once the target acceptance rate has been reached, we then check for chain convergence to determine when the burn-in period is complete. In the MCMC literature, there are many approaches to determine convergence (e.g. Gelman & Rubin 1992; Raftery & Lewis 1992; Gilks, Richardson & Spiegelhalter 1996); we opt for

SED Analysis Through Markov Chains
the Geweke diagnostic (Geweke 1992), which compares the average and variance of samples obtained in the first 10 per cent and last 50 per cent of a chain segment. If the two averages are equal (within the tolerance set by their variances), then the chain is deemed to be stationary and convergence is complete. We check both the acceptance rate and convergence, verifying the acceptance rate before checking for convergence, over a default period of 1000 steps until both have been satisfied. If necessary, the covariance matrix is modified to maintain the target acceptance rate. Once both criteria are fulfilled, burn-in is completed and the chain(s) are set to continue to provide sampling of the posterior. 2.2 Parallel tempering When constructing the posterior distribution, we must be careful to sample all of the nuances of parameter space as the posterior distribution is not guaranteed to be a smooth or even continuous function of the free parameters. As with all fitting approaches, the existence of local maxima in the posterior requires that either we known a priori the general location of the global maximum or we implement a sampling technique capable of increasing the probability of finding it. For example, one could choose to initialize a large number of chains that cover random locations throughout parameter space. While this approach is nearly guaranteed to find the global maximum, it is also extremely inefficient. SATMC utilizes a technique known as `Parallel Tempering' (see review by Earl & Deem 2005) which, in analogy to simulated annealing, uses several chains, each with progressive modifications to likelihood space parametrized as a statistical `temperature', to search parameter space and exchange information about the posterior at each chain's location. For a given chain at temperature T, the likelihood is `flattened' according to LT = L1/(1+T) . The tempered chains thus distort likelihood space, exchanging sensitivity to the details within posterior for the general shape, such that chains with higher temperatures will accept more steps, and thus sample larger regions of parameter space. By coupling these tempered chains, we allow `colder' chains to access areas of parameter space they may have otherwise been unable to reach. The process for handling tempered chains works as follows: (i) Temperatures are assigned to chains in a progressive manner with one chain designated the fiducial `cold' chain with T = 0 (e.g. T = [0, 10, 100, 1000]). (ii) Chains are allowed to progress for a set number of steps (SATMC uses three iterations of acceptance/convergence or 3000 steps by default). (iii) A `swap' of chain state information is proposed using the Metropolis acceptance algorithm = min(1,
1/(1+Tj ) (P ( x i | D)) (P ( x j | D))1/(1+Ti ) 1/(1+Tj ) (P ( x i | D))1/(1+Ti ) (P ( x j | D))


Downloaded from http://mnras.oxfordjournals.org/ at European Southern Observatory on July 17, 2015

Figure 1. Simple visualization of transitions between tempered chains. Transitions may occur either between the cold and tempered chains (top) or between adjacent pairs (bottom). The former allows for rapid sampling of parameter space to quickly reach global maxima while the latter improves sampling around a single location to avoid local maxima in the construction the posterior parameter distribution. In this example, chain transitions are proposed after every 100 steps; in practice, these transitions are proposed every 3000 steps.


(iv) The new set of chains is then allowed to run until the next swap of chain state information. This process of coupling individual chains with a Metropolis Hastings acceptance algorithm was initially proposed by Geyer (1991) and is referred to as a Metropolis Coupled Monte Carlo Markov Chain or MC3 . When applying temperatures to individual chains, the sampling of each chain will vary due to the distortions of likelihood space; i.e. a chain at a higher temperature will accept more steps than a chain with an identical proposal distribution at lower temperature. To maintain a reasonable sampling of the tempered likelihoods, each chain is treated individually for the purposes of acceptance and convergence testing. Note, however, that updating the tempered chains occurs only with the acceptance and

convergence testing of the cold chain. This ensures that we obtain a proper sampling of the posterior without being influenced by the otherwise distorted likelihood space. In parallel tempering techniques, one can specify whether the potential swaps occur only between the cold and any tempered chain (0, j) or between adjacent chains (i, i 1) (see Fig. 1). The former allows for rapid sampling and mixing of the chains to determine the global maximum and is thus used during the burn-in period. The latter passes state information down through the tempered chains so that the cold chain may access a more representative region of parameter space; this method is used after burn-in to properly sample parameter space around the maximum likelihood. In this case, chains progress for a set length (500 steps by default) whereafter the chain transitions are proposed. Should a swap be made with the cold chain after burn-in is complete, we re-compute the acceptance rate and convergence as outlined in the previous section to verify proper sampling. Due to the modified acceptance rate and misshapen likelihood space of the tempered chains, it is generally undesirable to use them in reconstructing the posterior parameter distribution. If there have been no swaps to the cold chain after 10 iterations, the temperatures of all chains are set to 0 and the MCMC is allowed to continue until


S. P. Johnson et al.

sufficient samples of the posterior distribution around the maximum have been obtained. 3


The methods presented in Section 2 are generalized for any MCMCbased algorithm. Here, we detail specific modifications and methods utilized in SATMC for SED fitting. 3.1 Observations and upper limits In order to perform a fit, SATMC requires at least two sets of information: (1) a file containing the M observations including their wavelengths (i ) or frequencies ( i ), observed fluxes (fobs,i ) and corresponding uncertainties (obs,i ) and (2) the model libraries. Following from equation (2), we define the likelihood for a given set of observations and models as

Downloaded from http://mnras.oxfordjournals.org/ at European Southern Observatory on July 17, 2015

i -1

exp -

(fobs,i - fmod,i ) 2 2fobs,i



or similarly

ln(L) =
i -1


(fobs,i - fmod,i )2 , 2 2fobs,i

Figure 2. Example likelihood distribution for upper limits. Here, the upper limit has a value of 3 mJy (dotted vertical line). The remaining lines represent the likelihood distribution for the simple step function (solid), one-sided Gaussians assuming the upper limit is at the 5 level with a cut-off at 1.8 mJy (dashed) and 0 mJy (dotdashed) and one-sided Gaussians with the upper limit at the 3 level and cut-offs at 2 mJy (triple dotdashed) and 0 (long dash).


where fmod,i are the model fluxes which for narrow passbands are interpolated model SED values at the observed frequency/wavelength. Alternatively, fmod,i may be calculated from a set of passbands with a given filter response [p( ) as provided by the user] as fmod,i =
max min

SED( )pi ( ) d
max min

pi ( )d



where min and max are taken from the filter files. In cases where there is no detection but only an upper limit, SATMC provides a parametrized probability distribution to allow the user to reflect their confidence in the underlying statistics of the observations. A simple approach for incorporating upper limits is to assume a step function accepting all models that fall below the limit and rejecting all that lie above. While attractive in its simplicity, this method harshly truncates regions of parameter space and improperly weights the upper limit's contribution to L. To partially alleviate these issues, SATMC implements a one-sided Gaussian distribution so that models with fluxes above the upper limit may contribute to the posterior with a small but non-zero probability (see also Feigelson & Nelson 1985; Isobe, Feigelson & Nelson 1986; Sawicki 2012). The one-sided Gaussian is defined by 1 if fmod < fco (8) L 2 exp - (fco -f2mod ) if fmod fco , 2

Bolzonella, Miralles & Pello 2000) though some recent codes have adopted the Bayesian formalism for using priors (e.g. BPZ, Benitez 2000; EAZY, Brammer, van Dokkum & Coppi 2008). The standard approach in photo-z estimation is to apply an SED library (often limited to a few galaxy types) to find the photo-z and then repeat the fitting using the same (or in some cases different) library fixed at the photo-z to estimate source properties. This approach underestimates the true errors on the photo-z and other fitted parameters. SATMC fits the redshift simultaneously with all other parameters and produces a direct determination of the redshiftparameter probability distributions in addition to the marginalized redshift probability distribution P(z).2 We have tested and improved the photo-z estimation with SATMC in collaboration with the CANDELS team. Dahlen et al. (2013) provides a complete analysis of various photometric redshift estimation techniques for samples of CANDLES galaxies in the GOODS-S field (Giavalisco et al. 2004). Out of the 13 participating groups, SATMC was the only MCMC-based code used to generate photo-z's. Since the tests reported in Dahlen et al., we have reduced the outlier fraction (fraction of sources with (zspec - zphot )/(1 + zspec ) > 0.15) from 914 per cent to 38 per cent through modification of our input templates, luminosity priors (e.g. Kodama, Bell & Bower 1999; Benitez 2000) and zero-point photometry corrections (e.g. Dahlen et al. 2010). 3.3 Template libraries and SED synthesis routines As MCMC samplers require a continuous parameter space, SATMC allows the user to incorporate SED synthesis routines (e.g. GALAXEV, GRASIL, CIGALE) to generate SEDs at candidate steps and compute the resulting likelihoods. Empirical templates are often defined with a scalable normalization factor as one of the few (if any) free parameters which remains continuous when used with SATMC. Template libraries, unfortunately, rarely offer a fully continuous parameter space, often mixing sets of continuous (e.g. model normalization) and discretely sampled parameters. To create a pseudo-continuous
2 Cosmology is presently fixed in H0 = 70, = 0.7 and M = 0.3.

where fco represents the cut-off transition from flat L = 1 to a Gaussian `tail' with standard deviation co . Parametrizing the probability distribution in this manner allows for a compromise between the simple step function, also available in SATMC, and overinterpreting the shape of the noise distribution at small fluxes, as demonstrated in Fig. 2. 3.2 Photometric redshift estimation Since redshift is just another parameter in a source's SED, SATMC is capable of providing photometric redshift (photo-z) determinations for sources in the context of the input SED models. Traditionally, photo-z codes implement least-squares methods (e.g. HYPERZ;


and assumes flat

CDM with

SED Analysis Through Markov Chains
space from such template libraries, SATMC computes the likelihoods of models bracketing the current step according to equations (5) and (8) and applies multilinear interpolation to determine the likelihood of the current proposed step. We emphasize that SATMC may be used for any class of SED models (empirical, template library or synthesis routine), so long as the appropriate interface is constructed by the user. One should note, however, that there is a trade-off between the discretely sampled template libraries or empirical templates and continuous parameter space offered by SED synthesis routines (see also Acquaviva et al. 2011; Acquaviva, Gawiser & Guaita 2012 since the run-time of the SED synthesis routines will dominate the MCMC calculation (e.g. 3+ days run-time with GRASIL versus 35 min with a template library). Regardless of which class of SED models one wishes to adopt, empirical and theoretical SED models are commonly derived for a particular physical process and/or over a particular wavelength regime (e.g. Bruzual & Charlot 1993; Efstathiou et al. 2000; Siebenmorgen et al. 2004). To fully reproduce a galaxy's SED, additional components may be required either to complete the wavelength coverage or to include a missing physical process (e.g. adding AGN emission to a star formation template set). SATMC will construct a linear combination of multiple input SED models under the assumption that the underlying physical processes are independent. The MCMC process itself does not change: the combination of two models with N1and N2 free parameters, respectively, is viewed as a single model with N1 + N2 parameters when calculating likelihoods and taking potential steps. 3.4 Inclusion of priors An added feature of SATMC is the ability to include additional information to provide additional weights and constraints to likelihood space. In the Bayesian formalism, this extra information forms the priors of equation (1). For our implementation, we expand the definition of priors from the traditional Bayesian definition to include options for limiting and inherently correlating parameter space. This was deemed necessary for circumstances of fitting multiple template libraries where parameters from each model have the same physical interpretation and thus are not independent (e.g. AV from one model library and optical depth in another) and cases where additional information not available to the fits is available (e.g. restricting the age of a galaxy at a given redshift). 3.5 Application of the posterior A final feature of SATMC lies in the determination of the posterior parameter distribution. As we store the likelihood and location in parameter space for each step, it becomes a simple task to construct parameter confidence intervals and even parameterparameter confidence contours to examine relative parameter degeneracies (see Fig. 3). Unfortunately, in order to visualize an nD -dimensional parameter space, we must project or `marginalize' parameter space into a one- or two-dimensional form. When marginalizing sets of parameters, the true shape of the posterior will be distorted which may not reveal correlations in higher dimensions. This also leads to a simplification when reporting the confidence intervals on individual parameters as traditional terms such as `1 ' imply Gaussianity in the posterior which is unlikely to exist. Instead, one-dimensional confidence levels are produced from the parameter range where 68 per cent of all accepted steps are contained, marginalized over all other free parameters. Throughout the text, `errors' quoted when derived from SATMC refer to these marginalized parameter ranges. 4 TESTING OF THE M CMC A LGORITHM


With the SATMC algorithm as outlined in Sections 2 and 3, we now set out to verify the fitting results by analysing sets of well-known sources and template libraries. Though the examples provided here are for galaxy SED modelling, SATMC makes no distinction between source types and may just as easily be applied to Galactic sources, localized regions within a particular galaxy or spectroscopy of individual sources. We begin with the well-known galaxies Arp 220 and M82 and the starburst SED library of Siebenmorgen & Krugel (2007, hereafter SK07) which has already been shown to provide reasonable fits to Arp 220 and M82. The SED library consists of over 7000 templates with emission from a starburst parametrized by its total luminosity Ltot , nuclear radius R, visual extinction AV ,ratio of luminosity from OB stars to the total luminosity and hotspot dust density n. The observed SEDs for Arp 220 and M82 were constructed from data available on the NASA/IPAC Extragalactic Database (NED)3 over the 11500 m wavelength range. To ensure a meaningful comparison to SK07, we utilize the same photometric data for M82 and Arp 220 (including multiple aperture JHK photometry for M82) such that the only difference lies in the fitting method. Table 1 provides the best-fitting models as obtained by SK07 and SATMC with the models and parameterparameter confidence contours shown in Fig. 3. Note that while SATMC will sample all of parameter space within that defined by the input template library, it is not possible to extrapolate likelihood information beyond the limits of the templates. This effect is responsible for the apparent truncation of parameter space seen in Fig. 3. Despite the irregular, non-uniform parameter sampling of the SK07 templates, SATMC closely recovers similar values to SK07 for the best-fitting model parameters. The parameterparameter constraints as shown in Fig. 3 highlight the uncertainty in applying template sets, particularly for Arp 220 where SK07 suggest two likely best-fitting templates. For a more realistic test in determining the physical properties of a galaxy, we turn to the bright, lensed submillimetre galaxy SMM J2135-0102 (Swinbank et al. 2010) and apply the template SED library of Efstathiou et al. (2000, hereafter ERS00). The ERS00 starburst library consists of 44 templates with emission parametrized simply by the age of the starburst and the optical depth of molecular clouds where the new stars are forming. A normalization factor is required with the ERS00 templates to scale the emission from a single giant molecular cloud (on which the templates were formulated) to the entire system. This normalization factor roughly translates into an SFR with the model age according to SFR Norm e-
t/20 Myr

Downloaded from http://mnras.oxfordjournals.org/ at European Southern Observatory on July 17, 2015

as ERS00 assumes an exponentially decaying SFH with an e-folding time of 20 Myr. Applying the templates to the lensingcorrected SED of SMM J2135-0102, we find best-fitting parameters of 1960+282 , 56.6+10..5 Myr and 199.9+0.1.4 for model normal-44 -250 -15 2 ization, age and optical depth, respectively; this model is shown in the left-hand panel of Fig. 4. Using the standard FIRSFR relation of Kennicutt (1998), the best-fitting parameters imply an SFR of 192+28 M yr-1 . Comparatively, Swinbank et al. fit SMM -25 J2135-0102 with a two-temperature modified blackbody and derived an SFR of 210 50 M yr-1 using Kennicutt (1998), fully consistent with our results.




S. P. Johnson et al.

Downloaded from http://mnras.oxfordjournals.org/ at European Southern Observatory on July 17, 2015

Figure 3. Best-fitting SEDs for M82 (left) and Arp 220 (right) using the SK07 template library. Top: observed photometry along with the best-fitting models from SK07 (dashed red line) and those obtained with SATMC (solid black). An unreddened blackbody of temperature T = 2500 K has been added to M82 templates for consistency with SK07. Bottom: smoothed parameterparameter likelihood surfaces. The location of the maximum likelihood is marked by the large `X'. Contours are placed at ln(L) intervals corresponding to the 68 and 90 per cent confidence contours. The contours show truncation due to the parameter limits of the SK07 template library. Table 1. Comparison of best-fitting SK07 parameters using standard techniques (SK07) and SATMC. For M82, we add an unreddened blackbody of temperature T = 2500 K to the SK07 library for consistency. Note also that SK07 suggest two models for Arp 220, one with R = 1and AV = 120 and another with R = 3and AV = 72 (the latter being the accepted best-fitting template), which are represented in our fitting through the large uncertainties in R and AV . For each model, we also report the log-likelihood (ln(L)) of the fit. Errors are estimated at the 68 per cent confidence level. Name Ltot M82 Arp 220 10

R 0.35 3

SK07 AV LOB /Ltot 36 72 0.4 0.4


n 10 000 10 000

Ltot 1010. 10
12. 5+0.2 -0.1 1+0.1 -0.2

R 0.38+ -
0. 0. 2.84+0. -1. 57 03 16 70

AV 35.2+ -
33. 14. 70.1+69. -32. 9 8 5 2

LOB /Ltot 0.418+ -
0. 0. 0.410+0. -0. 167 018 176 010

n 9710 190
+290 -3907 +6883 -90

ln(L) -1218.7 -132.4


As a final verification of SATMC, we create a simulated SED with the SED synthesis routine GRASIL and attempt to recover the input parameters. GRASIL allows for numerous different parameters to be specified that will then determine the physical scale, chemical evolution and various attributes including dust content and inclusion of a starburst component. For a full description of GRASIL and its parameters, we refer the reader to Silva et al. (1998). For our synthetic galaxy, we set the galaxy age and total mass to 2 Gyr and 1012 M , respectively. A moderate amount of dust was included in the form of a dust mass of 109 M . The optical depth of UV/optical photons depends on the ratio of the mass of the molecular clouds in which new stars are formed and their size; we parametrize this as

the `CLOUD_RATIO' with an input value of 1000 M pc-2 . We then allow the simulated SED to undergo a `merger' at 1.95 Gyr. In effect, we set two SFH, one where the simulated galaxy is folk lowing a standard Schmidt SFR (SFR = sch Mgas where k is fixed at 1 and sch = 0.3) and one for the `merger-triggered' starburst which will convert 1010 M of gas into stars over a 50 Myr period following an exponentially decaying SFH with e-folding time of 50 Myr. The two SFHs are then combined before producing the final SED. Table 2 provides a summary of our adopted parameters. The resulting GRASIL spectrum was then redshifted to z = 2 and sampled at simulated wavelengths of 1.25 m, 1.6 m, 2.2 m, 3.6 m, 4.5 m, 5.8 m, 8.0 m, 24 m, 250 m, 350 m, 500 m, 1.1 mm

SED Analysis Through Markov Chains


Figure 4. Left: best-fitting ERS00 model to the observed SED of SMM J2135-0102. Right: best-fitting Table 2. Comparison of input and recovered at the 68 per cent confidence level. Parameter Tgal CLOUD_RATIO M M M
final burst dust


model to our synthetic z = 2 SED.

Downloaded from http://mnras.oxfordjournals.org/ at European Southern Observatory on July 17, 2015

parameters with


Errors are estimated

Description Galaxy age (Gyr) Mass/Radius2 for molecular clouds M pc-2 Galaxy mass accumulated over Tgal (M ) Mass of stars formed during starburst (M ) Dust mass (M ) Efficiency of Schmidt SFR ln(L) of fit

Input 2.0 1000 1e12 1e10 1e9 0.3

Best-fitting values
0.01 0.01 +33.5 990.6-35.7 9.99+0.16 e11 -0.10 1.02+0.01 e10 -0.04 9.98+0.21 e8 -0.30 0.300+0.002 -0.006

2.00+ -




and 21 cm corresponding to the wavelength coverage of the nearIR JHK bands, Spitzer/IRAC+MIPS, Herschel, AzTEC 1.1 mm and Very Large Array (VLA) observations. Errors in each of these bands were chosen to be representative of similar observations, from 0.010.1 Jy in the near/mid-IR and radio bands to 0.11.0 mJy for the (sub)mm. A small amount of random noise generated from the uncertainty in each band was added to the simulated data points prior to fitting. The fitting process for the simulated SED is then the same as previous cases though we now use GRASIL to produce SEDs rather than refer to a template library. The results of our fitting are shown in the bottom-right panel of Fig. 4 and Table 2. In addition to the synthetic GRASIL template shown here, we have also tested SATMC using several different combinations of input parameters and error handling (e.g. errors proportional to the simulated observations or representative errors as used previously) which all show similar agreement between the input and fitted SEDs. The excellent agreement between standard least-squares fits, our synthetic input SEDs and the recovered parameters with SATMC shows that our MCMC-based SED fitting method is robust and reliable. 5 A PPLICATION O F

tion at (sub)mm wavelengths and relative faintness of their optical counterparts. In Johnson et al. (2013), we took samples of AzTEC 1.1 mm-detected SMGs in three of the most widely studied fields, GOODS-N, GOODS-S and COSMOS, and examined their X-ray and IR/radio counterparts for evidence of active galactic nuclei (AGNs). Our SED analysis, performed using SATMC, showed that the IR-to-radio SED was dominated by starburst activity whereas an AGN component, whose IR emission was defined by the Siebenmorgen et al. (2004) template library, contributed little to the overall emission but was necessary to match the X-ray luminosity prior. The starburst templates, provided by ERS00, generally underpredicted the submm flux which suggested the presence of an extended cold dust component. Presently, we show the results of modelling a subsample of the Johnson et al. (2013) X-ray-detected SMGs using GRASIL, selected from GOODS-S with additional optical and Herschel observations. 5.1 Sample selection For a full description of our X-ray-detected AzTEC SMG sample, we refer the reader to Johnson et al. (2013). In summary, the fields GOODS-N, GOODS-S and COSMOS contain the most comprehensive multiwavelength coverage from space and ground based facilities including XMMNewton, Chandra, Hubble Space Telescope (HST), Spitzer, Herschel and VLA. In Johnson et al., X-ray counterparts were determined for the AzTEC SMGs detected in each of the three fields in order to constrain the relative contribution from an obscured AGN. Our preliminary analysis using the template libraries of Siebenmorgen et al. (2004) and ERS00 showed that while the X-ray detections indicated the presence of an AGN, the AGN component has a negligible contribution to the mid-to-far-IR SED as the observed X-ray luminosities limit the


For an illustration of the astrophysical impact of SATMC, we next set off to model a sample of high redshift submillimetre galaxies (SMGs) using GRASIL. SMGs are extremely luminous in the IR, 1012 L (Blain et al. 2004; Chapman et al. 2005), with LIR a redshift distribution peaking at z > 2. Given their high luminosities and IR(sub)mm emission, SMGs are generally believed to be very dusty systems containing powerful starbursts (SFRs 1000 M yr-1 ; Hughes et al. 1998; Coppin et al. 2006; Wei et al. 2009; Scott et al. 2010); however, the exact nature of these sources is still uncertain given the low angular resolu-


S. P. Johnson et al.
Table 3. Counterpart IDs and redshifts to the X-ray-detected SMGs in our sample. Counterparts are referenced according to their IAU identifier or catalogue ID number if not available. Redshifts given are the adopted spectroscopic/photometric values from Xue et al. (2011) with photometric redshifts in parentheses. SMM ID AzGS7b AzGS10 AzGS1 AzGS13 AzGS7 AzGS11 AzGS17a AzGS17b AzGS16 AzGS18 AzGS25 AzGS9 Chandra ID J033205.34-274644.0 J033207.12-275128.6 J033211.39-275213.7 J033212.23-274620.9 J033213.88-275600.2 J033215.32-275037.6 J033222.17-274811.6 J033222.56-274815.0 J033238.01-274401.2 J033244.02-274635.9 J033246.83-275120.9 J033302.94-275146.9 HST ID J033205.35-27454431 J033212.22-274620.7 J033215.27-275039.4 J033222.14-244811.3 J033222.57-274814.8 J033238.04-274403.0 J033244.01-274635.0 J033246.84-275121.2 J033303.05-275145.8 Spitzer ID J033205.35-244644.07 J033207.09-275128.96 J033211.36-245213.01 J033212.23-274620.66 J033213.85-275559.93 J033215.30-275038.31 J033222.16-274811.35 J033222.54-274814.94 J033238.01-274400.61 J033244.01-274635.24 J033246.83-275120.90 J033303.00-275146.27 Herschel ID Herschel72 Herschel435 Herschel280 Herschel375 Herschel1196 Herschel353 Herschel1316 Herschel1316 Herschel973 Herschel403 Herschel961 Herschel559 VLA ID Redshift (2.808) (1.829) (3.999) 1.034 (6.071) (0.682) (1.760) (2.535) 1.404 2.688 (1.101) (4.253)

660 336 390

178 188

Downloaded from http://mnras.oxfordjournals.org/ at European Southern Observatory on July 17, 2015

amount of IR emission from the AGN. We should point out, however, that the Siebenmorgen et al. (2004) AGN models assume dust reprocessing on scales much larger than the obscuring torus which may result in colder effective dust temperatures. Regardless, even traditional AGN torus models (e.g. Nenkova et al. 2008) showed negligible IR emission when coupled with the observed X-ray luminosities. We may therefore assume that our sample and the corresponding UV-through-IR SEDs are free from AGN contamination. While there may be radio excess in our sample due to possible AGN contribution, we currently lack the information to accurately model the AGN, star formation and possible radio jets. Though we have selected a particular subsample of the full SMG population, we are confident that the results obtained with this sample may be applied to all SMGs as emission from star formation will continue to dominate in non-X-ray-detected SMGs. For the present analysis, we limit our sample to the GOODS-S field where we have the deepest Chandra data to date (integrated exposure time 4 Ms) and excellent redshift coverage. The Chandra data are presented in Xue et al. (2011) though we use the custom source catalogue of Johnson et al. (2013) produced through standard Chandra Interactive Analysis of Observations (CIAO) reprocessing. Though our catalogue was created with a more stringent detection criteria (false detection probability of 10-6 ), it shows excellent agreement with that of Xue et al. Utilizing the X-ray data provides a statistically robust counterpart selection due to the lower X-ray source density than comparable optical or near-IR catalogues; for reference, 1 source of our current sample is expected to result from the false overlap of X-ray and AzTEC sources. The X-ray data also allow for unique counterpart identification in the optical/IR catalogues thanks to its very small (2 arcsec) positional accuracy. X-ray counterparts to the AzTEC GOODS-S sample (Scott et al. 2010) of SMGs are initially determined using a fixed search radius of 10 arcsec which is comparable to the average search radius used in Yun et al. (2012). This results in 16 X-ray-detected AzTEC SMGs where two have two potential X-ray counterparts. Recently, Hodge et al. (2013) have provided ALMA follow-up to 126 SMGs in the LABOCA ECDFS Submillimeter Survey (LESS; Wei et al. 2009) which allows an accurate counterpart identification and examination of potential multiple sources. However, only 4 of our 16 X-ray detected AzTEC sources are present in the LESS catalogue. The ALMA follow-up to the four common sources indicates that they are single systems (none of the potential multiples were detected by LABOCA) whose position agrees with our X-ray identifications to within 2 arcsec (see below). For those AzTEC sources with potential

multiple counterparts, we treat each source separately and make no attempt to split the AzTEC flux as we do not have enough information to determine how the AzTEC flux may be distributed between multiple sources. Once the X-ray counterparts have been determined, we cross-match our X-ray catalogue with publicly available VLA (Kellermann et al. 2008), Herschel (Oliver et al. 2012), Spitzer SIMPLE (Damen et al. 2011) and FIDEL (Spitzer PID30948; PI: Dickinson), and HST (Giavalisco et al. 2004) catalogues using a 2 arcsec search radius, consistent with the typical positional uncertainty of the X-ray sources; a 10 arcsec search radius was used for the Herschel catalogue given the much larger beam-size. Spectroscopic and photometric redshifts were obtained by crossmatching our X-ray catalogue to that of Xue et al. which compiles spectroscopic and photometric redshifts from 13 different groups (see Xue et al. for more details). For sources with photometric redshifts, the SED fits are fixed at their photometric values. One can include errors on the photometric redshift as a prior (Section 3.4); however, we only consider fixed redshifts here in order to maintain a consistent parametrization of the SED for all sources regardless of spec-z or photo-z quality and to avoid introducing redshift related correlations into the fitted parameters. To improve on the SED fitting of Johnson et al. (2013) and accurately model the far-IR peak, we limit our sample to only those with Herschel counterparts, putting our final sample at 12 X-ray-detected AzTEC SMGs. A summary of these objects is listed in Table 3. 5.2

parametrization and SED fitting

contains over 50 different parameters that control various aspects of the output UV-to-radio SED. These parameters cover the dust content, geometry of dust, stars and gas, the star formation rate history (SFRH) and more. For our SED fitting, we choose a set of eight free parameters, listed in Table 4, that we find to be representative of the major physical processes and have the most impact in our fits. These differ slightly from the standard input format for GRASIL but are easier to interpret. In addition to these free parameters, we have set additional constraints for the SFH and dust content. Here, we provide a brief motivation for our parametrization followed by the application to our SMG sample. The standard GRASIL implementation assumes that the SFH has two components: a Schmidt law of the form SFR(t ) = sch Mgas (t )k and a burst which can either be constant or have an exponential decay. The parameters sch , inf (rate of new gas in-falling on to the galaxy) and Mfinal (total gas mass Mgas and stellar mass M ) control

SED Analysis Through Markov Chains
Table 4. Description of parameters used for fitting our sample of SMGs with Parameter Tgal CLOUD_RATIO Mfinal Mburst inf etastart sch Mdust Description Galaxy age when `observed' (Gyr) Density of molecular clouds, sets effective 1 m optical depth (M Total galaxy mass (gas+stars) at Tgal (M ). Total mass of stars formed during starburst (M ) In-fall time-scale of IGM gas on to galaxy (Gyr) Escape time-scale for stars from their molecular clouds (Gyr) Efficiency of Schmidt SFR law Total dust mass pc-2 )


. Limiting range [0.5,GALAGE(z)] [1.0e-3,1.0e7] [1.0e10,1.0e13] [1.0e8,1.0e11] [0.1,10] [0.001,10] [1.0e-4,4] [1.0e8,1.0e10]

the Schmidt portion of the SFH (we fix the exponent k = 1 following Iglesias-Paramo et al. 2007). For the burst component, we assume that a mass Mburst of stars will be created from a starburst that occurs 50 Myr prior to the final galaxy age Tgal with an exponential decay of 50 Myr. Both the burst and Schmidt components follow the same initial mass function (IMF) which has been set to a Salpeter IMF (Salpeter 1955). GRASIL treats the SFH as a `closed-box' with no additional gas inflow outside of the continuous in-fall set by inf . We therefore modify the standard GRASIL SFH to create a `mergerlike' SFH following the prescription outlined in Section 4. In effect, the galaxy is passively evolving through the Schmidt law until the `burst' occurs which deposits additional gas to be turned into stars. This modified SFH aids in re-creating major/minor mergers or sudden peaks in gas inflow/accretion where the magnitude of the merger/in-fall event is measured with the ratio of Mfinal to Mburst ,the effective `merger ratio'. The second part of GRASIL is the handling of dust and radiative transfer to produce the synthetic SED. The primary influence to the IR portion of the SED is that of dust reprocessing which depends on

the dust mass Mdust and the stellar populations. These stellar populations are best parametrized according to their SFH (from above), the properties of their molecular clouds (CLOUD_RATIO) and the rate at which stars are able to escape their birth-sites (etastart). In addition to these free parameters, we fix the dust-to-gas ratio to that of the Milky Way (1/110; Draine & Lee 1984) and assume that dust/gas/stars follow a King radial density profile (see Silva et al. 1998). One possible limitation to GRASIL is that it has no prescription for any emission from an AGN; however, as we showed in Johnson et al. (2013), the AGN contribution is likely negligible to the mid-IR-to-radio SED of our SMG sample. 5.3 Do SMGs form a single population? With our GRASIL parameter space as defined above, it is straightforward to fit our SMGs sample with SATMC. A sample of the final SED produced by SATMC is shown in Fig. 5 with a summary of the fitting results in Fig. 6 and Tables 5 and 6. An important aspect to note from the parameterparameter covariances is that many parameters

Downloaded from http://mnras.oxfordjournals.org/ at European Southern Observatory on July 17, 2015

Figure 5. Best-fitting GRASIL model to the SED of J033211.39-275213.7. Shown here are the smoothed parameterparameter likelihood distributions with contours placed at the 68 and 90 per cent confidence intervals. The best-fitting model and observed SED are inlaid in the upper-right corner. Though parameters in GRASIL take on a continuous nature, the SFRH is computed on small, discrete time-steps which results in the `lumpy' confidence contours of Tgal .


S. P. Johnson et al.

Downloaded from http://mnras.oxfordjournals.org/ at European Southern Observatory on July 17, 2015

Figure 6. Best-fitting value of the fit.


models to the SMG sample. Included in each panel are the source ID, redshift used during the fitting (see Table 3), and the ln(L)

Table 5. Fitted parameters for our sample of 12 AzTEC SMGs using Source ID J033205.34-274644.0 J033207.12-275128.6 J033211.39-275213.7 J033212.23-274620.9 J033213.88-275600.2 J033215.32-275037.6 J033222.17-274811.6 J033222.56-274815.0 J033238.01-274401.2 J033244.02-274635.9 J033246.83-275120.9 J033302.94-275146.9 Tgal (Gyr) 2.03+ - 1 0 2 0 1 2 1 2 1 4 1 .55+ - .74+ - .62+ - .82+ - .16+ - .79+ - .41+ - .02+ - .03+ - .01+ - .00+ -
0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 2. 5. 0. 0. 0. 0. 0. 0. 0. 0. 2. 5. 5. 5. 09 35 01 01 04 01 01 01 01 01 03 00 01 01 01 04 01 01 01 01 05 00 28 00


. Errors on fitted parameters estimated at the 68 per cent confidence level. Mburst (1 0 9 M ) 1.91+ -
25.00 1.81 +10.2 21.2-8.4 0.11+1.04 -0.01 0.15+1.00 -0.05 9.66+3.36 -9.11 0.10+0.01 -0.01 29.5+2.9 -1.4 4.60+5.75 -1.90 0.11+0.59 -0.01 15.8+1.6 -5.3 0.84+0.02 -0.02 2.24+1.01 -1.11

etastart (Gyr) 1.07+ - 0.01+ - 0.12+ - 1.95+ - 0.94+ - 7.48+ - 0.21+ - 0.50+ - 3.83+ - 0.28+ - 2.78+ - 1.27+ -
0. 0. 0. 0. 0. 0. 0. 0. 9. 0. 2. 6. 0. 0. 0. 0. 6. 2. 0. 0. 0. 0. 8. 0. 14 35 03 01 02 01 10 22 03 58 52 62 05 01 44 18 17 16 05 02 27 12 70 65


Mfinal (1011 M )
19 06 05 15 29 18 01 01 43 03 01 01 01 08 18 05 06 65 01 07 01 01 11 03 8.3 9.5 +31.42 0.46-0.24 26.6+6.5 -4.3 5.34+0.33 -0.29 93.0+7.04 -48. 16.8+0.8 -0.9 0.201+0.02 -0.01 24.0+10.2 -12.6 3.52+0.18 -0.40 . 29.8+1672 -1. +0.07 2.95-0.07 26.3+1.3 -3.7

inf (Gyr) 0.77+ - 0.48+ - 1.94+ - 0.15+ - 8.12+ - 0.15+ - 1.83+ - 2.19+ - 0.28+ - 0.22+ - 0.75+ - 0.19+ -
1. 0. 0. 0. 0. 1. 0. 0. 1. 5. 0. 0. 0. 0. 3. 0. 0. 0. 0. 0. 0. 0. 0. 0. 81 45 11 27 46 05 01 01 87 83 01 01 01 05 10 51 03 03 01 10 01 01 15 05

0.06 0.16 e4 +11.70 2.08-1.24 e4 5.04+0.77 e2 -0.55 7.40+0.65 e3 -0.55 1.52+0.20 e2 -0.10 13.76+1.01 -1.33 4.68+0.18 e2 -0.13 1.71+0.35 e2 -0.25 1.69+0.13 e3 -0.16 7.18+0.49 e2 -0.53 9.95+0.20 e2 -0.39 1.96+0.11 e2 -0.15

Mdust (1 0 9 M )
+ - 8.79+ - 4.00+ - 1.71+ - 1.11+ - 3.97+ - 5.75+ - 0.87+ - 9.76+ - 6.02+ - 3.16+ - 1.32+ -

ln(L) -149.4 -0.4 -13.3 -169.6 -80.3 -1816.4 -155.8 -16.9 -31.6 -69.3 -772.4 -283.3

3.81+ - 0.16+ - 1.08+ - 2.79+ - 0.32+ - 0.02+ - 1.63+ - 0.15+ - 3.48+ - 0.18+ - 0.54+ - 0.46+ -

0. 2. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.

26.4+ -

1.37+ -


0. 0. 1. 5. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 2. 0. 0. 0. 0. 0. 0.

93 67 20 41 72 44 46 25 32 27 41 30 72 05 25 27 24 17 81 94 00 09 28 35

show tight constraints and large degeneracies generally do not exist in our adopted parametrization. The strongest correlations exist between Mfinal and sch which together set the stellar mass M . In the observed SEDs, the stellar mass is measured through the stellar bump which is most prominent in the IRAC bands. As the IRAC observations have very small flux errors ( 1 per cent), this limits parameter space to the narrow range of Mfinal and sch capable of reproducing the observed fluxes. When examining a population of sources, it is often helpful to view the SEDs as a composite and construct an empirical median template. This allows one to easily see trends common in the population including the approximate flux ranges covered by the models

and dominate aspects of the SEDs. Fig. 7 shows our composite restframe GRASIL models along with the median SED of Michalowski et al. (2010), who fit the SEDs of SMGs using a grid developed from GRASIL, and the composite star-forming SEDs of Kirkpatrick et al. (2012), derived from z 1 and 2 sources in GOODS-N and ECDFS with Spitzer IRS spectroscopy. Though Michalowski et al. and Kirkpatrick et al. normalize their fits at specific wavelengths to construct the composite, we have applied no re-normalization to match our SED set. The agreement between our SEDs with those of Mickalowski et al. and Kirkpatrick et al. is immediately apparent. From Kirkpatrick et al., our composite SED best agrees with the z 2 star-forming galaxy; not surprising as SMGs peak around

SED Analysis Through Markov Chains
Table 6. Derived attributes from the fitted parameters of Table 5. SFR and stellar mass M are provided as outputs from the GRASIL SFH calculation. The SFRs reported are the SFRs averaged over the age of the burst (50 Myr prior to Tgal ) with stellar masses given at Tgal . For comparison, we also include the SFR as derived from the FIR luminosity (SFRFIR ) of the best-fitting model following Kennicutt (1998). Specific star formation rates (SSFRs) are simply SFR/M . Errors on SFR and M are given at the 68 per cent confidence level.
Source ID (M J033205.34-274644.0 J033207.12-275128.6 J033211.39-275213.7 J033212.23-274620.9 J033213.88-275600.2 J033215.32-275037.6 J033222.17-274811.6 J033222.56-274815.0 J033238.01-274401.2 J033244.02-274635.9 J033246.83-275120.9 J033302.94-275146.9 SFR yr-1 ) SFRFIR yr-1 ) 1230.0 182.5 1892.8 84.4 1966.0 18.1 187.1 308.2 64.0 458.4 50.0 1018.3 (1011 M 23.14+ - 0.23+ - 6.80+ - 4.86+ - 8.50+ - 0.21+ - 0.36+ - 1.86+ - 3.21+ - 3.03+ - 1.83+ - 6.09+ - M )
8 9 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 .39 .65 .38 .15 .52 .33 .29 .27 .97 .45 .01 .01 .02 .04 .11 .28 .17 .38 .40 .13 .05 .04 .11 .54


sch and Mtotal had shown tight correlations in individual fits, the composite fits show no signs of the original correlations. 5.4 SMGs and the `main sequence' of galaxy formation It is often believed that SMGs must result from powerful starbursts possibility triggered by major mergers. However, these results show a wide range in derived SFRs (202000 M yr-1 ) with many below the 1000 M yr-1 typically assumed for SMGs. Furthermore, the `merger ratio' (Mburst /Mfinal ) would also suggest that major mergers are uncommon. This is in agreement with Hayward et al. (2011) who find that starbursts are highly inefficient at boosting the submm flux (e.g. a 16 boost in SFR produces 2 increased submm flux). Indeed, the common trend in our sample is that the bulk of the stellar mass is built from the Schmidt component rather than a `merger' triggering a starburst. A more detailed study including a larger sample is needed to further clarify this point. The sources with low SFR could result from source blending giving rise to an enhanced submm flux as also suggested by Hayward et al. Karim et al. (2013) have shown that the brightest LABOCA sources (S850 m 12 mJy) are composed of multiple, fainter SMGs. For our AzTEC SMG sample, however, the corresponding bright flux limit is 6.7 mJy (assuming R = S850 m /S1.1mm = 1.8; Chapman et al. 2009) which only applies to the source J033211.39-275213.7, a source with no multiple X-ray detection and derived SFR 2000 M yr-1 . Furthermore, only 23 out of 12 of our SMGs have multiple potential counterparts (Table 3), none of which correspond to a low SFR source. While source blending is certainly problematic for correctly interpreting multiwavelength SED fits, it does not have a significant influence in our sample. Expanding on this point, we show the specific star formation rate (SSFR = SFR/M ) with respect to redshift for our SMG sample in Fig. 9. Elbaz et al. (2011) have suggested an evolution in SSFR for quiescent and starbursting galaxies which follow the general forms SSFR = 26 t-2.2 Gyr-1 and SSFR > 56 t-2.2 Gyr-1 , where t is the cosmic time, for quiescent or `main-sequence' galaxies and starbursts, respectively. We should point out, however, that the stellar masses used to derive the Elbaz et al. relation are of the order of 5 times lower than stellar masses in our sample. As a result, only three of our SMGs lie above the `main-sequence' relation of Elbaz et al. Despite this discrepancy, our SFRs and stellar masses agree with Michalowski et al. (2010) and the predictions of Hayward et al. (2011); differences in stellar mass derivations can be


SSFR (Gyr-1 )
0 0 +42 19.26-13 3.07+0 -0 0.11+0 -0 3.33+0 -0 1.46+0 -0 14.60+6 -0 2.04+1 -0 0.19+0 -0 2.46+0 -0 0.41+0 -0 1.57+0 -0

419.3 227.6 +235.5 428.8-140.7 2086.5+2026.2 -90. 53.2+196.6 -2. 2829.0+166..7 -304 7 30.3+0..6 -0 5 518.0+217.8 -0.4 379.2+1399.2 -13. 59.6+105.5 -4. 745.7+76..7 -75 9 74.9+0..2 -5 2 955.5+24..6 -18 1

950.8+ -

0.41+ -

.42 .18 .82 .20 .19 .13 .04 .01 .10 .42 .01 .01 .59 .40 .05 .08 .04 .01 .24 .36 .01 .03 .19 .06

Downloaded from http://mnras.oxfordjournals.org/ at European Southern Observatory on July 17, 2015

z 2.5. Our composite also shows excellent agreement with the median Michalowski et al. SED. There are some minor discrepancies between the composites though they are likely attributed to differences in parametrizing the SED. Compared to the Kirkpatrick et al. composites, our SFRs and stellar mass estimates are higher by a factor of 2 which can readily be explained by the SFH adopted. Regardless, our composite SED is not only able to recover the medians of Kirkpatrick et al. and Michalowski et al. but also encompass the scatter in the observations and dispersion in similar template sets. Despite the similar shapes and well-defined empirical median of our fitted SEDs, the SMGs in our sample are not heterogeneous in nature and, in fact, occupy large regions of parameter space. Fig. 8 shows that while some areas of parameter space are well constrained (e.g. inf versus Tgal ), there is no single set of parameters that may describe our sample as an individual population. Of all our adopted parameters, Mdust shows the tightest constraints (Table 4). This is primarily due to the inclusion of the (sub)mm observations which are highly sensitive to the total dust mass. Though the parameters

Figure 7. Compiled SEDs for the SMG sample. The median SED is given by the solid black line with the grey shaded region indicating the dynamic range between individual SEDs. Left: The composite SEDs overlaid with the average star-formation templates of Kirkpatrick et al. (2012). The photometry used in constructing the z 2 SF sample has been included to demonstrate the scatter in the data. Right: The composite SEDs overlaid with median SMG template of Michalowski et al. (2010) with the dashed lines encompassing the dispersion in the Mickalowski et al. models. Note that these models have not been re-normalized to fit our model set.


S. P. Johnson et al.

Downloaded from http://mnras.oxfordjournals.org/ at European Southern Observatory on July 17, 2015

Figure 8. Compiled histogram distribution of accepted steps from all fits. Note that the histogram distributions differ slightly from the likelihood distributions but maintain the same general characteristics due to the multidimensional nature of the data volume.

of a mixture of starburst and quiescent galaxies of various ages and SFRHs. 5.5 Comparing model parametrizations When developing SED models and deriving the associated parameters, there are certain attributes that must remain fixed in order to keep a consistent parametrization. Examples of these attributes include the IMF, spatial properties like the disc/bulge scaleheights and radial distribution, and the form of the SFH. Altering any of these properties changes the very nature of the models and is in no way continuous for MCMC samplers. In Bayesian statistics, it is possible to compare models of different parametrizations through the Bayes factor to determine which model offers the `best' fit (e.g. Jefferys 1961; Weinberg 2013; Weinberg, Yoon & Katz 2013). The Bayes factor follows from Bayes' Theorem and is derived from the ratio of the posteriors for the two models M1 and M2 : P ( M 1 ) P ( D |M 1 ) P (M 1 | D ) = . P (M2 | D ) P ( M 2 ) P ( D |M 2 ) (9)

Figure 9. SSFR as a function of redshift. The area within the dashed lines represents the main sequence of galaxies from Elbaz et al. (2011) with galaxies above this limit defined as starbursts. The lower long dashed line represents the approximate main-sequence relation for our sample accounting for the factor of 5 difference in stellar masses between the Elbaz et al. sample and those presented here.

attributed to differences in the parametrization of the SFRH, IMF (see Section 5.4) and stellar evolution tracks (see also Michalowski et al. 2012). Given the uncertainties in determining the stellar mass, we caution against overinterpretation of the SSFRs. However, adjusting for the relative differences in stellar masses, we see that the majority of our sample lie above or near the main-sequence relation while spanning a large range of SSFRs. Overall, these results suggest that there is not a single population that defines our SMG sample; instead, they are homogeneous in origin composed

The Bayes factor is defined as the second term on the right-hand side; the first term is the ratio of the prior probability of each model and is typically set to unity as the models are considered under equal weight. Explicit calculation of the Bayes factor (B12 )involves integrating the likelihood distributions of each model (L1 ( D|1 )and L2 ( D|2 )) with the prior parameter distributions (1 (1 )and 2 (2 )) according to B12 = P ( D |M 1 ) P ( D |M 2 ) 1 (1 |M1 )L1 ( D|1 ,M1 )d1 . 2 (2 |M2 )L2 ( D|2 ,M2 )d2 ) (10)

While there are MCMC and Bayesian based tools available that compute the Bayes factor (e.g. the Bayesian Inference Engine;

SED Analysis Through Markov Chains
Table 7. Comparison of fitted parameters for the source J033213.88-275600.2 using a Salpeter and Chabrier IMF. IMF Tgal (Gyr) 0.82
+ - 0.90+ - 0. 0. 0. 0. 01 01 01 01


etastart (Gyr) 0.94+ -
9.03 0.58 +9.01 0.94-0.38


Mfinal (1011 M ) 93.0+ -
7.0 48.4 +11.7 88.3-7.2

Mburst (1 0 9 M )
3.36 9.11 +23.69 1.16-1.06

inf (Gyr) 8.12+ -
1.87 5.83 +0.15 0.19-0.09

CLOUD_RATIO (M pc-2 ) 1.52+ -
0.20 0.10 +0.17 1.30-0.03

Mdust (1 0 9 M ) 1.11 1.44
+ - + - 0. 0. 0. 0. 32 27 20 52

ln(L) -80.3 -102.9

Salpeter Chabrier

0.32+ -

0.43 0.03 +0.06 0.52-0.11

9.66+ -

e2 e2

Table 8. Comparison of derived parameters for J033213.88-275600.2 using Salpeter and Chabrier IMFs. IMF (M Salpeter Chabrier SFR yr-1 )
7 7 1 5

M (1011 M )
0. 0. 13.47+1. -1.

SSFR (Gyr-1 ) 3.33 2.89

166. 304. 3893.4+192. -122.

2829.0+ -

8.50+ -

97 45 79 94

Figure 10. Best-fitting GRASIL SEDs for J033213.88-275600.2 using Salpeter (solid) and Chabrier (dashed) IMFs.

Weinberg 2013), such methods are not currently available in SATMC. Instead, we may follow a similar prescription by repeating the SED fits after changing the desired parametrization(s). We provide a brief example of this process by comparing the results obtained with the Salpeter IMF used in our initial fits and those obtained with a Chabrier IMF (Chabrier 2003). Tables 7 and 8 and show that for the source J033213.88-275600.2 the fitted parameters remain nearly constant except for inf which causes an increase in the derived SFRs and stellar mass while the fitted SEDs are nearly identical (Fig. 10). We caution, however, that without results from a full sample and direct computation of the Bayes factor, we cannot differentiate between the two IMFs at this time. Without the Bayes factors, one may include secondary observations, i.e. spectroscopyand SED-independent estimates of fitted and/or derived parameters, to derive priors and determine which, if any, of the fits best describe all observations. 5.6 Application of additional observations to SED fits As mentioned in Section 2.2, whenever an SED synthesis routine is used in SATMC, the SEDs for each step are saved along with its location in parameter space and likelihood. Using these `SED steps', one can examine the influence of additional or improved photometric data simply by referring back to equations (5) and (8)

without having to re-run the entire MCMC fit. This can be exceptionally useful given the long typical run time of a single fit and is best used for examining the influence on parameter space when including low signal-to-noise observations (e.g. LABOCA) or observations one believes will truncate regions of parameter space (e.g. priors on parameters). To examine the full impact of improved or additional observations, it is necessary to repeat the fits as the sampling of likelihood space is wholly dependent on the initial set of observations and relative uncertainties. As an example, we have created simulated GRASIL SEDs based on the best-fitting parameters of J033211.39-275213.7 (z 4). The simulated SED is produced using the same formalism and parametrization described in Sections 4 and 5.2 and has been sampled at wavelengths corresponding to the actual observations (3.6 m, 4.5 m, 5.8 m, 8.0 m, 24 m, 250 m, 350 m, 500 m, 1.1 mm). Each band is then assigned the measured flux uncertainties. For one simulation, we have increased the AzTEC signal-to-noise by a factor of 10 ( 1.1 mm 0.1 mJy) similar to the level of improvement expected with Large Millimetre Telescope (LMT), ALMA and future submm telescopes whereas the control simulation maintains the original AzTEC signal-to-noise ( 1.1 mm 1 mJy). Fig. 11 shows the 68 per cent contours for the GRASIL fits to the simulated SEDs where we immediately see tighter constraints on fitted parameters for the improved AzTEC photometry. Referring back to Fig. 6, we see that the AzTEC photometry is one of the few points beyond the thermal dust peak and along the RayleighJeans tail of the dust blackbody distribution. At lower redshifts, the Herschel/SPIRE bands become more prominent in the beyond the dusk peak which serves to diminish the impact of improved (sub)mm photometry. For the most distant galaxies, however, high signal-to-noise (sub)mm photometry is instrumental for fully constraining the SED shape and its driving parameters, namely its SFRH (Tgal , Mtotal , Mburst , sch , inf ) and dust content Mdust . The future observations of SMGs with LMT and ALMA will therefore be paramount to our understanding of high-redshift galaxies. 6 S UMMAR Y We have presented the general purpose MCMC-based SED fitting tool SATMC. Utilizing MCMC algorithms, the code is able to take any set of SED templates or models of the user's choice and return the requested best-fitting parameter estimates in addition to the posterior parameter distribution which can be used to construct confidence levels and unveil parameter correlations. In testament to SATMC's flexibility, we have provide a series of test cases comparing SATMC with traditional least-squares methods where SATMC recovers best-fitting values similar to those from traditional methods. Furthermore, we show that SATMC is capable of reproducing the set of input parameters from a simulated SED which serves to prove SATMC's adaptability and reliability. Highlighting SATMC's scientific value, we also provide the bestfitting GRASIL models to a sample of AzTEC SMGs in the GOODS-S field. These fits indicate that while some parameters exhibit obvious and strong correlations in individual SED fits, such correlations

Downloaded from http://mnras.oxfordjournals.org/ at European Southern Observatory on July 17, 2015


S. P. Johnson et al.

Downloaded from http://mnras.oxfordjournals.org/ at European Southern Observatory on July 17, 2015

Figure 11. 68 for the original are marked by signal-to-noise

per cent error contours for a simulated GRASIL SED based on the fitted spectrum of J033211.39-275213.7 at z 4. The contours are shown AzTEC photometric errors (blue dashed) and improved by a factor of 10 (red solid). The location of the maximum likelihoods for each case the X's with values used in construction of the simulated SEDs marked by the dotted lines. For the high-z simulation, improving the AzTEC provides significantly greater constraints on parameter space.

are absent when constructing the population. Furthermore, these sources appear to span a wide dynamic range of fitted and derived parameters, specifically when considering their derived SFRs (30 3000 M yr-1 ) and stellar masses (1010 1012 M ); these results agree both with the previous GRASIL-based SEDs of Michalowski et al. (2010) and simulations of SMGs (Hayward et al. 2011). We also provide a brief demonstration of the impact of improved photometric data on the best-fitting results. Increasing the signal-tonoise of the AzTEC photometry provides significantly greater constraints on fitted parameters and highlights the need to obtain deeper (sub)mm data. As observations, the complexity of models and the desire to obtain more detailed information from model fits increase, the approach of tools like SATMC will be required in both observational and theoretical astronomy in order to refine our understanding of the Universe and its constituents.

Acquaviva V., Gawiser E., Guaita L., 2011, ApJ, 737, 47 Acquaviva V., Gawiser E., Guaita L., 2012, in Tuffs R. J., Popescu C. C., eds, Proc. IAU Symp. 284, The Spectral Energy Distribution of Galaxies. Cambridge Univ. Press, Cambridge, p. 42 Atchade Y. F., Rosenthal J. S., 2005, Bernoulli, 11, 815 Benitez N., 2000, ApJ, 536, 571 Blain A. W., Chapman S. C., Smail I., Ivison R., 2004, ApJ, 611, 52 Bolzonella M., Miralles J.-M., Pello R., 2000, A&A, 363, 476 Brammer G. B., van Dokkum P. G., Coppi P., 2008, ApJ, 686, 1503 Bruzual A., Charlot S., 1993, ApJ, 405, 538 Bruzual A., Charlot S., 2003, MNRAS, 344, 1000 Calzetti D., Armus L., Bohlin R. C., Kinney A. L., Koornneef J., StorchiBergmann T., 2000, ApJ, 533, 682 Chabrier G., 2003, PASP, 115, 763 Chapman S. C., Blain A. W., Smail I., Ivison R. J., 2005, ApJ, 622, 772 Chapman S. C., Blain A. W., Ibata R., Ivison R. J., Smail I., Morrison G., 2009, ApJ, 691, 560 Chib O., Jeliazkov I., 2001, J. Am. Stat. Assoc., 96, 270 Coppin K. E. et al., 2006, MNRAS, 372, 1621 Dahlen T. et al., 2010, ApJ, 724, 425 Dahlen T. et al., 2013, ApJ, ApJ, 775, 93 Damen M. et al., 2011, ApJ, 727, 1 Draine B. T., Lee H. M., 1984, ApJ, 258, 89 Earl D. J., Deem M. W., 2005, Phys. Chem. Chem. Phys., 7, 3910 Efstathiou A., Rowan-Robinson M., Siebenmorgen R., 2000, MNRAS, 313, 734 (ERS00) Elbaz D. et al., 2011, A&A, 533, 119 Feigelson E. D., Nelson P. I., 1985, ApJ, 293, 192 Gawiser E., 2000, New Astron. Rev., 53, 50 Gelman A., Rubin D., 1992, Stat. Sci., 7, 457

A C KNO W LEDGEMENTS We would like to thank V. Acquaviva, the CANDELS collaboration and A. Kirkpatrick for their assistance in developing and providing data sets for testing and improving our MCMC algorithm. We would also like to thank the anonymous reviewer for their helpful comments on improving the manuscript. This work has been funded in part by the North East Alliance under the National Science Foundation (NSF) grant HRD 0450339, NSF grants AST-0907952 and AST-0838222 and CXO grant SAO SP1-12003X. KSS is supported by the National Radio Astronomy Observatory, which is a facility of the NSF operated under cooperative agreement by Associated Universities, Inc.

SED Analysis Through Markov Chains
Gelman A., Roberts G. O., Gilks W. R., 1996, in Bernado J. M., Berger J., Dawid A. P., Smith A. F. M., eds, Bayesian Statistics 5. Oxford Univ. Press, Oxford, p. 599 Geman S., Geman D., 1984, IEEE Trans. Pattern Anal. Mach. Intell., 6, 721 Geweke J., 1992, in Bernardo J. M., Berger J., Dawid A. P., Smith A. F. M., eds, Bayesian Statistics 4. Oxford Univ. Press, Oxford, p. 169 Geyer C. J., 1991, in Keramidas E. M., ed., Proceedings of the 23rd Symposium on the Interface, Computing Science and Statistics. Am. Stat. Assoc., New York, p. 156 Geyer C. J., 1992, Stat. Sci., 7, 473 Giavalisco M. et al., 2004, ApJ, 600, 93 Gilks W. R., Richardson S., Spiegelhalter D., 1996, in Gilks W. R., Richardson S., Spiegelhalter D., eds, Markov Chain Monte Carlo in Practice. Chapman and Hall, London, p. 64 Hastings W. K., 1970, Biometrika, 57, 97 Hayward C. C., Keres D., Jonsson P., Narayanan D., Cox T. J., Hernquist L., 2011, ApJ, 743, 159 Hodge J. A. et al., 2013, ApJ, 768, 91 Hughes D. H. et al., 1998, Nat, 394, 241 Iglesias-Paramo J. et al., 2007, ApJ, 670, 279 Isobe T., Feigelson E. D, Nelson P. I., 1986, ApJ, 306, 490 Jefferys H., 1961, Theory of Probability, 3rd edn. Oxford Univ. Press, Oxford Johnson S. P. et al., 2013, MNRAS, 431, 662 Karim A. et al., 2013, MNRAS, 432, 2 Kellermann K. I., Fomalont E. B., Mainieri V., Padovani P., Rosati P., Shaver P., Tozzi P., Miller N., 2008, ApJS, 179, 71 Kennicutt R. C., 1998, ARA&A, 36, 189 Kirkpatrick A. et al., 2012, ApJ, 759, 139 Kodama T., Bell E. F., Bower R. G., 1999, MNRAS, 302, 152 Mackay D., 2003, Information Theory, Inference, and Learning Algorithms. Cambridge Univ. Press, Cambridge Metropolis N., Rosenbluth A. W., Rosenbluth M. N., Teller A. H., Teller E., 1953, J. Chem. Phys., 21, 1087 Michaloski M. J., Watson D., Hjorth J., 2010, ApJ, 712, 942


Michalowski M. J., Dunlop J. S., Cirasuolo M., Hjorth J., Hayward C. C., Watson D., 2012, A&A, 541, 85 Nenkova M., Sirocky M. M., Nikutta R., Ivezic Z., Elitzur M., 2008, ApJ, 685, 160 ~ Noll S., Burgarella D., Giovannoli E., Buat V., Marcillac D., Munoz-Mateos J. C., 2009, A&A, 507, 1793 Oliver S. J. et al., 2012, MNRAS, 424, 1614 Pirzkal N., Rothberg B., Nilsson K. K., Finkelstein S., Koekemoer A., Malhotra S., Rhoads J., 2012, ApJ, 748, 122 Raftery A. E., Lewis S., 1992, in Bernardo J. M., Berger J., Dawid A. P., Smith A. F. M., eds, Bayesian Statistics 4. Oxford Univ. Press, Oxford, p. 763 Roberts G. O., Rosenthal J. S., 2001, Stat. Sci., 16, 351 Roberts G. O., Gelman A., Gilks W., 1997, Ann. Appl. Probab., 7, 110 Sajina A., Scott D., Dennefeld M., Dole H., Lacy M., Lagache G., 2006, MNRAS, 369, 939 Salpeter E. E., 1955, ApJ, 121, 161 Sawicki M., 2012, PASP, 124, 1208 Scott K. S. et al., 2010, MNRAS, 405, 2260 Serra P., Amblard A., Temi P., Nurgarella D., Giovannoli E., Buat V., Noll S., Im S., 2011, ApJ, 740, 22 Siebenmorgen R., Krugel E., 2007, A&A, 461, 445 (SK07) Siebenmorgen R., Freudling W., Krugel E., Haas M., 2004, A&A, 421, 129 Silva L., Granato G. L., Bressan A., Danese L., 1998, ApJ, 509, 103 Swinbank A. M. et al., 2010, Nat, 464, 733 Verde L. et al., 2003, ApJS, 148, 195 Wei A. et al., 2009, ApJ, 707, 1201 Weinberg M. D., 2013, MNRAS, 434, 1736 Weinberg M. D., Yoon I., Katz N., 2013, preprint (arXiv:1301.3156) Xue Y. Q. et al., 2011, ApJS, 195, 10 Yun M. S. et al., 2012, MNRAS, 420, 957

Downloaded from http://mnras.oxfordjournals.org/ at European Southern Observatory on July 17, 2015

A This paper has been typeset from a TEX/L TEX file prepared by the author.