Planck Observations of M33
C. T. Tibbs F. P. Israel R. J. Laureijs J. A. Tauber B. Partridge M. W. Peel L. Fauvet




firstpage–lastpage 

 Improved Silbey-Harris polaron ansatz for the   spin-boson model
    Shu He^1,2, Liwei Duan^1, and Qing-Hu Chen^1,3,*
    
=================================================================









We have performed a comprehensive investigation of the global integrated flux density of M33 from radio to ultraviolet wavelengths, finding that the data between ∼100 GHz and 3 THz are accurately described by a single modified blackbody curve with a dust temperature of T_dust = 21.67 ± 0.30 K and an effective dust emissivity index of β_eff = 1.35 ± 0.10, with no indication of an excess of emission at millimeter/sub-millimeter wavelengths. However, sub-dividing M33 into three radial annuli, we found that the global emission curve is highly degenerate with the constituent curves representing the sub-regions of M33. We also found gradients in T_dust and β_eff across the disk of M33, with both quantities decreasing with increasing radius. Comparing the M33 dust emissivity with that of other Local Group members, we find that M33 resembles the Magellanic Clouds rather than the larger galaxies, i.e., the Milky Way and M31. In the Local Group sample, we find a clear correlation between global dust emissivity and metallicity, with dust emissivity increasing with metallicity. A major aspect of this analysis is the investigation into the impact of fluctuations in the Cosmic Microwave Background (CMB) on the integrated flux density spectrum of M33. We found that failing to account for these CMB fluctuations would result in a significant over-estimate of T_dust by ∼5 K and an under-estimate of β_eff by ∼0.4. 






galaxies: individual: M33 – galaxies: ISM – galaxies: photometry – infrared: galaxies – submillimetre: galaxies – radio continuum: galaxies







§ INTRODUCTION



In the region between high-frequency radio waves (ν ≳ 10 GHz) and long-wavelength infrared (IR) emission (λ ≳ 100 μm), thermal radiation from interstellar dust and ionized gas, as well as non-thermal synchrotron radiation, all contribute to the emission from cosmic objects. By unravelling the various contributions, we may obtain information on the ionising stars and the properties of interstellar dust in a variety of galactic environments. The observations provided by the Planck mission <cit.> allow us to sample the poorly observed far-IR to millimetre (mm) gap in the continuum emission spectrum of objects such as entire galaxies. 

In the past, attempts have been made to extrapolate the incomplete IR continuum flux density spectrum (frequently, but incorrectly, referred to as a spectral energy distribution or SED[An SED is a plot of energy as a function of frequency or wavelength, i.e., ν S_ν vs ν, or λ S_λ vs λ, while a flux density spectrum is a plot of flux density as a function of frequency or wavelength i.e., S_ν vs ν, or S_λ vs λ.]) cutting off somewhere between 100 and 160 μm (IRAS, Spitzer Space Telescope) by assuming a single effective integrated dust emissivity index of β_eff = 2 for the Rayleigh-Jeans extrapolation. Often, the values actually measured at wavelengths around 1 mm significantly exceed such extrapolated flux densities. This so-called “millimetre excess” was readily interpreted as evidence for a large mass of colder dust <cit.>. However, both the increased far-IR wavelength coverage (up to 500 μm) of the Herschel Space Observatory and the results of terrestrial laboratory experiments have subsequently indicated that the actual value of β_eff is generally < 2. As a consequence, both the historic millimetre excess and the implied large mass of colder dust, can be reduced to an artefact of interpretation, and effectively disappear.

Thus far, reliable and complete continuum flux density spectra ranging from the radio to the mid-IR or even optical wavelengths, well-sampling the mm to far-IR range, have been published for a variety of Milky Way sources, but only for a few galaxies beyond. The wavelength coverage of Planck renders extrapolation superfluous; the value of β_eff can be measured directly. This measurement is complicated by the degeneracy between β_eff and the dust temperature, T_dust, derived from flux densities that have finite instrumental noise: the two parameters are inversely correlated <cit.>. Nevertheless, the Planck spectra of the Local Group galaxies, the Large and Small Magellanic Clouds <cit.>, and M31 <cit.>, imply galaxy-wide effective dust emissivities well below two. Similar emissivities have been found for other nearby galaxies <cit.>.

Surprisingly, the complete flux density spectra of the LMC and SMC, incorporating WMAP and COBE data, published by <cit.> and interpreted by <cit.>, do show a pronounced excess of emission at mm to cm wavelengths. This “new” excess emission is not to be confused with the apparent historical millimetre excess discussed earlier as it does not result from an arbitrary assumption on the dust emissivity, but is a well-sampled spectral feature. Its existence was confirmed by <cit.>, who explained the observed excess in the LMC as a fluctuation of the Cosmic Microwave Background (CMB), but admitted to the presence of a significant intrinsic excess in the SMC. <cit.> proposed that this intrinsic excess in the SMC could be explained if the interstellar dust includes magnetic nanoparticles, emitting magnetic dipole radiation resulting from the thermal fluctuations in the magnetisation.

Perhaps relatedly, an excess of emission at longer (cm) wavelengths has also been observed in many environments within the Milky Way <cit.>. This cm excess, more commonly known as anomalous microwave emission (AME), is typically observed at frequencies around 30 GHz (or wavelengths of 1 cm), is observed to be highly correlated with the IR dust emission <cit.>, and is believed to be due to electric dipole radiation from very small rapidly spinning dust grains <cit.>.

In this paper we present a study of the small Local Group spiral galaxy M33, using the most recent Planck data along with data from the literature, to produce a comprehensive continuum flux density spectrum from radio to ultraviolet (UV) wavelengths. We profit from the fact that, due to its proximity <cit.> and modest dimensions (approximately 70 arcmin × 40 arcmin - see Fig. <ref>), M33 is an exceedingly well-studied object. In this analysis we will specifically address: (a) the shape of the Rayleigh-Jeans spectrum, (b) the magnitude of the effective dust emissivity spectral index, β_eff, and (c) possible differences between the inner and outer regions of M33. Since the flux density spectra of both individual interstellar clouds and entire galaxies have a minimum close to the peak of the CMB, at these frequencies the CMB emission typically exceeds the interstellar contribution. Thus, our results depend critically on the reliability of the CMB subtraction. For this reason, we will pay special attention to an analysis of the CMB fluctuations, as these dominate the M33 spectrum at mm wavelengths. 

This paper is organised as follows. In Section <ref> we describe the data used in this analysis, while in Section <ref> we produce a global continuum flux density spectrum for M33, accounting for contributions from both CMB fluctuations and CO line emission. We also spatially decompose M33 into three annuli, producing a flux density spectrum for each. In Section <ref> we discuss the results of our work, and we present our conclusions in Section <ref>.






§ DATA












 §.§ Planck



The Planck mission <cit.> was the third cosmological satellite mission to observe the entire sky in a series of wide spectral passbands (Δν/ν∼0.3–0.6) designed to sample the CMB. It measured the emission from the sky with the Low Frequency Instrument (LFI) at 28.4, 44.1, and 70.4 GHz (1.0–0.4 cm) with amplifiers cooled to 20 K between 2009 and 2013, and the High Frequency Instrument (HFI) at 100, 143, 217, 353, 545, and 857 GHz (3.0–0.35 mm), with bolometers cooled to 0.1 K between 2009 and early 2012 (cf., Table <ref>). 

In this paper, we use the most recent Planck 2015 “full” data release <cit.>. These data cover the full mission from 12 August 2009 to 23 October 2013 and are available to download from the Planck Legacy Archive.[http://pla.esac.esa.int/pla/] The full LFI/HFI data processing and calibration procedures are described in <cit.> with an overview provided in <cit.>. The Planck full-sky maps are provided in HEALPix format <cit.>, but for this analysis we extracted 2D projected maps centred on M33 using the Gnomdrizz package <cit.>, which accurately conserves the photometry during the data re-pixelization. After performing this extraction for each of the nine Planck frequency maps, the maps were converted from units of CMB temperature to MJy sr^-1 using the coefficients described in <cit.>. The nine Planck frequency maps are displayed in the left column of Fig. <ref>.

The 2015 Planck data have been calibrated on the orbital modulation of the “cosmological dipole”, resulting in extremely high (sub-percent) calibration accuracy <cit.>. However, it is important to note that the quoted accuracies are appropriate for diffuse emission at large angular scales, where the calibration signal appears. Additional uncertainties apply at smaller angular scales. For relatively compact sources such as M33, the main additional contributors are related to colour correction and to beam uncertainty. Colour correction uncertainties depend on the spectral shape of the source <cit.>, while the beam uncertainties depend on angular scale <cit.>; for M33 we conservatively assume the entire solid angle uncertainty to be applicable. Combining these uncertainties in quadrature, we conservatively assume a photometric uncertainty of 7 % for the 545 and 857 GHz bands, and 1 % for the other seven bands. Ultimately, the uncertainty on the flux determination of compact sources is limited by fluctuations in both the physical backgrounds and foregrounds rather than photometry errors.




 §.§ Herschel



M33 was mapped with Herschel within the framework of the open-time key programme as part of the HerM33es KPOT_ckrame01_1 <cit.> and OT2_mboquien_4 <cit.> proposals. The <cit.> observations (observation ID 1342189079 and 1342189080) were performed simultaneously with the PACS (100 and 160 μm) and SPIRE (250, 350, and 500 μm) instruments in parallel mode in two orthogonal directions to map a region of approximately 90 arcmin × 90 arcmin. The <cit.> observations (observation ID 1342247408 and 1342247409) were performed solely with the PACS instrument in two orthogonal directions at 70 and 160 μm, and covered a smaller area of approximately 50 arcmin × 50 arcmin. These maps have spatial resolutions ranging from approximately 6 to 11 arcsec for the PACS maps and approximately 18 to 37 arcsec for the SPIRE maps. 

The fully reduced maps were made publicly available by the HerM33es team as Herschel User Provided Data Products, and we downloaded these data from the Herschel Science Archive.[http://archives.esac.esa.int/hsa/whsa/] Full details of the PACS and SPIRE data reduction and map-making are described in detail by <cit.> and <cit.>, respectively. Following <cit.> and <cit.>, we assume a 15 and 10 % photometric uncertainty on the extended emission in these PACS and SPIRE maps, respectively. 




 §.§ IRAS



The original IRAS measurements of M33 were presented and discussed by <cit.>. In this analysis we use the Improved Reprocessing of the IRAS Survey <cit.> data for all four IRAS bands at 12, 25, 60, and 100 μm. These data have been reprocessed resulting in an improvement in the zodiacal light subtraction, destriping, and calibration. We use the photometric uncertainties estimated by <cit.> of 5.1, 15.1, 10.4, and 13.5 % for the 12, 25, 60, and 100 μm bands, respectively.




 §.§ Spitzer



Spitzer mapped M33 as part of the Guaranteed Time Observations (PID 5, PI. R. Gehrz) and we use the MIPS 24 μm data along with the IRAC 8.0, 5.8, 4.5, and 3.6 μm data. The IRAC observations have spatial resolutions of ∼2.0 arcsec, while the MIPS 24 μm map has a resolution of 6 arcsec. A detailed description of the Spitzer photometry of M33 is provided by <cit.>. For this analysis, we downloaded the data from the Spitzer Heritage Archive[http://sha.ipac.caltech.edu/applications/Spitzer/SHA/] and reprocessed the data by subtracting the contribution from the zodiacal light, applying the extended emission correction, mosaicking the data, performing an overlap correction, and subtracting the brightest point sources. This processing was performed using mopex and apex in a similar manner to that discussed in <cit.>, and we assume a calibration uncertainty of 10% on these maps.	







§ ANALYSIS AND RESULTS









 §.§ Aperture photometry



In order to determine the integrated flux densities for M33 we applied aperture photometry to the datasets summarised in Section <ref>. For our aperture photometry analysis we used an elliptical aperture with a semi-major axis of 45 arcmin (11 kpc at the distance of M33), a semi-major to semi-minor axis ratio of 10^0.23 <cit.>, and a position angle of 22.5 degrees with respect to an equatorial reference frame <cit.>, centred on M33 (ℓ = 133.60^∘, b = -31.34^∘) as indicated in Fig. <ref>. The size of the aperture was selected based on computing the integrated flux density in apertures of increasing size to determine when the computed flux density stopped growing. The unrelated background and foreground emission was estimated within an elliptical annulus with inner and outer semi-major axes of 1.15 and 1.50 times the aperture semi-major axis, respectively.

The estimated uncertainty on the computed flux densities is a combination of the photometric uncertainty, ϵ_phot, for which we have adopted the values listed in Table <ref>, and the flux measurement uncertainty, ϵ_bg, which contains two terms: the first term is the variance in the aperture flux, and the second term is due to the background/foreground subtraction. Both terms are estimated from the variance in the annulus surrounding the source, computed as


    ϵ_bg = σ_bg[ N_aper + π N_aper^2/2 N_bg]^0.5 ,


<cit.>, where σ_bg is the standard deviation of the pixels within the annulus, and N_aper and N_bg are the number of pixels within the aperture and the annulus, respectively. The total uncertainty on the computed flux densities is then estimated to be


    ϵ = √(ϵ_phot^2 + ϵ_bg^2) .


The aperture photometry described above takes into account the average CMB contribution, but not the effect of the CMB fluctuations, which may have a considerable impact on estimates of the M33 integrated flux density. Therefore, before determining the intrinsic M33 flux density spectrum, we need to consider the effect of the CMB on the shape of the spectrum.






 §.§ CMB contribution







At the galactic latitude of M33, Milky Way foreground emission is relatively weak, even at the higher frequencies, and the main background affecting the source flux measurements is the CMB itself, especially at the lower observed frequencies. As can be seen from the left column in Fig. <ref>, M33 only starts to become clearly visible above the CMB at frequencies ≳ 143 GHz. Planck has extracted the CMB over the whole sky using four different methods: , a non-parametric method that computes the CMB map by linearly combining all of the Planck maps with weights that vary with multipole in the spherical harmonic domain; Needlet Internal Linear Combination (), which produces a CMB map using the Planck maps between 44 and 857 GHz by applying the Internal Linear Combination technique in the needlet (wavelet) domain; , which estimates a CMB map based on linear template fitting in the map domain using internal templates constructed using the Planck data; and , which is a Bayesian parametric method that models all of the astrophysical signals in the map domain <cit.>. 

It is important to note that in these CMB estimations, bright sources in the input maps are masked and the resulting CMB maps contain “inpainted” values within the masked areas. These inpainted values are good estimations of the CMB near the edge of the masked area (to preserve continuity), but are only statistically representative of the CMB within the mask. For this reason, the CMB masks employed by each of the four CMB separation techniques, along with the resulting CMB maps for the vicinity of M33, are displayed in Fig. <ref>, where it is apparent that M33 was only masked for the  analysis. In addition to providing maps of the CMB, the Planck Legacy Archive also provide maps containing only foreground emission (i.e., the frequency maps with the CMB contribution subtracted). Since there are four different CMB maps, this produces four sets of CMB-subtracted maps. As recommended by <cit.>, it is not advisable to produce an analysis that depends solely on a single component separation algorithm, and therefore we investigate the impact of the CMB by incorporating all four of the CMB-subtracted datasets in our analysis. As for the standard Planck maps discussed in Section <ref>, we used Gnomdrizz to extract 2D projected maps of the CMB-subtracted maps and subsequently converted them into MJy sr^-1. These CMB-subtracted maps are displayed in Fig. <ref>, where it is clear to see variations between the maps, reflecting the different approaches used to estimate the CMB. The differences are most obvious at the lower frequencies ≲ 200 GHz, where the CMB and its fluctuations are strongest. 

The Planck consortium recommends using the  CMB map at small angular scales, as it results in the lowest foreground residuals. However, in this case, since M33 has been masked, it is the least reliable of the four CMB-subtraction techniques. Nonetheless, since all four of the CMB-subtracted maps are in principle statistically indistinguishable, we keep all four of them in our analysis.

In order to quantify the effect of the CMB as a contaminant of the integrated emission within M33, we compared the flux densities (estimated using aperture photometry as described in Section <ref>) in each of the nine Planck bands for the five different datasets (the standard Planck frequency maps and the , , , and  CMB-subtracted maps), which we plot in Fig. <ref>. These plots clearly illustrate how the CMB impacts the estimate of the integrated flux density. Not only are there differences between the flux density computed using the standard (non CMB-subtracted) maps, but there are also variations between the different CMB estimation methods. These variations are systematic as they only depend on the map frequency and input CMB amplitude, with  and  producing the highest and lowest flux densities, respectively, and they also reflect that M33 is masked for the  analysis, but not for the other three methods. Fig. <ref> also quantitatively shows what can be seen in Fig. <ref>, which is that the CMB contamination is non-negligible at ν ≲ 217 GHz, while at higher frequencies the emission from M33 is dominant and hence the CMB contribution becomes increasing negligible.






 §.§ Contamination by CO line emission



The wide passbands of the Planck instruments cover the rest frequencies of various molecular lines. Emission from these lines will contaminate the measured broadband continuum flux densities, causing the latter to be overestimated. As discussed by <cit.>, Galactic line emission from carbon monoxide (CO) is strongly detected in the Planck bands. Specifically, only the J=1-0, J=2-1, and J=3-2 transitions of ^12CO need to be considered, as only they are strong enough to have a significant effect on the Planck HFI100, HFI217, and HFI353 bands, respectively. Although the Planck data themselves have been used to produce maps of the ^12CO emission <cit.>, these maps were produced for the velocity range of the Milky Way, which is substantially different to that of M33. We have inspected these maps and determined that they are largely unsuitable for this analysis. However, we can still investigate the magnitude of possible CO contamination by using the ^12CO maps obtained with ground-based telescopes. M33 has been observed by <cit.>, who used the Five College Radio Astronomy Observatory (FCRAO) 14 m telescope to map the J=1-0 CO emission, while <cit.> used the J=2-1 CO maps observed using the Institut de Radioastronomie Millimétrique (IRAM) 30 m telescope. Based on these data, we computed a flux density of the J=1-0 CO emission in the HFI100 band to be ∼0.1 Jy, which when compared to the total flux densities estimated in the Planck 100 GHz band accounts for ≲ 9 % of the emission. Likewise, from the IRAM data <cit.> we estimated a J=2-1 CO flux density in the HFI217 band of ∼0.7 Jy, which accounts for ≲ 4 % of the emission in that band. Finally, the compilation of CO line ratios in spiral galaxies (Israel et al., in prep.) suggests that the integrated brightness temperature ratio between the J=3-2 and J=1-0 lines is ∼0.7, predicting a corresponding J=3-2 CO flux density of ∼1.9 Jy, which is ≲ 3 % of the emission in the HFI353 band.

Therefore, we use these flux densities to correct for the contribution from CO line emission. Throughout the rest of this analysis, we include small flux density corrections by subtracting the estimated CO flux density from the observed flux densities at 100, 217, and 353 GHz in order to determine the intrinsic M33 continuum flux density spectrum as accurately as possible.

 




 §.§ Global continuum flux density spectrum and spectral energy distribution



To produce the global continuum flux density spectrum for M33, we performed aperture photometry on the Planck, Herschel, IRIS, and Spitzer maps at full angular resolution. The resulting flux densities are listed in Table <ref>. As discussed in Section <ref>, to account for the effect of the CMB, we computed the flux densities for the standard Planck maps along with the CMB-subtracted Planck maps. The flux densities listed in Table <ref> have also been corrected for CO contamination as described in Section <ref>.

For a proper analysis of the M33 continuum spectrum, we expand its frequency (wavelength) range by adding results available in the literature, in addition to those listed in Table <ref>. These include the radio flux densities from the compilation by <cit.> and from <cit.>, as well as the 2MASS J, H, and K_S band flux densities by <cit.>, the U, B, V values from <cit.>, and the GALEX far- and near-UV flux densities by <cit.>. The resulting flux density spectra, with and without CMB-subtraction, are displayed in Fig. <ref>, while the corresponding SEDs, obtained by multiplying each flux density by its frequency, are shown in Fig. <ref>.







 §.§ Decomposition of the continuum flux density spectrum



In order to quantify the resulting flux density spectra, we fitted each of them with a model simultaneously fitting contributions representing thermal dust emission, as a combination of two modified blackbodies at different temperatures (the use of two modified blackbodies is to insure an accurate fit to the peak of the cold dust component), but with identical dust emissivity indices, 


    S_ν^dust = ∑_i=1^2 C_dust,i( ν/ν_0)^β_eff B_ν(T_dust,i) ,


non-thermal synchrotron emission,


    S_ν^sync = C_sync( ν/ν_0)^α_sync ,


free-free emission, 


    S_ν^ff = C_ff( ν/ν_0)^α_ff ,


and AME, assuming that it is due to spinning dust emission,


    S_ν^AME = C_AME j_ν ,


where j_ν is the spinning dust emissivity for the warm ionised medium computed using the spinning dust model, spdust <cit.>. For each flux density spectrum we used the IDL fitting routine mpfit <cit.>, which employs the Levenberg-Marquardt least-squares minimisation technique, to fit the data between 1.4 GHz and 24 μm for C_dust,i, T_dust,i, β_eff, C_sync, α_sync, C_ff, α_ff, and C_AME. During the fitting process, C_dust,i, T_dust,i, C_sync, C_AME were constrained to be physically realistic (i.e., ≥0), β_eff and α_sync were unconstrained, while α_ff was fixed to -0.1 and C_ff, which as discussed below, was constrained based on additional observations. The estimated uncertainties on these fitted parameters were computed from the resulting covariance matrix.








In order to derive reliable thermal dust parameters, the unrelated contributions of the thermal (free-free) and non-thermal (synchrotron) emission components of the gas must be accurately determined. Unfortunately, the decomposition of the low-frequency radio continuum of galaxy flux density spectra is usually not clear-cut because of the degeneracy between the free-free and synchrotron contributions, especially since the intrinsic spectral index of the synchrotron emission is not known. This degeneracy specifically hampers the determination of any AME contribution to the observed emission spectrum and additional constraints are desirable. In the case of M33, such constraints exist. The sum of the directly measured Hii region flux densities corresponds to 235 mJy at 10 GHz <cit.>. Unfortunately, this does not include any contribution from the diffuse emission and the tally of Hii regions is incomplete, especially at large radii. It thus provides us only with a useful lower limit. However, a more accurate estimate of the total thermal radio emission may be obtained from the integrated Hα emission (I_Hα = 3.6×10^-13 W m^-2) measured by <cit.>, after first correcting for global extinction. The M33 SED shown in Fig. <ref> exhibits a peak at optical wavelengths (∼5×10^5 GHz) representing the directly observed integrated starlight, and another peak in the far-IR (∼2×10^3 GHz) representing absorbed and re-emitted starlight. Since the second peak is less than the first one, a relatively minor fraction of all starlight is intercepted by dust, and the (small) global extinction can be estimated from the ratio of the peak fluxes. The M33 SED resembles those of the SMC and, in particular, the LMC <cit.>, with an optical luminosity exceeding the IR luminosity by a factor of ∼2.5. From the luminosity ratio of the far-IR and optical peaks in Fig. <ref>, we estimate a visual extinction A_V = 0.4 ± 0.1 mag, of which 0.1 mag is due to the Milky Way foreground <cit.>. Assuming A_Hα = 0.81A_V, which is a typical Milky Way extinction curve <cit.>, and at these wavelengths is very similar to typical SMC and LMC extinction curves <cit.>, this corresponds to an Hα extinction A_Hα = 0.24 ± 0.08 mag internal to M33. Hence, the corrected Hα flux is (4.5 ± 0.5) × 10^-13 W m^-2. Using


    S^ff_ν/I_Hα = 1.15×10^-14 [1-0.21×log(ν/GHz)] Hz^-1 ,
 

we find that the free-free flux density at 10 GHz is S^ff_10 GHz = 410 ± 45 mJy. This is higher than the value of 280 mJy inferred from the work by <cit.>, but in agreement with the value of 400 mJy that follows from the determination by <cit.>. Therefore, during the fitting process we constrained C_ff to be 410 ± 45 mJy at 10 GHz, and fixed α_ff = -0.1. This is an important element in the decomposition of the observed continuum spectrum, essential for a reliable evaluation of the AME contribution in the 5–50 GHz frequency range.

The full results of our fitting analysis are shown in Fig. <ref>, including the fitted parameters (which are also listed in Table <ref>), the individual fitted components, the overall flux density spectrum fit, and the normalise residuals of the fits. It is clear that there is significant difference between the fit to the standard Planck data compared to the CMB-subtracted Planck data. Although there is a spread in the fitted T_dust and β_eff values estimated from the CMB-subtracted data (blue, pink, red, and green curves in Fig. <ref>), the fit to the standard data (black curve in Fig. <ref>) is significantly outside this range. Focusing solely on the fits to the CMB-subtracted data, and combing these four fits, we find that the dust emission spectrum of M33 between ∼100 GHz and 3 THz is adequately described by a single modified blackbody curve, with a peak of ∼2000 Jy, a mean dust temperature T_dust = 21.67 ± 0.30 K, and a mean effective dust emissivity β_eff = 1.35 ± 0.10. There is also a warm dust component with a mean temperature of 61.89 ± 0.67 K that was forced to have the same effective dust emissivity as the cold dust component. Since this warm component was simply included to insure an accurate fit to the peak of the cold dust component, we do not interpret this any further. The mean synchrotron radio continuum spectral index is α_sync = -1.03 ± 0.03. The relevant mean values are also listed in Table <ref>. For the individual entries we list the internal errors, whereas the errors given for the mean values reflect the rather larger dispersion of the individual values. Comparing the mean values to the standard values (i.e., comparing the first and last rows in Table <ref>) we find that not correcting for the CMB contribution would result in a significant over-estimate of T_dust (by ∼5 K) and under-estimate of β_eff (by ∼0.4), clearly highlighting the importance of correcting for the CMB. 

We computed the mean fraction of thermal radio emission at 20 cm and 3.6 cm to be ∼16 % and ∼49 %, respectively, which are consistent with the estimates from <cit.>, confirming that our estimate of the thermal emission from the Hα emission is reasonable. Not surprisingly, our estimated synchrotron spectral index is also consistent with the results obtained by <cit.>. Although our fitted free-free emission amplitudes match the limits of the estimate obtained from the Hα observations, we confirmed that even without imposing this constraint on the free-free emission, we find consistent results, with a mean dust temperature of T_dust = 22.36 ± 0.69 K and a mean spectral index of β_eff = 1.31 ± 0.10, but the fraction of free-free emission is decreased to ∼11 % and ∼34 % at 20 cm and 3.6 cm, respectively. Ignoring the radio data, and simply fitting the data at frequencies ≥ 100 GHz, we find a slightly lower value of β_eff = 1.28 ± 0.10. The fact that these different approaches all yield consistent results confirms that our fit is not biased by the radio data, the decomposition between the thermal and non-thermal radio emission, nor the degeneracy between the free-free emission and the AME.





We find that the AME in M33 is at best a minor component, both in an absolute sense and when compared to the free-free and synchrotron emission. Even though we modelled the AME using a spinning dust model for the warm ionised medium (as was used by <cit.> for their M31 analysis), we obtained consistent results using the cold neutral medium, warm neutral medium, or molecular cloud spinning dust models. Based on observations of AME in the Milky Way, the ratio between the AME emission at 30 GHz and the thermal dust emission at 100 μm (often incorrectly referred to as an AME emissivity) is of the order of ∼2×10^-4 <cit.>. Therefore, since we find a 100 μm flux density of ∼1350 Jy, this would lead us to expect ∼0.3 Jy of AME at 30 GHz, while we only estimate an AME flux density of ≲ 0.04 Jy. Although this ratio between the AME and thermal dust emission is sensitive to dust temperature <cit.>, these results indicate that there is significantly less AME in M33 compared to our own Galaxy, which is consistent with what has been observed in M82, NGC253, and NGC4945 <cit.>. On the other hand, in M31 the AME appears to be much more prominent, with a tentative detection that is comparable to what would be expected based on the AME level observed in our own Galaxy <cit.>. However, we note that for M31, the lack of observations between ∼1 and 20 GHz could bias the fit, which is not the case for M33, M82, NGC253, and NGC4945, where the wavelength coverage is more complete.

Finally, we emphasise that a single dust temperature and a single effective emissivity index, i.e., a single curve, provides a reasonable fit to the observed data between ∼100 GHz and 3 THz. There is no spectral break, and there is no indication of an “excess” of any kind in the global (spatially integrated) flux density spectrum of M33.










 §.§ Radial variations



In view of the radial changes in the brightness of the M33 disk, it is of interest to establish whether or not the flux density spectrum changes with radial distance from the center of M33. For this purpose, we now use the Planck, Herschel, and IRAS maps at the best common angular resolution (10 arcmin), ignoring the three Planck LFI bands, which have beams larger than 10 arcmin. Even so, the limited extent of M33 allows only three fully independent concentric ellipses to be constructed with semi-major axes of 8, 5.33, and 2.67 kpc. As before, the background/foreground emission was estimated within an elliptical annulus with inner and outer semi-major axes of 1.15 and 1.50 times the semi-major axis of the 8 kpc aperture, respectively (see Fig. <ref>).

<cit.> have shown that across the disk of M33, the radial profiles of the far-IR and radio emission are very similar. Therefore, in addition to fitting the far-IR/sub-mm emission within each annulus, we also approximated the radio emission by scaling the curves from Fig. <ref>. To do this, we computed the mean ratio between the 160 μm flux density and the 4.8 GHz flux density from the analysis in Section <ref>, and assuming that this is constant, we estimated the level of the radio emission within each annulus. We also assume that the fraction of free-free emission at 4.8 GHz is fixed across M33, and fit the total radio emission with a synchrotron spectral index identical to the mean found for the entire galaxy. As before, we fit for synchrotron, free-free, AME, and thermal dust emissions. However, for this analysis we combine the four flux density measurements at each wavelength using the scatter as a measure of the uncertainty and perform a single fit, rather than a separate fit for each of the four CMB-subtracted maps. The resulting flux density spectra for each of the three elliptical annuli representing the inner, middle, and outer regions of M33 are displayed in Fig. <ref>, along with the corresponding normalised residuals of the fits.

The radial dependence of both the dust temperature, T_dust, and the effective emissivity, β_eff, can be inferred from the plots, and it is clear that both dust temperature and the effective dust emissivity decrease with increasing radius. In the center of M33, T_dust = 22.36 ± 0.16 K and β_eff = 1.53 ± 0.01. As the radius increases, both decrease to T_dust = 21.42 ± 0.16 K and β_eff = 1.38 ± 0.01 in the middle annulus (semi-major axis between 2.67 and 5.33 kpc) and to T_dust = 19.24 ± 0.15 K and β_eff = 1.29 ± 0.02 in the outer annulus (semi-major axis between 5.33 and 8.00 kpc). We emphasize that this result cannot be caused by the T_dust, β_eff degeneracy referred to earlier as that would require β_eff to increase as T_dust decreases.

Our finding that both T_dust and β_eff decrease with radius is qualitatively similar to the result obtained by <cit.>, who used Monte Carlo simulations on more limited Herschel flux density spectra between 100 and 500 μm to deduce a simultaneous decrease of dust temperature (from ∼24 K to ∼18 K) and emissivity (from ∼1.8 to ∼1.2) going from the center of M33 out to a radius of 6 kpc. The more extended spectral coverage presented here allow quantitatively more robust results even though the spatial resolution is lower.






§ DISCUSSION



The global continuum flux density spectrum of M33 is characterised by an overall emissivity β_eff = 1.35 ± 0.10, which is below the value of 1.5 estimated from combined Herschel and Spitzer observations down to 600 GHz (500 μm) by <cit.>. The difference illustrates the bias introduced by the lack of low-frequency flux densities that most tightly constrain the Rayleigh-Jeans slope of the flux density spectrum. Even though we can fit the Rayleigh-Jeans part of the M33 global flux density spectrum with a single modified blackbody, it is a priori not likely that all of the dust in M33 radiates at a single temperature. However, a superposition of modified blackbodies representing grains with emissivities, β_g, radiating at a range of temperatures may create a profile that is observationally hard to distinguish from a single-temperature modified blackbody profile with an apparent emissivity β_eff≤β_g, especially when dust temperature and emissivity are negatively correlated as originally suggested by <cit.> and <cit.>, and later confirmed by <cit.>. Our analysis clearly shows that this is the case. The global flux density spectrum, whose Rayleigh-Jeans part is well-defined by a single modified blackbody is shown to be the sum of at least three different flux density spectra representing the inner, middle, and outer regions of M33, each with Rayleigh-Jeans sections equally well fitted by a single modified blackbody. As the number of sub-spectra is only limited by the available angular resolution, we expect that each of these in turn could be decomposed further.




 §.§ Dust mass



Using


    M_dust = S_ν d^2/κ_ν B_ν(T_dust) ,
 

where κ_ν is the dust opacity, we estimated the global dust mass of M33, along with the dust mass in each of the three annuli. It is known that values of the dust opacity in the literature can vary by orders of magnitude <cit.>, and in this work we adopt a value of κ_ν = 1.4 m^2kg^-1 at 160 μm taken from the “standard model” dust properties from <cit.>. Incorporating our results from Sections <ref> and <ref> into Equation <ref>, we estimated a global dust mass for M33 of (2.3±0.4)×10^6 M_⊙, and (0.8±0.1)×10^6 M_⊙, (1.1±0.2)×10^6 M_⊙, and (0.7±0.1)×10^6 M_⊙ for the inner, middle, and outer regions of M33, respectively. We find that our global dust mass estimated assuming a single modified blackbody is consistent, within the uncertainties, with the sum of the three dust masses estimated for the sub-regions, suggesting that fitting the entirety of M33 is degenerate with fitting the three sub-regions.




 §.§ Local Group sample



The global dust emissivity, β_eff = 1.35 ± 0.10, for M33 may also be compared to those derived from Planck observations of other Local Group galaxies: β_eff = 1.62 ± 0.10, 1.62 ± 0.11, 1.48 ± 0.25, and 1.21 ± 0.27 for the Milky Way <cit.>, M31 <cit.>, the LMC <cit.>, and the SMC <cit.>, respectively. We find that the M33 emissivity is significantly lower than that observed in the Milky Way and M31, and is more consistent with the values found in the Magellanic Clouds, with M33 actually falling between the LMC and the SMC values. Interestingly, these dust emissivities closely follow the mean metallicities of the Local Group galaxies: 12+log[O/H] = 8.32 ± 0.16, 8.67 ± 0.04, 8.72 ± 0.19, 8.43 ± 0.05, and 8.11 ± 0.03, for M33, the Milky Way, M31, the LMC, and the SMC, respectively <cit.>. To illustrate this, in Fig. <ref> we plot the dust emissivity as a function of metallicity for these 5 galaxies (filled symbols), which clearly shows that the dust emissivity increases with increasing metallicity. This trend is also observed across the M33 disk itself, where our observed emissivity gradient follows the metallicity gradient <cit.>, as can be seen when we plot the results from our three annuli within M33 (open squares) in Fig. <ref>.




 §.§ T_dust and β_eff radial variations



The apparent decrease in both T_dust and β_eff with increasing M33 radius was discussed in some detail by <cit.>. Without fully subscribing to their conclusion, we note that there are, in principle, two possible physical explanations for the observed radial decreases. The first involves dust grain composition and di-electric properties. For instance, the dust emissivity may decrease with the average interstellar energy density. Mechanical and radiative erosion of dust grains should be stronger in the more energetic inner regions than in the more quiescent outer regions. This would favour more delicate carbon/ice dust grains in the outer regions and more robust silicate-rich grains in the inner regions. The intrinsic dust grain composition may also undergo radial changes following radial gradients in the population of stellar dust producers. The second explanation involves large-scale dust cloud properties. Dust cloud heating and effective emissivity may decrease with the average radiation field, more specifically the mix of dust cloud temperatures within a specific temperature range many change as a function of irradiation. For instance, consider the possibility that each of the profiles in Fig. <ref> actually represent a collection of dust clouds and filaments with identical emissivities but different temperatures. In a radially decreasing average radiation field, clouds with temperatures at the high end would occur less frequently and have a smaller filling factor, resulting in a more skewed composite profile with a downward shift of apparent mean temperature and a consequent flattening of the Rayleigh-Jeans slope. However, the results presented in this analysis do not allow us to distinguish between these possibilities. 







 §.§ Comparison to previous studies



Our analysis is not the first to produce a full flux density spectrum of M33, as <cit.> used the Planck 2013 “nominal” mission data along with a single CMB-subtraction method () to derive the full flux density spectrum for M33. Using complex models <cit.>, they fitted their M33 flux density spectrum, deriving an excess of emission at mm/sub-mm wavelengths. However, in this analysis, we not only use the most recent Planck 2015 “full” mission data, but we also incorporate and evaluate four different CMB-subtraction techniques. As we have discussed in some depth, the contribution from the CMB fluctuations is significant and must be accurately accounted for. From Fig. <ref> and Table <ref> it is clear that there is scatter in both T_dust and β_eff between each of the four CMB-subtraction methods, highlighting the dangers of adopting a single method. 

In this analysis, we have chosen to fit our M33 flux density spectrum with a relatively simple model (i.e., a modified blackbody) and find no indication of any emission excess. A major advantage of fitting a modified blackbody to the data, rather than a complex model, is that this approach has been adopted by previous Planck analyses <cit.>, allowing direct comparisons to be made with other galaxies.






§ CONCLUSIONS



We have performed a comprehensive analysis of the global continuum flux density spectrum of M33 over a very large wavelength range from radio to UV wavelengths. In the course of this analysis, we have demonstrated the importance of accurately accounting for the contribution of CMB fluctuations to the flux density spectrum, which if neglected, results in an over-estimate of T_dust of ∼5 K and and under-estimate of β_eff of ∼0.4.

Surprisingly, we find that the global integrated emission of M33 between ∼100 GHz and 3 THz is adequately described by a single modified blackbody curve, with a mean dust temperature T_dust = 21.67 ± 0.30 K and a mean effective dust emissivity β_eff = 1.35 ± 0.10, even though such constancy of emission over all of the galaxy and throughout the line of sight is physically unlikely. In order to investigate this, we split M33 into the three independent annuli that the available resolution allows. We find that both T_dust and β_eff decrease from the centre to the outskirts of M33. This correlation is not due to any observational effect and it confirms in a direct manner an earlier conclusion reached by <cit.>. The dust emission spectrum between ∼100 GHz and 3 THz of each of the three sub-regions can be fitted with a single (but not identical) modified blackbody curve, and as the sum of the three curves (the global emission curve) itself is well-fitted with a single modified blackbody curve, we conclude that the combination of individual flux density spectra representing different parts of M33 is highly degenerate with the mean flux density spectrum.

Comparing the global far-IR emission of M33 to that of the other Local Group galaxies for which coverage is available, we find that M33 resembles the Magellanic Clouds rather than the larger spirals, the Milky Way and M31. Within this limited sample, there is a good correlation between the observed Rayleigh-Jeans slope, β_eff, and the metallicity of these galaxies, with β_eff increasing with increasing metallicity. This global correlation for the Local Group galaxies is further strengthened by the finding that it also holds within M33. The internal M33 β_eff gradient follows its metallicity gradient, with the inner part of M33 much like the LMC, and the outer part much like the lower-metallicity SMC.






§ ACKNOWLEDGMENTS

We thank the anonymous referee for providing detailed comments that have improved the content of this paper.
MWP acknowledges grant #2015/19936-1, São Paulo Research Foundation (FAPESP).







[Ali-Haïmoud et al.(2009)]Ali-Haimoud:09 
Ali-Haïmoud, Y., Hirata, C. M., & Dickinson, C. 2009, , 395, 1055 

[Bot et al.(2010)]Bot:10 
Bot, C., Ysard, N., Paradis, D., et al. 2010, , 523, A20 

[Boquien et al.(2011)]Boquien:11 
Boquien, M., Calzetti, D., Combes, F., et al. 2011, , 142, 111 

[Boquien et al.(2015)]Boquien:15 
Boquien, M., Calzetti, D., Aalto, S., et al. 2015, , 578, A8 

[Buczilowski(1988)]Buczilowski:88 
Buczilowski, U. R. 1988, , 205, 29 

[Casassus et al.(2006)]Casassus:06 
Casassus, S., Cabrera, G. F., Förster, F., et al. 2006, , 639, 951 

[Clark et al.(2016)]Clark:16 
Clark, C. J. R., Schofield, S. P., Gomez, H. L., & Davies, J. I. 2016, , 459, 1646 

[Désert et al.(2008)]Desert:08 
Désert, F.-X., Macías-Pérez, J. F., Mayet, F., et al. 2008, , 481, 411 

[de Vaucouleurs et al.(1991)]DeVaucouleurs:91 
de Vaucouleurs, G., de Vaucouleurs, A., Corwin, H. G., Jr., et al. 1991, Third Reference Catalogue of Bright Galaxies. Volume I: Explanations and references.  Volume II: Data for galaxies between 0^h and 12^h.  Volume III: Data for galaxies between 12^h and 24^h., by de Vaucouleurs, G.; de Vaucouleurs, A.; Corwin, H. G., Jr.; Buta, R. J.; Paturel, G.; Fouqué, P.. Springer, New York, NY (USA), 1991, 2091 p., ISBN 0-387-97552-7, Price US$ 198.00. ISBN 3-540-97552-7, Price DM 448.00. ISBN 0-387-97549-7 (Vol. I), ISBN 0-387-97550-0 (Vol. II), ISBN 0-387-97551-9 (Vol. III)., I,  

[Draine & Lazarian(1998)]Draine:98 
Draine, B. T., & Lazarian, A. 1998, , 508, 157 

[Draine & Hensley(2012)]Draine:12 
Draine, B. T., & Hensley, B. 2012, , 757, 103 

[Druard et al.(2014)]Druard:14 
Druard, C., Braine, J., Schuster, K. F., et al. 2014, , 567, A118 

[Dupac et al.(2003)]Dupac:03 
Dupac, X., Bernard, J.-P., Boudet, N., et al. 2003, , 404, L11 

[Fitzpatrick & Massa(2007)]Fitzpatrick:07 
Fitzpatrick, E. L., & Massa, D. 2007, , 663, 320 

[Freedman et al.(1991)]Freedman:91 
Freedman, W. L., Wilson, C. D., & Madore, B. F. 1991, , 372, 455 

[Galliano et al.(2005)]Galliano:05 
Galliano, F., Madden, S. C., Jones, A. P., Wilson, C. D., & Bernard, J.-P. 2005, , 434, 867 

[Galliano et al.(2011)]Galliano:11 
Galliano, F., Hony, S., Bernard, J.-P., et al. 2011, , 536, A88 

[Gordon et al.(2003)]Gordon:03 
Gordon, K. D., Clayton, G. C., Misselt, K. A., Landolt, A. U., & Wolff, M. J. 2003, , 594, 279 

[Górski et al.(2005)]Gorski:05 
Górski, K. M., Hivon, E., Banday, A. J., et al. 2005, , 622, 759 

[Groves et al.(2008)]Groves:08 
Groves, B., Dopita, M. A., Sutherland, R. S., et al. 2008, , 176, 438-456 

[Hermelo et al.(2016)]Hermelo:16 
Hermelo, I., Relaño, M., Lisenfeld, U., et al. 2016, , 590, A56 

[Heyer et al.(2004)]Heyer:04 
Heyer, M. H., Corbelli, E., Schneider, S. E., & Young, J. S. 2004, , 602, 723 

[Hoopes & Walterbos(2000)]Hoopes:00 
Hoopes, C. G., & Walterbos, R. A. M. 2000, , 541, 597 

[Israel(1980)]Israel:80 
Israel, F. P. 1980, , 90, 246 

[Israel et al.(1992)]Israel:92 
Israel, F. P., Mahoney, M. J., & Howarth, N. 1992, , 261, 47 

[Israel & van der Kruit(1974)]Israel:74 
Israel, F. P., & van der Kruit, P. C. 1974, , 32, 363 

[Israel et al.(2010)]Israel:10 
Israel, F. P., Wall, W. F., Raban, D., et al. 2010, , 519, A67 

[Jarrett et al.(2003)]Jarrett:03 
Jarrett, T. H., Chester, T., Cutri, R., Schneider, S. E., & Huchra, J. P. 2003, , 125, 525 

[Juvela & Ysard(2012a)]Juvela:12a 
Juvela, M., & Ysard, N. 2012a, , 539, A71 

[Juvela & Ysard(2012b)]Juvela:12b 
Juvela, M., & Ysard, N. 2012b, , 541, A33 

[Kramer et al.(2010)]Kramer:10 
Kramer, C., Buchbender, C., Xilouris, E. M., et al. 2010, , 518, L67 

[Laher et al.(2012)]Laher:12 
Laher, R. R., Gorjian, V., Rebull, L. M., et al. 2012, , 124, 737 

[Lee et al.(2011)]Lee:11 
Lee, J. C., Gil de Paz, A., Kennicutt, R. C., Jr., et al. 2011, , 192, 6 

[Markwardt(2009)]Markwardt:09 
Markwardt, C. B. 2009, Astronomical Data Analysis Software and Systems XVIII, 411, 251 

[Miville-Deschênes & Lagache(2005)]Miville-Deschenes:05 
Miville-Deschênes, M.-A., & Lagache, G. 2005, , 157, 302 

[Pagel(2003)]Pagel:03 
Pagel, B. E. J. 2003, CNO in the Universe, 304, 187 

[Paradis et al.(2012)]Paradis:12 
Paradis, D., Dobashi, K., Shimoikura, T., et al. 2012, , 543, A103 

[Paturel et al.(2003)]Paturel:03 
Paturel, G., Petit, C., Prugniel, P., et al. 2003, , 412, 45 

[Peel et al.(2011)]Peel:11 
Peel, M. W., Dickinson, C., Davies, R. D., Clements, D. L., & Beswick, R. J. 2011, , 416, L99 

[Planck Collaboration I.(2011)]Planck_Early_Results_I:11 
Planck Collaboration I. 2011, , 536, A1 

[Planck Collaboration XVI.(2011)]Planck_Early_Results_XVI:11 
Planck Collaboration XVI. 2011, , 536, A16

[Planck Collaboration XVII.(2011)]Planck_Early_Results_XVII:11 
Planck Collaboration XVII. 2011, , 536, A17 

[Planck Collaboration XXV.(2015)]Planck_Intermediate_Results_XXV:15 
Planck Collaboration XXV. 2015, , 582, A28 

[Planck Collaboration XV.(2014)]Planck_Intermediate_Results_XV:14 
Planck Collaboration XV. 2014, , 565, A103 

[Planck Collaboration IX.(2014)]Planck_2013_Results_IX:14 
Planck Collaboration IX. 2014, , 571, A9 

[Planck Collaboration XI.(2014)]Planck_2013_Results_XI:14 
Planck Collaboration XI. 2014, , 571, A11 

[Planck Collaboration XIII.(2014)]Planck_2013_Results_XIII:14 
Planck Collaboration XIII. 2014, , 571, A13 

[Planck Collaboration XXVIII.(2014)]Planck_2013_Results_XXVIII:14 
Planck Collaboration XXVIII. 2014, , 571, A28 

[Planck Collaboration I.(2016)]Planck_2015_Results_I:16 
Planck Collaboration I. 2016, , 594, A1
 
[Planck Collaboration II.(2016)]Planck_2015_Results_II:16 
Planck Collaboration II. 2016, , 594, A2 

[Planck Collaboration III.(2016)]Planck_2015_Results_III:16 
Planck Collaboration III. 2016, , 594, A3

[Planck Collaboration IV.(2016)]Planck_2015_Results_IV:16 
Planck Collaboration IV. 2016, , 594, A4

[Planck Collaboration V.(2016)]Planck_2015_Results_V:16 
Planck Collaboration V. 2016, , 594, A5

[Planck Collaboration VI.(2016)]Planck_2015_Results_VI:16 
Planck Collaboration VI. 2016, , 594, A6

[Planck Collaboration VII.(2016)]Planck_2015_Results_VII:16 
Planck Collaboration VII. 2016, , 594, A7

[Planck Collaboration VIII.(2016)]Planck_2015_Results_VIII:16 
Planck Collaboration VIII. 2016, , 594, A8

[Planck Collaboration IX.(2016)]Planck_2015_Results_IX:16 
Planck Collaboration IX. 2016, , 594, A9

[Planck Collaboration X.(2016)]Planck_2015_Results_X:16 
Planck Collaboration X. 2016, , 594, A10

[Popescu et al.(2011)]Popescu:11 
Popescu, C. C., Tuffs, R. J., Dopita, M. A., et al. 2011, , 527, A109 

[Rice et al.(1990)]Rice:90 
Rice, W., Boulanger, F., Viallefond, F., Soifer, B. T., & Freedman, W. L. 1990, , 358, 418 

[Shetty et al.(2009)]Shetty:09 
Shetty, R., Kauffmann, J., Schnee, S., & Goodman, A. A. 2009, , 696, 676 

[Silsbee et al.(2011)]Silsbee:11 
Silsbee, K., Ali-Haïmoud, Y., & Hirata, C. M. 2011, , 411, 2750 

[Tabatabaei et al.(2007a)]Tabatabaei:07a 
Tabatabaei, F. S., Krause, M., & Beck, R. 2007a, , 472, 785 

[Tabatabaei et al.(2007b)]Tabatabaei:07b 
Tabatabaei, F. S., Beck, R., Krügel, E., et al. 2007b, , 475, 133 

[Tabatabaei et al.(2014)]Tabatabaei:14 
Tabatabaei, F. S., Braine, J., Xilouris, E. M., et al. 2014, , 561, A95 

[Tauber et al.(2010)]Tauber:10 
Tauber, J. A., Mandolesi, N., Puget, J.-L., et al. 2010, , 520, A1 

[Tibbs et al.(2010)]Tibbs:10 
Tibbs, C. T., Watson, R. A., Dickinson, C., et al. 2010, , 402, 1969 

[Tibbs et al.(2011)]Tibbs:11 
Tibbs, C. T., Flagey, N., Paladini, R., et al. 2011, , 418, 1889 

[Tibbs et al.(2012)]Tibbs:12b 
Tibbs, C. T., Paladini, R., & Dickinson, C. 2012, Advances in Astronomy, 2012, 124931 

[Tibbs et al.(2013)]Tibbs:13 
Tibbs, C. T., Scaife, A. M. M., Dickinson, C., et al. 2013, , 768, 98 
	
[Todorović et al.(2010)]Todorovic:10 
Todorović, M., Davies, R. D., Dickinson, C., et al. 2010, , 406, 1629 	
	
[Toribio San Cipriano et al.(2016)]Toribio:16 
Toribio San Cipriano, L., García-Rojas, J., Esteban, C., Bresolin, F., & Peimbert, M. 2016, , 458, 1866 

[Verley et al.(2007)]Verley:07 
Verley, S., Hunt, L. K., Corbelli, E., & Giovanardi, C. 2007, , 476, 1161 

[Verley et al.(2009)]Verley:09 
Verley, S., Corbelli, E., Giovanardi, C., & Hunt, L. K. 2009, , 493, 453 

[Xilouris et al.(2012)]Xilouris:12 
Xilouris, E. M., Tabatabaei, F. S., Boquien, M., et al. 2012, , 543, A74 













