Exploring the physical properties of local star-forming ULIRGs from the ultraviolet to the infrared
Key Words.:Galaxies: evolution – galaxies: fundamental parameters – galaxies: starburst – galaxies: ISM.
We present an application of the da Cunha, Charlot & Elbaz (2008) model of the spectral energy distribution (SEDs) of galaxies from the ultraviolet to far-infrared, to a small pilot sample of purely star-forming Ultra-Luminous Infrared Galaxies (ULIRGs). We interpret the observed SEDs of 16 ULIRGs using this physically-motivated model which accounts for the emission of stellar populations from the ultraviolet to the near-infrared and for the attenuation by dust in two components: an optically-thick starburst component and the diffuse ISM. The infrared emission is computed by assuming that all the energy absorbed by dust in these components is re-radiated at mid- and far-infrared wavelengths. This model allows us to derive statistically physical properties including star formation rates, stellar masses, as well as temperatures and masses of different dust components and plausible star formation histories. We find that, although the ultraviolet to near-infrared emission represents only a small fraction of the total power radiated by ULIRGs, observations in this wavelength range are important to understand the properties of the stellar populations and dust attenuation in the diffuse ISM of these galaxies. Furthermore, our analysis indicates that the use of mid-infrared spectroscopy from the Infrared Spectrograph on the Spitzer Space Telescope is crucial to obtain realistic estimates of the extinction to the central energy source, mainly via the depth of the 9.7-m silicate feature, and thus accurately constrain the total energy balance. Our findings are consistent with the notion that, in the local Universe, the physical properties of ULIRGs are fundamentally different from those of galaxies with lower infrared luminosities and that local ULIRGs are the result of merger-induced starbursts. While these are well-established ideas, we demonstrate the usefulness of our SED modelling in deriving relevant physical parameters which provide clues to the star formation mode of galaxies.
Ultra-Luminous Infrared Galaxies (ULIRGs) are galaxies with infrared luminosities higher than . They emit most of their energy (over 90 per cent) in the infrared, which implies that they are heavily dust-obscured. To power their huge infrared luminosities, these galaxies must be undergoing intense star formation, forming new stars at rates of the order of yr. In the local Universe, such large starbursts are thought to result from major interactions or mergers of gas-rich galaxies (e.g. Murphy et al. 2001; Sanders & Mirabel 1996 and references therein). This is supported by the fact that local ULIRGs have typically disturbed morphologies (e.g. Clements et al. 1996; Duc et al. 1997; Rigopoulou et al. 1999). Another possible source of energy in ULIRGs is the accretion of large quantities of gas onto a central supermassive black hole in an active galactic nucleus (AGN; e.g. Sanders et al. 1988). In such cases, dust in an optically-thick torus surrounding the AGN absorbs large amounts of energy and re-radiates it at longer wavelengths, contributing to the large infrared luminosities of ULIRGs. In most cases, there is likely to be a combination of both star formation and AGN contributing to the emission from ULIRGs, and the relative contribution by the two processes has been extensively investigated (e.g. Lutz et al. 1998; Genzel et al. 1998; Laurent et al. 2000; Armus et al. 2007; Spoon et al. 2007; Veilleux et al. 2009).
Understanding the formation and evolution of ULIRGs is key to understanding the cosmic evolution of galaxies. ULIRGs play a significant role in galaxy evolution at high redshifts and they are major contributors to the infrared luminosity and star formation density at redshifts (e.g. Chary & Elbaz 2001; Le Floc’h et al. 2005). Local ULIRGs may also provide analogs to dusty galaxy populations at high redshifts (e.g. Blain et al. 2002; Daddi et al. 2007; Dey et al. 2008). To understand the evolution of ULIRGs in a cosmological context, the study of the physical properties of these system both at low and high redshifts is required.
With the advent of a number of sensitive space observatories and new instrumentation on ground-based telescopes, a wealth of new data on local samples of ULIRGs have been compiled over recent years. The analysis of such data has focused mainly on the infrared range, where ULIRGs emit most of their energy. Several models have been proposed to interpret the infrared spectral energy distributions of ULIRGs (e.g. Klaas et al. 2001; Farrah et al. 2003; Vega et al. 2008; Marshall et al. 2007). However, such models are not easily applicable in the statistical study of large samples of galaxies and they often neglect the optical emission of ULIRGs, which may carry important information about their intermediate-age and old stellar populations and hence their star formation histories.
Recently, da Cunha et al. (2008) have presented a simple, physically motivated model that allows the interpretation of the integrated ultraviolet, optical and infrared spectral energy distributions (SEDs) of galaxies in terms of fundamental physical parameters. This model was used in da Cunha et al. (2008) to derive median-likelihood estimates of parameters such as the star formation rate, stellar mass, dust attenuation and dust mass from the observed SEDs of 66 galaxies in the Spitzer Infrared Nearby Galaxy Survey (SINGS, Kennicutt et al. 2003). In another study, da Cunha et al. (2010) have used this model to investigate the relation between star formation and dust content of a sample of 3258 low-redshift IRAS galaxies detected by Sloan Digital Sky Survey (SDSS) with complementary observations by GALEX (ultraviolet) and 2MASS (near-infrared). These studies have provided valuable insight into the physical properties of local galaxies with moderate total infrared luminosities ().
The method of da Cunha et al. (2008) was designed to interpret SEDs with sufficiently wide wavelength coverage (from the UV to the far-IR) to allow for accurate extinction corrections and energy balance. The main assumptions are that the energy input originates only from stars and all energy absorbed by dust in the ultraviolet and optical is re-radiated in the infrared. In this paper, we investigate how the model of da Cunha et al. (2008) can be calibrated to interpret multi-wavelength observations of star-forming ULIRGs, which have much higher infrared luminosities and for which only a small fraction of the energy is detected in the ultraviolet and optical. This is the first step of a larger study in which we aim at expanding our methodology to interpret the physical properties of samples of galaxies covering a wide range of infrared luminosities, star formation and AGN activities, such as the Great Observatories All-Sky LIRG Survey (GOALS; Armus et al. 2009) and the local ULIRG sample of Desai et al. (2007).
In the present study, we start by analysing the SEDs of 16 purely star-forming ULIRGs from the primary sample of Desai et al. (2007). In a forthcoming paper, we will extend our analysis to galaxies with AGN contribution.
This paper is organized as follows. In Section 2, we describe our sample of local ULIRGs and the multi-wavelength data used in our study. We detail our method to derive statistical constraints on the physical parameters of our galaxies from fits to their observed spectral energy distributions in Section 3. In Section 4, we provide the fits to the SEDs of the galaxies and the estimates of their physical parameters. In Section 5 we show that our results are consistent with previous findings regarding the formation and evolution of ULIRGs, and compare the physical properties of local galaxies in two bins of total infrared luminosity: and . Our summary and conclusions are presented in Section 6. When necessary, we adopt the following cosmology: km s Mpc, and .
2 The sample
We select our galaxies from a large sample of 107 local ULIRGs with observed in the region with the Spitzer/IRS (Armus et al. 2007; Desai et al. 2007). The sources were primarily selected from the IRAS 1 Jy survey (Kim & Sanders 1998), the IRAS 2 Jy survey (Strauss et al. 1992) and the FIRST/IRAS radio-far-infrared sample of Stanford et al. (2000). We note that preliminary theoretical interpretation of infrared observations of this sample was preformed by Armus et al. (2007) using the spectral decomposition method of Marshall et al. (2007).
2.1 Selection of our star-forming ULIRGs
Given the theoretical tools we have in hand, in the present study, our goal is to explore how the physical properties of purely star-forming ULIRGs can be characterized using the model of da Cunha et al. (2008). Modelling the SEDs of ULIRGs presents a challenge since theses systems have infrared luminosities a factor of 10 greater than what had been studied so far (da Cunha et al. 2008, 2010) and very high infrared-to-optical ratios. Typically more than 90% of the energy output of ULIRGs is emitted in the infrared. This energy originates not only from young massive, heavily obscured starbursts, but in some cases from an AGN. Identifying the presence of an AGN in the dust-enshrouded nucleus of a ULIRG is often challenging, in particular in the optical, where the extinction effects are more prevalent. The best diagnostic is the detection of hard X-ray emission, or high ionization lines which cannot be excited by stellar photospheric emission. The task is somewhat easier in the infrared, and various methods have been developed to reveal an AGN and quantify its contribution to the bolometric emission of galaxy (see e.g. Genzel et al. 1998; Laurent et al. 2000; Imanishi et al. 2006; Armus et al. 2007; Charmandaris 2008; Veilleux et al. 2009). In the present study, we focus on the modelling of a sample of ULIRGs with no indication of significant AGN contribution. The inclusion of galaxies with AGN activity will be the subject of a forthcoming study.
We select the star-forming ULIRGs to be studied here from the primary sample of Desai et al. (2007) . Our selection process of the 16 systems is not designed to be complete. Following the motivations described above we have identified galaxies for which none of the well-known optical, near-infrared and mid-infrared spectroscopic diagnostics show indication of an AGN dominating in the optical or mid-infrared. Our 16 galaxies have strong polycyclic aromatic hydrocarbon (PAH) emission, in particular at 6.2 m, and no [NeV] emission, or a mid-IR continuum indicative of an AGN. Additionally, all the ULIRGs in our sample can be classified as ‘cool’ according to their IRAS 25- to 60-m flux ratios (i.e. , see Table 1), which has been taken as an indication of no significant AGN contribution to the mid-IR emission (e.g. Sanders et al. 1988; Farrah et al. 2007). We note that even though two systems (IRAS12112+0305 and IRAS16334+4630) are classified as LINERs from optical spectroscopy in the literature (Kim et al. 1998; Veilleux et al. 2002), this does not necessarily imply the presence of an AGN but rather evidence of starburst-driven shocks (see Kim et al. 1998; Lutz et al. 1999). Furthermore, detailed analysis of the infrared emission of these galaxies using the method of Marshall et al. (2007) shows that the contribution of an AGN to their infrared luminosity would be negligible (in the case of IRAS12112+0305, this is further supported by the SED analysis of Farrah et al. 2003, who estimate that an AGN would contribute by at most 0.1% to the total infrared luminosity of this galaxy).
In Fig. 1, we compare the redshift distribution of the primary sample of Desai et al. (2007) with that of our sub-sample of 16 star-forming ULIRGs. We note that our galaxies are at typically lower redshift than the bulk of the galaxies from the primary sample due to the availability of higher signal-to-noise data for nearby systems which allowed us to select them as purely star-forming. The basic properties of our ULIRGs are listed in Table 1.
|Galaxy||R.A. (J2000)||Decl. (J2000)||log()|
2.2 Mid-infrared Spitzer/IRS spectroscopy
All ULIRGs in the sample of Desai et al. (2007) were observed in staring mode with both subslits of the Short-Low (SL) and Long-Low (LL) modules of the IRS. The final spectra have a spectral resolution typically between 60 and 120 over the wavelength range from 5 to 38 m. More detail on the IRS observations and data reduction can be found in Desai et al. (2007).
In Fig. 2, we show the mid-infrared IRS spectra of the 16 ULIRGs studied in this paper. The well-known PAH emission features, as well as silicate absorption, molecular hydrogen lines and neon fine structure lines are clearly visible. It is interesting to compare the mid-infrared spectra of our ULIRGs with those of typical local starbursts of lower infrared luminosity. Therefore, we also plot the average starburst IRS spectrum of Brandl et al. (2006) in Fig. 2. In general, the shapes of the mid-IR emission from our star-forming ULIRGs and the starburst template are remarkably similar, except for two striking differences: (i) most of our ULIRGs present much higher silicate absorption, meaning higher dust optical depths; (ii) the rising continuum from hot dust at m is steeper for the ULIRGs, indicating that dust is heated to higher temperatures (see Desai et al. 2007 for a more detailed analysis of the mid-infrared spectra of ULIRGs).
Given the small spatial scales in which star formation takes place, our ULIRGs are essentially unresolved (e.g. Soifer et al. 2001, Díaz-Santos et al. 2010), and therefore our IRS spectra are representative of the total emission of the galaxies, and ca be used to characterize the global physical properties of our galaxies. We note that we do not aim at fitting the detailed shape of the mid-infrared spectra with the model used in this paper, but rather the global SEDs from the ultraviolet to the far-infrared. We use the IRS spectra to help select star-forming galaxies using mid-IR AGN diagnostics, and also to extract essential information such as the depth of the silicate feature at 9.7 m (Section 3.1.2) and mid-infrared broad-band colours (Section 2.3) to provide general information about the shape of the mid-infrared continuum of the galaxies.
2.3 Ancillary multi-wavelength photometry
We compile available photometric data in the ultraviolet, optical and infrared in order to obtain a sampling of the observed spectral energy distribution of each galaxy as complete as possible, which is crucial to fully exploit the potential of the da Cunha et al. (2008) model. We include observations obtained with GALEX, SDSS, 2MASS, HST, ISO, IRAS and Spitzer found in the NASA/IPAC Extragalactic Database (NED). Optical fluxes for IRAS17208-0014, IRAS19297-0406 and IRAS22491-1808 were kindly provided ahead of publication by D. Sanders and V. U. For consistency, all the selected fluxes correspond to the total flux from the galaxies in each band. Additionally, we compute mid-infrared fluxes in the IRAS 12-m and 25-m, ISO 15-m, Spitzer/IRAC 5.8-m and 8.0-m, and MIPS fluxes 24-m bands by convolving the observed IRS spectra with the respective filter response functions.
Several of our ULIRGs appear double in the optical and near-IR imaging (e.g. Veilleux et al. 2002). However, in all cases we have Spitzer IRAC and/or IRS 16-m images which allow us to identify multiple mid-infrared emission components more clearly. These images reveal that, in the case of double systems, either the two components are too close and both are included within the 3.6 arcsec width of the IRS slits, or one of the two clearly dominates, contributing with over 90% of the mid-infrared flux. We have checked that in these cases the IRS spectra are taken on the mid-infrared dominant source, which continues to remain unresolved at longer infrared wavelengths. In addition, for consistency, we use ultraviolet to near-infrared fluxes extracted using large apertures which include the total flux from both sources in the case of double systems.
In Appendix A, we list all the fluxes used in this work from the ultraviolet to the far-infrared.
3.1 Spectral energy distribution modelling
We use the simple, physically motivated model of da Cunha et al. (2008) to interpret the mid- and far-infrared spectral energy distributions of galaxies consistently with the emission at ultraviolet, optical and near-infrared wavelengths. This model was previously calibrated using samples of star-forming galaxies with moderate star formation rates and infrared luminosities. Here we investigate how well this prescription can account for the SEDs of local star-forming ULIRGs, which have much higher infrared-to-optical ratios. In this section, we briefly recall the main features of the model of da Cunha et al. (2008), and describe the modifications we have performed to the original model in order to make it more suitable for the interpretation of the multi-wavelength emission by ULIRGs.
We compute the emission by stars in galaxies using the latest version of the Bruzual & Charlot (2003) population synthesis code (Charlot & Bruzual, in preparation). This code predicts the spectral evolution of stellar populations in galaxies from far-ultraviolet to far-infrared wavelengths and at ages between and yr, for different metallicities, initial mass functions (IMFs) and star formation histories. We adopt the Chabrier (2003) Galactic-disc IMF.
The ultraviolet, optical and near-infrared emission from stars is attenuated by dust. We account for this using the two-component model of Charlot & Fall (2000). This model includes the fact that stars are born in dense molecular clouds with finite lifetimes ; after that time, these birth clouds dissipate as a result of the pressure caused mainly by strong stellar winds and supernovae explosions, and stars migrate to the ambient (diffuse) ISM. Thus, the light produced by stars younger than is attenuated by dust in the birth clouds and in the ambient ISM, while the light produced by older stars is attenuated only by dust in the ambient ISM. Charlot & Fall (2000) have calibrated this model using a sample of local, ultraviolet-selected starburst galaxies of low infrared luminosities, and show that a birth cloud lifetime of yr can account for the observed line and continuum emission properties of normal star-forming galaxies.
Observations of local ULIRGs show that these galaxies are forming stars at high rates in a very concentrated central region with typical scale of at most a few kpc (e.g. CO interferometer observations by Bryant & Scoville 1996, Downes & Solomon 1998; also observations of Arp220 by Sakamoto et al. 2008). In this scenario, the central region of the ULIRG would be similar to a huge, optically-thick molecular cloud with a lifetime typically larger than yrs. Further evidence by Lahuis et al. (2007) suggests the presence of high abundances of warm, dense gas, associated with deeply embedded star formation in ULIRGs. In this environment, H ii regions are prevented from expanding by large pressure gradients of gravity, thus increasing the lifetime of the star formation process. Therefore, in the framework of the attenuation model of Charlot & Fall (2000), we adopt yr for ULIRGs. The interstellar medium is still described using two components. The first is the ‘birth cloud’ component, which accounts for the central, optically-thick starburst powering most of the infrared emission from these galaxies. The second component, the ‘diffuse ISM’, has smaller effective optical depths and is heated by stars older than . This component is necessary to account for the small but still present ultraviolet and optical emission from ULIRGs, which cannot arise from the central starburst due to the extremely high optical depths (see also spectroscopic evidence supporting this scenario in Soto & Martin 2010).
We use an ‘effective absorption’ curve for each component, , with the slope reflecting both the optical properties and the spatial distribution of the dust. Following Charlot & Fall (2000), we adopt for the diffuse ISM
where is the total effective -band absorption optical depth of the dust seen by young stars inside birth clouds, and is the fraction of this contributed by dust in the surrounding diffuse ISM. For the stellar birth clouds, we adopt:
The slopes of the effective absorption curve in the birth clouds and in the diffuse ISM are different in order to account for the different spatial distributions of dust. For the birth cloud component, we adopt , which corresponds to a foreground dust screen. This should be adequate to characterize the attenuation of optically-thick dust in the concentrated dense central starburst region of ULIRGs. The diffuse ISM is better described by a random distribution of discrete dust clouds (Charlot & Fall 2000), which can be accounted for adopting (i.e. a ‘greyer’ effective absorption curve).
We use the prescription described above to compute the total energy absorbed by dust in the ‘birth cloud’ component (i.e. the dense central starburst) and in the surrounding ‘ambient ISM’; this energy is then re-radiated by dust at infrared wavelengths (assuming conservation of energy). The total luminosity emitted by dust in the galaxy is then
where and are the total luminosity re-radiated by dust in the birth cloud and in the ambient ISM components, respectively.
This approach, developed in Charlot & Fall (2000) and da Cunha et al. (2008), allows us to interpret the infrared emission from galaxies consistently with the emission at shorter wavelengths using a simple energy balance argument (given that the main energy source in galaxies are stars, i.e. there is no contribution by an AGN). For normal star-forming galaxies with moderate amounts of dust (optically-thin case), the effective -band optical depth seen by stars younger than inside birth cloud is constrained by the attenuation suffered by emission lines, e.g. , while the effective -band optical depth seen by older stars in the diffuse ISM is mainly constrained by the attenuated ultraviolet and optical continuum emission from the galaxy.
In starburst-dominated ULIRGs, practically all the radiation from stars in the birth cloud component is absorbed by the optically-thick dust, and therefore it is not possible to reliably constrain using observations in the ultraviolet or optical. In such cases, where the optical depths are very high, dust will absorb effectively even mid-infrared radiation, in particular in the silicate band located around 9.7 m, where the grain cross-section is maximum (and also to a lesser extent at 18 m). The Spitzer/IRS mid-infrared spectra of our ULIRGs allow us to measure the strength of the 9.7 m silicate feature for each one of them. This feature strength can be converted to an apparent dust optical depth , which can be used to compute the -band optical depth by assuming dust optical properties and geometry. With this in mind, for each galaxy, we compute the strength of the 9.7-m silicate absorption feature, , from the observed IRS spectra, as:
where is the observed flux density of the feature and that of the continuum, both evaluated at the wavelength, , where the absorption presents its maximum (typically about 9.7 m). To compute the we fit the IRS spectrum of each galaxy to a simple linear function with anchors at 5.5 and 13.2 m (rest-frame) and interpolate it through the absorption feature. This method is similar to that used by Sirocky et al. (2008) for measuring the silicate strength in PAH-dominated IRS spectra, but here we use a long-wavelength anchor at 13.2 m instead of at 14.5 m and a straight line instead of a power-law to fit the spectra. We note that both methods yield similar results for . We have checked that our values of the silicate absorption strength, shown in Table 2, are consistent with those found by Hao et al. (2007) for a larger sample of local ULIRGs and are significantly larger than the silicate absorption strength of local starbursts of lower infrared luminosity of Brandl et al. (2006), as expected from Fig. 2.
We assume that the observed silicate feature of our ULIRGs originates in the birth cloud component (which can be approximated by an optically-thick dust screen), and use the ULIRG silicate model of Sirocky et al. (2008) to estimate the -band optical depth in the birth clouds from the apparent 9.7-m optical depth inferred from the silicate absorption feature strength. The adopted model assumes a Mathis et al. (1977) grain size distribution and a cold Ossenkopf et al. (1992) dust screen, which imply , where (see more details in Sirocky et al. 2008). We include a 30% uncertainty associated with derived in this way to account for uncertainties related to the measurement of (about 10%) and model assumptions. Table 2 shows the silicate absorption feature strengths of our ULIRGs, and the -band optical depth predicted using the model of Sirocky et al. (2008). We will use these values to help constrain our model , as described in Section 3.2.2.
We distribute and in wavelength between 3 and 1000 m using four main components (see da Cunha et al. 2008 for details): (i) emission from polycyclic aromatic hydrocarbons (PAHs; i.e. mid-infrared emission features), (ii) mid-infrared continuum emission from hot dust with temperatures in the range 130–250 K, (iii) emission from warm dust in thermal equilibrium with adjustable temperature, and (iv) emission from cold dust in thermal equilibrium with adjustable temperature.
In stellar birth clouds, the relative contributions to by PAHs, the hot mid-infrared continuum and warm dust are kept as adjustable parameters. These clouds are assumed not to contain any cold dust. In the ambient ISM, the contribution to by cold dust is kept as an adjustable parameter. The relative proportions of the other three components are fixed to the values reproducing the mid-infrared cirrus emission of the Milky Way (more detail in section 2.2 of da Cunha et al. 2008).
3.2 Median-likelihood estimates of physical parameters
We use the model described in the previous section to interpret the observed spectral energy distributions of our ULIRGs and derive median-likelihood estimates of their physical parameters. To do so, we adopt a Bayesian approach similar to that used by da Cunha et al. 2008 (see also Kauffmann et al. 2003, Brinchmann et al. 2004, da Cunha et al. 2010). This method relies on a large library of random models encompassing all plausible parameter combinations (e.g. star formation histories, metallicities, dust optical depths, dust masses, dust temperatures). Such a library was already built in da Cunha et al. (2008, 2010) to derive statistical constraints on the physical parameters of local star-forming galaxies of moderate infrared luminosities. For each galaxy, we build the likelihood distribution of any physical parameter by evaluating how well each model in the library accounts for the observed SED of the galaxy. The underlying assumption of our method is that the model library represents the distribution from which the data were randomly drawn, so it is important that the prior distributions of the parameters sample well the observational space and do not give too much weight to a priori implausible corners of the parameter space.
To reproduce the observations of our ULIRGs, we find that some of the parameter priors have to be changed with respect to the priors used to model galaxies with in da Cunha et al. (2008, 2010). This confirms the idea that ULIRGs are probing a very different region of the parameter space of galaxies. In particular, we include higher dust optical depths and temperatures and increase the probability and strengths of bursts of star formation. We describe the priors in detail in the next section.
We first build a large library of 50 000 stellar population models with random star formation histories, metallicities and dust contents. We distribute the galaxy age uniformly between 2 and 13.5 Gyr. We parametrize the star formation history of each model in terms of two components: an exponentially declining star formation rate of the form , where we distribute the time-scale parameter between 0 and 1 in the same way as da Cunha et al. (2008), and superimpose random bursts. Since our ULIRGs are highly star-forming galaxies, likely to be experiencing star formation bursts, we assign a 75 per cent probability of a burst occurring in the past 2 Gyr. The amplitude of each burst is defined as the ratio between the mass of stars formed during the burst and the total mass of stars formed by the underlying continuous model, . We distribute logarithmically between 0.1 and 10. The duration of each burst, , is distributed uniformly between and yr. The metallicity is uniformly distributed between 0.2 and 2 times solar.
To account for the high optical depths of ULIRGs, we have to change significantly the priors of the attenuation parameters (effective -band optical depth seen by stars younger than in the birth cloud component) and (fraction of contributed by the diffuse ISM component), compared to the previous study of da Cunha et al. (2008). This is motivated by the very high infrared-to-optical ratios observed in the galaxies for which optical observations are available, and also by the high -band optical depths derived from the IRS spectra using the silicate feature (Table 2). The strength of the silicate feature in our ULIRGs implies that the birth cloud component has to be optically-thick, and the fact that ultraviolet and/or optical emission is observed for several of our galaxies means that the diffuse ISM component should have much smaller optical depths. Therefore, we distribute uniformly in log between 2 and 50. These high optical depths imply that the optical emission from ULIRGs comes from a more diffuse, optically-thin component – the diffuse ISM component. An optically-thin ISM component is the only way the optical fluxes can be fitted consistently without the diffuse ISM contributing too much to the total infrared luminosity (i.e. too high ). This is the motivation for drawing the -band optical depth in the diffuse ISM, from a Gaussian prior centered at 1 with 0.5 dispersion, between 0 and 2. We note that this is in agreement with studies of the extinction in spatially-resolved ULIRGs (García-Marín et al. 2009), where values of optical extinction as low as 0.2 mag are found for the external regions (in typical scales of 1 – 3 kpc).
We also generate a library of 50 000 infrared emission models in a similar way to da Cunha et al. (2008). We distribute the fraction of total dust luminosity contributed by the diffuse ISM, , uniformly between 0 and 1. We randomly draw the fractional contributions by PAHs, hot mid-infrared continuum and warm dust in thermal equilibrium to the total infrared luminosity of birth clouds, from the same prior distributions as da Cunha et al. (2008). The equilibrium temperature of warm dust in birth clouds, , is distributed uniformly between 30 and 60 K. For the equilibrium temperature of cold dust in the diffuse ISM (), we extend the prior of da Cunha et al. (2008) to include higher values in order to account for the observed hotter far-infrared colours of our ULIRGs; is uniformly distributed between 15 and 30 K.
Finally, we combine the library of stellar population models and the library of dust emission models by joining models with similar parameter and scaling to the total infrared luminosity as detailed in da Cunha et al. (2008). For each model spectrum in our library, we compute the synthetic photometry in photometric bands from the ultraviolet to the far-infrared at the redshifts of our ULIRGs.
We perform spectral fits by comparing the observed spectral energy distributions and also the -band optical depths of the birth cloud component derived from the silicate strengths measured using the Spitzer/IRS spectra in Section 3.1.2 to every model in our stochastic library. Specifically, for each ULIRG, we compute the goodness-of-fit of each model in the library:
where and are the luminosities in the band of the observed galaxy and the model, respectively, is the observational uncertainty in , is the model scaling factor that minimises [see eq. (33) of da Cunha et al. (2008)], is the -band optical depth derived from the observed 9.7 m silicate feature and is the uncertainty associated with this measure. As discussed in Section 3.1.2, we explicitly fit to the effective -band optical depth in the birth clouds because they are so optically-thick this information is not accessible from ultraviolet and optical observations.
We fit all the available observed broad-band fluxes of our galaxies from the ultraviolet to the far-infrared, plus the -band optical depth of the birth clouds inferred from the mid-infrared silicate absorption strength (eq. 5). We do not attempt at reproducing the detailed shape of the observed IRS spectra of our galaxies because our SED model is not designed to reproduce the detailed emission in this spectral range (so that the number of free parameters is kept minimal). However, we aim at reproducing the rough shape of the mid-infrared SED of our galaxies, since this reflects the emission by hot, stochastically-heated dust, and also by PAHs. To do so, we compute broad-band fluxes in the mid-infrared from the observed IRS spectra and include them as observables in the fitting procedure.
We build the likelihood distribution of any given physical parameter for our observed galaxies by weighting the value of that parameter by the probability . Our best estimate of the parameter is the median of the resulting probability density function and our confidence interval the 16th–84th percentile range. The results of our SED fitting for the sample of 16 star-forming ULIRGs using this method are presented in the next section.
4 UV-to-IR SEDs and the physical parameters of star-forming ULIRGs
In Fig. 3, we show the results of our SED fitting for each of our 16 ULIRGs. We compare, in the top panels, the best-fit models (in black) to the observed multi-wavelength SEDs (in red). The residuals are shown in the bottom of each spectrum. The blue line shows the unattenuated emission from stars. For all our galaxies, this emission is very strong and blue, indicating that the stellar component is dominated by young stellar populations resulting from current or very recent star formation. Most of this stellar radiation is absorbed by dust and re-emitted in the infrared, leading to the large observed infrared-to-optical ratios in these galaxies. The infrared emission of our galaxies peaks typically at wavelengths shorter than 100 m, indicating overall high dust temperatures. For each galaxy, as a reference, we also plot the observed Spitzer/IRS spectrum (green line). We note that, even if we did not intend to fit the detailed shape of the mid-IR emission, our best-fit SED and the IRS spectrum are in good agreement in most cases. Both the shape of the PAH emission and the slope of the hot mid-infrared continuum are reproduced remarkably well considering the simplicity of the model. However, the fits do not reproduce the observed silicate absorption feature because dust self-absorption is not included in our infrared SEDs. This should have a negligible contribution in terms of the energy balance performed here.
Even if we do not explicitly include the silicate absorption in our infrared SEDs, the silicate optical depth measured from the observed IRS spectra provides valuable information about the optical depth, as discussed in Section 3. We find that in heavily dust-enshrouded systems as our ULIRGs, constraining well the dust optical depth is fundamental in our SED fitting and has a high impact in the statistical constraints of star formation rates, stellar masses, and contribution by the diffuse ISM and the birth clouds components to the dust heating. To illustrate how well the -band optical depth in the birth clouds is recovered by our best-fit model, in Fig. 4, we plot the distribution of the difference between the observational () and the best-fit () -band optical depth for the galaxies in our sample. In most cases, the -band optical depth is recovered within (i.e. 60% of ). This is a satisfactory result considering the uncertainties in deriving from the observed IRS spectra.
An advantage of our method is that it enables us not only to find a best fit model, but also to compute the full likelihood distribution of each model parameter and identify possible degeneracies between parameters. This allows us to determine how well the parameters are constrained and how the available set of observations affects these constraints (see more details in section 3.2 of da Cunha et al. 2008). The bottom panels of Fig. 3 show the likelihood distribution of several physical properties of our ULIRGs derived from the comparison of the observed SEDs with our stochastic library of models as described in Section 3.2.2. We list the median-likelihood estimates and confidence intervals of some important parameters for each galaxy in Table 3. The likelihood distributions of Fig. 3 show that all the parameters are reasonably well constrained for most of our galaxies. However, we can see some differences in how well some parameters are constrained from galaxy to galaxy, and this is primarily due to the different set of available photometry for each galaxy. Clearly, the more available data points, the better constrained are the parameters.
It is important to note how essential it is to include ultraviolet to near-infrared observations, even if the emission from ULIRGs is largely dominated by longer wavelengths. UV to near-IR observations probe stellar populations older than yrs which heat the diffuse ISM (i.e. stellar population which do not belong to the central, heavily obscured starburst). These observations help constrain the relative amount of dust luminosity heated by these stars (), the dust attenuation in the diffuse ISM (), the stellar mass () and, as a consequence, the specific star formation rate (). This is visible in the case of IRASF10156+3705, for which no observations short-ward of 3.6 m are available, leading to typically more extended likelihood distributions (and hence wider confidence ranges) for these parameters. We note, however, that even in this case, the star formation rate is still relatively well-constrained because, given the input optical depth derived from the IRS spectrum, the model determines what is the necessary amount of star formation to power the very large infrared luminosity of the galaxy (constrained by the infrared observations).
In the infrared, given the wealth of available data at least up to 160 m, the total dust luminosity () and the contribution to the total emission by PAHs (), the hot mid-IR continuum (), and warm dust in thermal equilibrium (), as well as the temperature of this dust (), are well constrained in most of the cases. Cold dust emits in the far-infrared to sub-millimetre, and hence is constrained by the available observations from 160 to 850 m. In general, the contribution of this component to the total emission () is well constrained by these observations. However, the equilibrium temperature of the cold dust component in our model () is relatively hard to constrain for our ULIRGs even with far-infrared data, because the contribution of this dust to the total infrared emission is small and hard to isolate from the dominating warm dust component. We note that this contributes to the broadening of the likelihood distributions of the total dust mass, since cold dust, even while contributing little to the total infrared luminosity, contributes significantly to the dust mass.
We have compared our median-likelihood estimates of the total dust luminosity with the infrared luminosities empirically derived using the IRAS fluxes, , presented in Table 1. The two luminosities agree within typically 0.07 dex (with no systematic offset) for most of the galaxies. For only four galaxies, the difference between the two luminosities is bigger: IRASF08208+3211, IRASF10156+3705, IRAS17208-0014 and IRAS19458+0944, although even in these cases the maximum offset between the two luminosities is at most 0.14 dex (i.e. 38%).
It is important to check how our stellar mass estimates compare with the results of previous studies, specially considering the very large optical depths of ULIRGs. Rodríguez Zaurín et al. (2010) derived the stellar masses of 36 local ULIRGs using optical spectroscopy and found that ULIRGs have moderate stellar masses (using a Salpeter 1955 IMF and Bruzual & Charlot 2003 models spectral synthesis models; see also, e.g. Genzel et al. 2001; Tacconi et al. 2002). We note that the use of a Salpeter IMF and the Bruzual & Charlot (2003) models can make the estimates systematically higher by 0.15–0.3 dex compared to the use of a Chabrier IMF and the latest version of the Bruzual & Charlot models used in our study (Bruzual 2007; da Cunha et al. 2010). Nevertheless, we find that our estimates of the stellar masses of our ULIRGs are consistent with the previous result of Rodríguez Zaurín et al. (2010), except for IRASF10156+3705, for which we find a significantly larger , but this may be attributed to the more uncertain estimate of for this galaxy due to the lack of optical and near-infrared data (as mentioned above).
5.1 Comparison between star formation rate and infrared luminosity
To test the validity of our approach, we compare our star formation rate estimates with the ones obtained from the infrared luminosity using the formula of Kennicutt (1998). This formula assumes that all the infrared emission from the galaxy is powered by recent star formation, i.e. all the energy emitted by young stars in the galaxy is absorbed and re-radiated by dust, and no dust is heated by old stellar populations (optically-thick starburst). In our model, this corresponds to having very small values of the parameter, i.e. the contribution by the diffuse ISM (heated by stars older than 100 Myr) to the total dust luminosity is very small (see Table 3). Also, to power the huge infrared luminosities of ULIRGs with young stars alone, the effective attenuation in the birth clouds must be very high, as well as the star formation rates. Therefore we expect the star formation rates of our ULIRGs to be similar to the ones derived using the Kennicutt formula.
In Fig. 5 and Table 4, we compare our median-likelihood estimates of the star formation rates of our ULIRGs with the star formation rates predicted from the total infrared luminosity using the Kennicutt (1998) formula for optically-thick starbursts (corrected for a Chabrier IMF). The star formation rate is computed as:
where is the total infrared luminosity between 8 and 1000 m, which we approximate with the total infrared luminosity derived from our fits, .
We verify that our star formation rates agree remarkably well with the ones predicted using the infrared luminosities. The difference between and is typically less than (except for IRASZ02376-0054 and IRAS17208-0014, but even in those cases is still in the confidence range of ) . This confirms that the model is correctly attributing most of the infrared emission of the ULIRGs to the starburst phase (i.e. the ‘birth cloud’ component in our model). This was only achieved through a correct estimation of the optical depths enabled by the use of the IRS spectra.
We caution that in more quiescent galaxies, dust heated by old stars in the ambient ISM contributes more significantly to the total infrared luminosity, and the star formation rates obtained using da Cunha et al. (2008) may differ significantly than those derived using the formula of Kennicutt (1998). In such cases the star formation rates would be overestimated using the Kennicutt law.
5.2 The star formation histories of ULIRGs
The model described in Section 3.1 allows us to explore possible star formation histories of galaxies, and the effect they have on observable properties. One interesting result of the SED fitting of our ULIRGs is that, for all 16 of them, the star formation history corresponding to the best-fit model (Fig. 3) shows a strong burst of star formation occurring at the present age, i.e. all the best-fit models are clearly ones with a current starburst (see example of Fig. 6). We find that, on average, this burst of star formation is responsible for forming of the measured stellar mass of our ULIRGs over the last 100 Myr (this corresponds to about formed during the burst).
The observed spectral energy distributions of our ULIRGs can only be well fitted with a simple, exponentially-declining star formation history with superimposed random bursts (as we use in our stochastic library of models; see Section 3.2.1) if the star formation history contains a strong burst of star formation at the present age. This star formation history, while not necessarily being a unique solution, is a plausible scenario consistent with the main result that can be inferred from the SED fitting – that the galaxy must be experiencing a strong burst of star formation where a significant fraction of its stellar mass is being formed. To constrain the detailed star formation history of a galaxy, we would require optical spectroscopic data to characterize the stellar populations of different ages (e.g. Rodríguez Zaurín et al. 2008, 2010). Nevertheless, our results show that fitting the broad-band SEDs can at least provide some clues on the star formation mode of these galaxies. In particular, our fits clearly show that, in order to reproduce the observed properties of our ULIRGs, a star formation history with a strong current burst of star formation is required.
5.3 Comparison with local galaxies of infrared luminosities between and
One of the main advantages of the method used in this paper is that our results are directly comparable to results found for other samples of galaxies using the same framework. The physical properties of a sample of 3258 local star-forming SDSS galaxies () detected with IRAS, for which GALEX and 2MASS observations are also available (SDSS-IRAS sample from now on), were derived in da Cunha et al. (2010) using the technique used this paper. In this section, we compare the physical properties of local galaxies in two bins of infrared luminosity: (1538 galaxies with high S/N observations selected from the SDSS-IRAS sample of da Cunha et al. 2010), and (15 ULIRGs from the present study). This is done in order to examine whether our findings based on SED modelling are consistent with the current understanding of ULIRGs. In Fig. 7, we plot the distributions of the median-likelihood estimates of several physical parameters related to the star formation activity and the dust properties of the galaxies. We describe the main similarities and differences between the two samples, even considering the low number of ULIRGs compared to the number of lower-luminosity galaxies.
Fig. 7(b) shows that the stellar masses of galaxies in the two infrared luminosity bins are similar, with the ULIRGs tending to be slightly more massive. We note, however, that the stellar mass estimates of ULIRGs are typically more uncertain than those of galaxies of lower infrared luminosities (and lower dust attenuation). In Fig. 7(c) we show that ULIRGs have typically much higher specific star formation rates than the SDSS-IRAS galaxies, implying doubling timescales of their stellar mass (assuming that they continue forming stars at the present rate) of the order of less than 1 Gyr [Fig. 7(d)]. This is another confirmation of the well-known starbursting nature of our local ULIRGs (see, e.g. Genzel et al. 1998).
We find that galaxies with also present very different properties of the dust emission than galaxies with . Fig. 7(f) shows that, in ULIRGs, the contribution of the diffuse ISM component to the total infrared luminosity, , is typically very small (about 0.10), compared with what is found for the SDSS-IRAS galaxies. da Cunha et al. (2010) find that, for these galaxies, dust heated by old stars in the diffuse ISM contributes significantly to the total infrared luminosity, with most values of typically between 0.40 and 0.70.
Furthermore, Figs. 7(g) and (i) show that the typical equilibrium temperatures of dust grains in the ‘birth cloud’ component, , and in the ‘diffuse ISM’ component, , are typically higher in ULIRGs than in the lower infrared luminosity galaxies. Additionally, the contribution of warm dust to the total infrared luminosity is very large for ULIRGs [Fig. 7(h)], while the contribution of cold dust is much smaller than in the other sample [Fig. 7(j)]. Warm dust is the dominant contributor to the total dust emission in ULIRGs, and therefore, the overall temperature of dust in our star-forming ULIRGs is higher than that of the SDSS-IRAS galaxies (see also Clements et al. 2010).
Based on the model fits we may also explore how the relations between the star formation rate, stellar mass and dust mass can be useful to study the star formation mode and evolution of our galaxies. In Fig. 8(a), we plot our median-likelihood estimates of the star formation rate against the stellar mass for our 15 local ULIRGs with , and the 1538 SDSS-IRAS galaxies with . The SDSS-IRAS galaxies present a strong correlation between and (the Spearman test indicates a 18- significance level for this correlation). Previous studies of local and intermediate-redshift galaxies have found that star-forming galaxies exhibit a strong correlation between stellar mass and star formation rate, forming a tight ‘main sequence’ in which the range of possible star formation rates of a galaxy is defined by its stellar mass (Brinchmann et al. 2004; Noeske et al. 2007b; Elbaz et al. 2007). A similar correlation has also been found for galaxies at (Elbaz et al. 2007) and at (Daddi et al. 2007), with the overall normalization of the correlation increasing with redshift, i.e. typical galaxies present higher star formation rates for a given stellar mass at high redshifts. Noeske et al. (2007a) explain this correlation in terms of a mass-dependent gradual decline of the star formation in disk galaxies as a result of gas exhaustion as they evolve. It is interesting to analyze the distribution of local ULIRGs in the star formation vs. stellar mass plot. The ULIRGs are clearly above the relation defined by the SDSS-IRAS galaxies, i.e. at fixed stellar mass, ULIRGs have higher star formation rates (by about one order of magnitude). Additionally, our ULIRGs do not show any statistically significant correlation between and .
In Fig. 8(b), we plot the median-likelihood estimates of star formation rate against dust mass for the two samples. da Cunha et al. (2010) find a very tight correlation between and for the SDSS-IRAS galaxies. Interestingly, our ULIRGs do not follow this relation. For a given dust mass, the ULIRGs have higher star formation rates than the SDSS galaxies – they are forming stars about 10 times more efficiently for a given . If we assume a constant dust-to-gas ratio and rely on the Schmidt-Kennicutt law (Schmidt 1959; Kennicutt 1998), this implies denser star-forming gas in ULIRGs (as observed by e.g. Gao & Solomon 2004; Bryant & Scoville 1999; Iono et al. 2009; Juneau et al. 2009). These results are consistent with the scenario in which gas and dust in local ULIRGs are concentrated in a compact central region (e.g. Soifer et al. 2001; Díaz-Santos et al., 2010) following a merger of two gas-rich galaxies, and creating a strong starburst that powers the huge infrared luminosities of these systems (as found by simulations of major mergers of spirals; e.g. Mihos & Hernquist 1994; Mihos 1999; Hopkins et al. 2006).
While these findings are not new, we show here how our SED modelling can be useful to derive relevant physical parameters which provide clues to the star formation mode of galaxies by using simple photometric observations. Additionally to a higher star formation efficiency, we also show dust in ULIRGs is typically warmer, as shown in the histograms of Fig. 7. This is again consistent with the idea that the star formation in these galaxies is more spatially concentrated, and dust grains remain closer to the young, newly-formed stars, reaching higher equilibrium temperatures (see also, e.g. Yang et al. 2007; Groves et al. 2008).
6 Summary and Conclusions
In this paper, we have studied the physical properties of a sample of 16 local star-forming ULIRGs observed in the mid-infrared with Spitzer/IRS for which we have compiled additional photometry from the ultraviolet to the infrared.
We have interpreted the observed ultraviolet-to-infrared spectral energy distributions of these galaxies using the simple and physically-motivated model of da Cunha et al. (2008). Using this prescription, we describe our ULIRGs using two main components: (i) an ‘optically-thick starburst’ with a lifetime of 100 Myr, analogous to the ‘birth cloud’ component in the original model of da Cunha et al. (2008), where dust is heated by young stars; (ii) a diffuse ISM component, with smaller dust effective optical depths, where dust is heated by stars older than 100 Myr. The model accounts for the emission of stellar populations from the ultraviolet to the near-infrared and for the attenuation by dust in the optically-thick starburst component and in the diffuse ISM. The infrared emission is simply computed by assuming that all the energy absorbed by dust in these components is re-radiated at mid- and far-infrared wavelengths, using four main dust components. We have combined this model with a Bayesian technique to derive statistical constraints on parameters such as the star formation rate, stellar mass, dust mass, dust temperature, and contribution of each dust component to the total infrared luminosity from their observed multi-wavelength spectral energy distributions.
We show that the model of da Cunha et al. (2008) is versatile enough to successfully account for the observed spectral energy distributions of star-forming ULIRGs. Although the ultraviolet to near-infrared emission represents only a small fraction of the total power radiated by ULIRGs, observations in this wavelength range are important to constrain the dust attenuation and the fraction of infrared luminosity contributed by the diffuse ISM, the total stellar mass, and the specific star formation rate. Additionally, our method can also provide clues to the star formation mode of our galaxies. The best fit models to the spectral energy distributions of our ULIRGs indicate that the observations can only be optimally matched if the star formation history includes a strong burst of star formation at the present age. While this is not surprising, given the known starburst nature of these galaxies, it demonstrates the consistency of our method.
To put the physical properties of local star-forming ULIRGs in context, we have compared them with a sample of SDSS galaxies of lower infrared luminosities () detected in the infrared with IRAS (da Cunha et al. 2010). Our results are consistent with a well-known fundamental difference in the star formation mode of local ULIRGs and other local galaxies of lower infrared luminosities. ULIRGs have star formation rates and dust temperatures significantly higher than the other galaxies, consistent with a merger-induced starburst, in which gas is driven to the central region, forming stars at very high rates and thus powering the huge infrared luminosities of these galaxies.
The method used in this paper provides a framework to better understand local star-forming galaxies spanning a wide range of infrared luminosities in a statistical manner. This has the additional advantage that it can also be used to preform quantitative comparisons to high-redshift galaxies, for which often much more sparsely sampled SEDs are available.
Acknowledgements.We thank the anonymous referee for comments and suggestions that significantly helped improve this paper. EdC, VC and TDS acknowledge partial support from the EU ToK grant 39965 and FP7-REGPOT 206469. We thank David Sanders and Vivian U for providing us the optical photometry for some of our galaxies. We are grateful to Bernhard Brandl and Emanuele Daddi for useful discussions, and to Brent Groves, Nick Kylafis, Vivienne Wild and Stephanie Juneau for comments on the manuscript. This research has made use of the NASA/IPAC Extragalactic Database (NED) which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration.
Appendix A Multi-wavelength photometric data
Notes: – 2MASS Point Source Catalog; – computed from convolving the observed IRS spectra with the corresponding IRAC bandpass; – ; – Scoville et al. (2000); – 2MASS XSC.