# A multiwavelength study of the Magellanic-type galaxy NGC 4449-I. Modelling the spectral energy distribution, the ionization structure and the star formation history

O. L. Karczewski, M. J. Barlow, M. J. Page, N. P. M. Kuin, I. Ferreras, M.

Baes, G. J. Bendo, A. Boselli, A. Cooray, D. Cormier, et al.

<figure><img src="image_1.png" /></figure>

► To cite this version:

O. L. Karczewski, M. J. Barlow, M. J. Page, N. P. M. Kuin, I. Ferreras, et al. A multiwavelength study of the Magellanic-type galaxy NGC 4449-I. Modelling the spectral energy distribution, the ionization structure and the star formation history. Monthly Notices of the Royal Astronomical Society, 2013, 431 (3), pp.2493-2512. 10.1093/mnras/stt345. hal-01439482

HAL Id: hal-01439482

https://hal.science/hal-01439482

Submitted on 20 Sep 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>
<figure><img src="image_3.png" /></figure>

# A multiwavelength study of the Magellanic-type galaxy NGC 4449 - I. Modelling the spectral energy distribution, the ionization structure and the star formation history

O. Ł. Karczewski,<fnref n="1" />* M. J. Barlow,<fnref n="1" /> M. J. Page,<fnref n="2" /> N. P. M. Kuin,<fnref n="2" /> I. Ferreras,<fnref n="2" /> M. Baes,<fnref n="3" /> G. J. Bendo,<fnref n="4" /> A. Boselli,<fnref n="5" /> A. Cooray,<fnref n="6" /> D. Cormier,<fnref n="7" /> I. De Looze,<fnref n="3" /> M. Galametz,<fnref n="8" /> F. Galliano,<fnref n="7" /> V. Lebouteiller,<fnref n="7" /> S. C. Madden,<fnref n="7" /> M. Pohlen,<fnref n="9" /> A. Rémy-Ruyer,<fnref n="7" /> M. W. L. Smith<fnref n="9" /> and L. Spinoglio<fnref n="10" />

Accepted 2013 February 21. Received 2013 February 20; in original form 2013 January 8

## ABSTRACT

We present an integrated photometric spectral energy distribution (SED) of the Magellanic-type galaxy NGC 4449 from the far-ultraviolet (UV) to the submillimetre, including new observations acquired by the Herschel Space Observatory. We include integrated UV photometry from the Swift Ultraviolet and Optical Telescope using a measurement technique which is appropriate for extended sources with coincidence loss. In this paper, we examine the available multiwavelength data to infer a range of ages, metallicities and star formation rates for the underlying stellar populations, as well as the composition and the total mass of dust in NGC 4449. Our analysis of the global optical spectrum of NGC 4449 fitted using the spectral fitting code STARLIGHT suggests that the majority of stellar mass resides in old ($\gtrsim$1 Gyr old) and metal-poor ($Z/Z_{\odot} \sim 0.2$) populations, with the first onset of star formation activity deduced to have taken place at an early epoch, approximately 12 Gyr ago. A simple chemical evolution model, suitable for a galaxy continuously forming stars, suggests a ratio of carbon to silicate dust mass comparable to that of the Large Magellanic Cloud over the inferred time-scales. We present an iterative scheme, which allows us to build an in-depth and multicomponent representation of NGC 4449 'bottom-up', taking advantage of the broad capabilities of the photoionization and radiative transfer code MOCASSIN (MOnte CArlo SimulationS of Ionized Nebulae). We fit the observed SED, the global ionization structure and the emission line intensities, and infer a recent star formation rate of $0.4 \ M_{\odot} \ yr^{-1}$ and a total stellar mass of $\approx 1 \times 10^9 \ M_{\odot}$ emitting with a bolometric luminosity of $5.7 \times 10^9 \ L_{\odot}$. Our fits yield a total dust mass of $2.9 \pm 0.5 \times 10^6 \ M_{\odot}$ including 2 per cent attributed to polycyclic aromatic hydrocarbons. We deduce a dust to gas mass ratio of 1/190 within the modelled region. While we do not consider possible additional contributions from even colder dust, we

<footnote n="1">¹Department of Physics and Astronomy, University College London, Gower Street, London WC1E 6BT, UK</footnote>

<footnote n="2">²Mullard Space Science Laboratory, University College London, Holmbury St. Mary, Dorking RH5 6NT, UK</footnote>

<footnote n="3">³Sterrenkundig Observatorium, Universiteit Gent, Krijgslaan 281-S9, 9000 Gent, Belgium</footnote>

<footnote n="4">⁴UK ALMA Regional Centre Node, Jodrell Bank Centre for Astrophysics, University of Manchester, Oxford Road, Manchester M13 9PL, UK</footnote>

<footnote n="5">⁵LAM, Université d'Aix-Marseille & CNRS, UMR7326, 38 rue F. Joliot-Curie, 13388 Marseille, France</footnote>

<footnote n="6">⁶Department of Physics & Astronomy, University of California, Irvine, CA 92697, USA</footnote>

<footnote n="7">⁷Laboratoire AIM, CEA/DSM-CNRS-Université Paris Diderot, IRFU/Service d'Astrophysique, CEA Saclay, L'Orme des Merisiers, Bât. 709, 91191 Gif-sur-Yvette, France</footnote>

<footnote n="8">⁸Institute of Astronomy, University of Cambridge, Madingley Road, Cambridge CB3 0HA, UK</footnote>

<footnote n="9">⁹School of Physics and Astronomy, Cardiff University, The Parade, Cardiff CF24 3AA, UK</footnote>

<footnote n="10">¹⁰Istituto di Astrofisica e Planetologia Spaziali, INAF-IAPS, Via Fosso del Cavaliere 100, I-00133 Roma, Italy</footnote>

<footnote>*E-mail: olk@star.ucl.ac.uk</footnote>
note that including the extended H I envelope and the molecular gas is likely to bring the ratio down to as low as ~1/800.

**Key words:** methods: numerical – dust, extinction – galaxies: dwarf – galaxies: individual: NGC 4449 – galaxies: ISM – galaxies: stellar content.

# 1 INTRODUCTION

The conditions in the interstellar medium (ISM) in a galaxy and the environment-dependent processes governing the formation of stars and dust can be studied on many physical scales. Considering an entire galaxy, the interplay between the underlying stellar emission, the degree of radiation reprocessing taking place in the ISM and the thermal or non-thermal emission due to dynamical or evolutionary processes give rise to a unique pattern of emission and absorption features in the observed spectral energy distribution (SED).

The spectroscopic measurements alone provide a wealth of information about the ionized (H II) regions or photodissociation regions (PDRs) in galaxies (e.g., Guseva et al. 2004; Vasta et al. 2010). Studies constrained by broad-band photometric measurements can allow a better understanding of the main components of the thermal continuum, namely, the stellar emission and the emission due to interstellar dust. Observed emission lines and diagnostic line ratios provide important further constraints (Martínez-Galarza et al. 2011). Synthetic single stellar populations (SSPs) can be used to characterize star-forming regions within galaxies (e.g. Groves et al. 2008; Ferreras et al. 2012), or to decompose the observed SEDs of entire galaxies (e.g. Amorín et al. 2012). The composition of dust and the masses of dust species can be inferred from observations by assuming measured laboratory properties of dust (e.g. Galliano, Dwek & Chanial 2008a and references therein). Together, these components can provide an in-depth view of a galaxy (e.g. Cormier et al. 2012; Hermelo et al. 2013).

Numerical models reproducing all the available photometric and spectroscopic measurements across a wide range of wavelengths can provide the most complete picture of the integrated properties of a galaxy. We wish to examine a multicomponent model treating the stellar content, the gaseous phase and the dust within a galaxy in a self-consistent way. Our initial goals include finding a robust numerical scheme, which could be applied to a sample of dwarf galaxies in order to study the details of their star formation histories (SFHs) and dust content.

In this paper, we describe the first multicomponent observation-driven model of a galaxy generated with the photoionization and radiative transfer code MOnte CArlo SimulationS of Ionized Nebulae (MOCASSIN; Ercolano et al. 2003, 2008; Ercolano, Barlow & Storey 2005), which includes a simultaneous and self-consistent treatment of the stellar component, the gaseous phase, its ionization structure and multiple dust species. MOCASSIN is a fully three-dimensional code, whose recent applications include modelling the extreme bipolar planetary nebula NGC 6302 (Wright et al. 2011) and studying dust emission by supernova (SN) ejecta (Wesson et al. 2010). Applying MOCASSIN to simulate entire galaxies is more resource consuming since both the gas phase and the dust phase need to be treated self-consistently over relatively large physical scales. The small intrinsic size of dwarf galaxies make them ideal candidates for such multiwavelength studies. Indeed, dwarf galaxies provide small self-contained but dynamic environments. As such, they are useful in developing our understanding of mechanisms underlying star formation, including chemical enrichment and feedback

processes, especially for low-metallicity objects resembling those formed at earlier epochs in the evolution of the Universe.

In a similar context, MOCASSIN was first used to model Mrk 996 (James 2009). In this paper, we describe a model of the well-known starburst galaxy NGC 4449, which is a low-metallicity [log(O/H) + 12 = 8.23, $1/3 Z_\odot$; Engelbracht et al. 2008] actively star forming (SFR ~ 0.5 $M_\odot$ yr$^{-1}$; Hill et al. 1998; Hunter et al. 1998) barred Magellanic-type irregular galaxy (de Vaucouleurs et al. 1991; Corwin, Buta & de Vaucouleurs 1994) seen face-on. At a distance of 3.8 Mpc (Annibali et al. 2008) its H I envelope extends to 12.9 kpc (11.6 arcmin), equal to approximately three Holmberg radii (Hunter, van Woerden & Gallagher 1999; Swaters et al. 2002), and forms a pronounced system of streamers with a counter-rotating H I core, which may be indicative of past mergers (Hunter et al. 1999; Theis & Kohle 2001). NGC 4449 is particularly suitable for multiwavelength studies because of its wide range of existing photometric data extending to the far-infrared (FIR). The new observations acquired with the Spitzer Space Telescope (Werner et al. 2004) and the Herschel Space Observatory (Pilbratt et al. 2010) are essential in studying gas cooling processes and the properties of dust in the ISM of NGC 4449.

In the following sections, we will discuss in detail a model of NGC 4449 based on multiple observational constraints and including a treatment of polycyclic aromatic hydrocarbons (PAHs). Thus, in Section 2 we describe the observations and data reduction techniques, and in Section 3 we describe the modelling method, as well as the assumptions and the convergence criteria adopted. The results and discussion are presented in Section 4, and the conclusions follow in Section 5. In a companion paper, Karczewski et al. (in preparation, hereafter Paper II), we study the FIR cooling lines and discuss properties of PDRs inferred from spatially resolved spectroscopic observations of NGC 4449.

# 2 PHOTOMETRIC OBSERVATIONS AND DATA REDUCTION

In constructing the observed SED of NGC 4449 we have used data from previous studies (e.g. Hunter et al. 1986; Kennicutt 1992; Böttner, Klein & Heithausen 2003; Engelbracht et al. 2008; Bendo, Galliano & Madden 2012b) and from systematic surveys [Sloan Digital Sky Survey (SDSS); York et al. 2000, Two Micron All Sky Survey (2MASS); Skrutskie et al. 2006, *WISE*; Wright et al. 2010 and *Planck*; Planck Collaboration 2011]. We present archival data from *GALEX* [Martin et al. 2005, ultraviolet (UV) wavebands] and *Swift* (Gehrels et al. 2004, UV and optical wavebands), as well as new FIR photometry from *Herschel* (Pilbratt et al. 2010). A journal of observations is given in Table 1.

In Table 2, we list all global photometric measurements used to construct the observed SED for NGC 4449. Existing submillimetre SCUBA (Submillimetre Common-User Bolometer Array) measurements at 450 and 850 µm (Böttner et al. 2003) were omitted from the list of photometric measurements due to their incomplete field of view (FOV) and the availability of space borne *Herschel* and *Planck* observations offering a similar wavelength coverage.
<table><caption>**Table 1.** List of observations.</caption><thead><tr><th>Instrument</th><th>Bandpass/λ<sub>eff</sub></th><th>Time (s)</th><th>Date observed</th><th>Observation ID</th></tr></thead><tbody><tr><td colspan="5" style="text-align: center;">New observations</td></tr><tr><td rowspan="4">Herschel/PACS</td><td>70 + 160 µm</td><td>3497</td><td>2011 May 16</td><td>1342221125</td></tr><tr><td>70 + 160 µm</td><td>3497</td><td>2011 May 16</td><td>1342221126</td></tr><tr><td>100 + 160 µm</td><td>3497</td><td>2011 May 16</td><td>1342221127</td></tr><tr><td>100 + 160 µm</td><td>3497</td><td>2011 May 16</td><td>1342221128</td></tr><tr><td>Herschel/SPIRE</td><td>250, 350, 500 µm</td><td>1035</td><td>2010 June 12</td><td>1342198243</td></tr><tr><td colspan="5" style="text-align: center;">Archival observations</td></tr><tr><td rowspan="7">Swift/UVOT</td><td>uvw2/2030 Å</td><td>1544</td><td>2007 Mar 27</td><td>00035873001</td></tr><tr><td>uvm2/2231 Å</td><td>1091</td><td>2007 Mar 27</td><td>00035873001</td></tr><tr><td>uvw1/2634 Å</td><td>771</td><td>2007 Mar 27</td><td>00035873001</td></tr><tr><td>u/3501 Å</td><td>385</td><td>2007 Mar 27</td><td>00035873001</td></tr><tr><td>b/4329 Å</td><td>385</td><td>2007 Mar 27</td><td>00035873001</td></tr><tr><td>v/5402 Å</td><td>385</td><td>2007 Mar 27</td><td>00035873001</td></tr><tr><td>FUV/1539 Å</td><td>835</td><td>2006 Mar 15</td><td>tile 5228</td></tr><tr><td rowspan="2">GALEX</td><td>NUV/2316 Å</td><td>765</td><td>2006 Mar 15</td><td>tile 5228</td></tr><tr><td>u/3551 Å</td><td>54</td><td>2003 Mar 24</td><td>3813/1/41/237/241</td></tr><tr><td rowspan="5">SDSS</td><td>g/4686 Å</td><td>54</td><td>2003 Mar 24</td><td>3813/1/41/237/245</td></tr><tr><td>r/6165 Å</td><td>54</td><td>2003 Mar 24</td><td>3813/1/41/237/237</td></tr><tr><td>i/7481 Å</td><td>54</td><td>2003 Mar 24</td><td>3813/1/41/237/239</td></tr><tr><td>z/8931 Å</td><td>54</td><td>2003 Mar 24</td><td>3813/1/41/237/243</td></tr></tbody></table>

Five illustrative false-colour images of NGC 4449 over selected broad-band ranges are presented in Fig. 1.

## 2.1 Ultraviolet observations with GALEX and Swift

At UV wavelengths NGC 4449 has been observed by two recent space-borne missions, GALEX (in two UV bands) and the Swift Ultraviolet/Optical Telescope (UVOT; in three UV and three optical bands). Standard aperture photometry using the in-flight calibration of Morrissey et al. (2007) was performed on tile 5228 acquired as part of GALEX Nearby Galaxy Survey and accessible via General Release 6. The overall uncertainties include photometric repeatability measurements by Morrissey et al. (2007).

Swift/UVOT (Roming et al. 2005) offers four times the angular resolution of GALEX in three narrow UV bands. The telescope is equipped with a photon-counting detector and was originally designed to detect and observe gamma-ray bursts. As a consequence, standard data reduction procedures cannot be applied straightforwardly to observations of extended sources. In the Appendix we discuss this problem in detail and present a method of obtaining global photometry, which is appropriate for extended sources.

## 2.2 Optical observations with SDSS

In addition to the optical fluxes derived from the Swift/UVOT observations, SDSS observations in five optical bands (u, g, r, i and z) are available as part of Data Release 7. NGC 4449 is contained entirely within one field, thus, minimizing deblending or sky subtraction issues which may arise in the case of extended sources. Standard aperture photometry was performed using elliptical apertures and the counts were converted into flux densities as described by West et al. (2010). The uncertainties in the derived fluxes are dominated by uncertainties in sky determination and subtraction.

## 2.3 Far-infrared observations with Herschel

NGC 4449 has been observed by Herschel (Pilbratt et al. 2010) as part of the Dwarf Galaxy Survey (Madden et al. 2013) with its two imaging photometers: the Photodetector Array Camera and Spectrometer (PACS; Poglitsch et al. 2010) at 70, 100 and 160 µm, and the Spectral and Photometric Imaging REceiver (SPIRE; Griffin et al. 2010) at 250, 350 and 500 µm. The full width at half-maximum (FWHM) of the point spread function (PSF) is 5.2, 7.7, 12.0, 18.2, 24.9, 36.3 arcsec in these bands, respectively. The details of all PACS and SPIRE data reduction steps, including error estimation, can be found in Rémy-Ruyer et al. (in preparation).

The PACS observations were performed as four pairs of orthogonal scans covering an area of 24 arcmin × 24 arcmin. We used an adapted version of the standard script of v7.0 of the Herschel Interactive Processing Environment (HIPE; Ott 2010). The basic processing includes flagging bad or saturated pixels, converting the signal into Jy pixel⁻¹ and applying the flat-field correction. Additionally, we systematically masked column 0 of all of the constituent 16 × 16 matrices in the PACS array to avoid electronic crosstalk and we performed second level deglitching. The resulting Level 1 products were converted into maps with the pixel size of 2, 2 and 4 arcsec for the three PACS bands and processed with Scanamorphos (Roussel 2012), which is particularly suitable for extended sources with low frequency noise. The integrated PACS fluxes of NGC 4449 are well represented by a 40 K blackbody and suitable colour corrections were applied accordingly (Müller, Okumura & Klaas 2011b), where no additional correction factors are required to correct for the extended nature of the source (Sauvage 2011). The quoted uncertainties include the calibration uncertainties at 5 per cent in all bands (Müller et al. 2011a; Rémy-Ruyer et al., in preparation).

The SPIRE observations were performed as two orthogonal scans covering an area of 24 arcmin × 24 arcmin. The corresponding maps were reduced using a modified version of the SPIRE pipeline in HIPE. The steps up to Level 1 were identical to those in the
<table><caption>**Table 2.** Summary of the derived photometric data for NGC 4449. All fluxes have been corrected for the foreground extinction using $E(B - V) = 0.019$ (Schlegel, Finkbeiner & Davis 1998) and the extinction law of Cardelli, Clayton & Mathis (1989) with $R_V = 3.1$. References: (1) this work; (2) 2MASS All-Sky Extended Source Catalog (Skrutskie et al. 2006); (3) WISE All-Sky Source Catalog (Wright et al. 2010); (4) Engelbracht et al. (2008); (5) Hunter et al. (1986); (6) Bendo et al. (2012b); (7) Rémy-Ruyer et al. (in preparation); (8) Planck Early Release Compact Source Catalogue (Planck Collaboration 2011) and (9) Böttner et al. (2003).</caption><thead><tr><th>Survey/instrument</th><th>λ<sub>eff</sub></th><th>F<sub>v</sub></th><th>Aperture<sup>a</sup></th><th>Reference</th></tr></thead><tbody><tr><td rowspan="2">GALEX</td><td>1539 Å</td><td>152 ± 12 mJy</td><td>11.5 arcmin</td><td rowspan="2">(1)</td></tr><tr><td>2316 Å</td><td>183 ± 9 mJy</td><td>7.5 arcmin</td></tr><tr><td rowspan="6">Swift/UVOT</td><td>1991 Å</td><td>175 ± 12 mJy</td><td rowspan="6">8.5 arcmin × 6 arcmin</td><td rowspan="6">(1)</td></tr><tr><td>2221 Å</td><td>189 ± 14 mJy</td></tr><tr><td>2486 Å</td><td>202 ± 15 mJy</td></tr><tr><td>3442 Å</td><td>249 ± 22 mJy</td></tr><tr><td>4321 Å</td><td>463 ± 41 mJy</td></tr><tr><td>5410 Å</td><td>611 ± 57 mJy</td></tr><tr><td rowspan="5">SDSS</td><td>3551 Å</td><td>242 ± 5 mJy</td><td rowspan="5">6.5 arcmin × 4.7 arcmin</td><td rowspan="5">(1)</td></tr><tr><td>4686 Å</td><td>457 ± 14 mJy</td></tr><tr><td>6165 Å</td><td>573 ± 40 mJy</td></tr><tr><td>7481 Å</td><td>622 ± 50 mJy</td></tr><tr><td>8931 Å</td><td>683 ± 34 mJy</td></tr><tr><td rowspan="3">2MASS</td><td>1.24 µm</td><td>916 ± 21 mJy</td><td rowspan="3">8 arcmin<sup>b</sup></td><td rowspan="3">(2)</td></tr><tr><td>1.66 µm</td><td>1070 ± 30 mJy</td></tr><tr><td>2.16 µm</td><td>839 ± 31 mJy</td></tr><tr><td rowspan="4">WISE</td><td>3.4 µm</td><td>391 ± 7 mJy</td><td rowspan="4">4.0 arcmin × 2.8 arcmin</td><td rowspan="4">(3)</td></tr><tr><td>4.6 µm</td><td>251 ± 4 mJy</td></tr><tr><td>12 µm</td><td>954 ± 17 mJy</td></tr><tr><td>22 µm</td><td>2860 ± 60 mJy</td></tr><tr><td rowspan="4">Spitzer/IRAC</td><td>3.6 µm</td><td>493 ± 15 mJy</td><td rowspan="4">–</td><td rowspan="4">(4)</td></tr><tr><td>4.5 µm</td><td>317 ± 10 mJy</td></tr><tr><td>5.8 µm</td><td>615 ± 19 mJy</td></tr><tr><td>8 µm</td><td>1420 ± 40 mJy</td></tr><tr><td rowspan="3">IRAS</td><td>12 µm</td><td>2.1 ± 0.2 Jy</td><td rowspan="3">8 arcmin</td><td rowspan="3">(5)</td></tr><tr><td>25 µm</td><td>4.7 ± 0.5 Jy</td></tr><tr><td>60 µm</td><td>36 ± 3 Jy</td></tr><tr><td rowspan="3">Spitzer/MIPS</td><td>24 µm</td><td>3.29 ± 0.13 Jy</td><td rowspan="3">9.3 arcmin × 6.6 arcmin</td><td rowspan="3">(6)</td></tr><tr><td>70 µm</td><td>43.8 ± 4.4 Jy</td></tr><tr><td>160 µm</td><td>78.1 ± 9.4 Jy</td></tr><tr><td rowspan="3">Herschel/PACS</td><td>70 µm</td><td>49.3 ± 2.5 Jy</td><td rowspan="3">8.3 arcmin</td><td rowspan="3">(7)</td></tr><tr><td>100 µm</td><td>75.9 ± 3.8 Jy</td></tr><tr><td>160 µm</td><td>79.5 ± 4.0 Jy</td></tr><tr><td rowspan="3">Herschel/SPIRE</td><td>250 µm</td><td>35.7 ± 2.5 Jy</td><td rowspan="3">8.3 arcmin</td><td rowspan="3">(7)</td></tr><tr><td>350 µm</td><td>16.2 ± 1.1 Jy</td></tr><tr><td>500 µm</td><td>5.6 ± 0.4 Jy</td></tr><tr><td rowspan="3">Planck</td><td>350 µm</td><td>16.14 ± 0.34 Jy</td><td>6.1 arcmin × 5.2 arcmin<sup>c</sup></td><td rowspan="3">(8)</td></tr><tr><td>550 µm</td><td>4.97 ± 0.24 Jy</td><td>5.8 arcmin × 4.8 arcmin<sup>c</sup></td></tr><tr><td>850 µm</td><td>1.45 ± 0.15 Jy</td><td>6.0 arcmin × 4.4 arcmin<sup>c</sup></td></tr><tr><td>IRAM</td><td>1.2 mm</td><td>260 ± 40 mJy</td><td>Two fields 2.3 arcmin</td><td>(9)</td></tr></tbody></table>

original pipeline provided by the SPIRE Instrument Control Center. Additionally, residual baseline subtraction was performed by subtracting the median of the time lines for each bolometer over the entire observation (Pohlen et al. 2010). This was followed by an iterative process to completely remove residual signals that appear as stripes in the maps (Bendo et al. 2010). The final map was constructed using Naive Mapper available in HIPE and calibrated for an extended source. A modified blackbody fit of the form $S_v \propto v^\beta B_v(T)$ yields $\beta \approx 2$ globally for NGC 4449 (Rémy-Ruyer et al., in preparation), which was used for colour correction of the integrated fluxes (SPIRE Observers' Manual; Valtchanov 2011). The uncertainties

include the revised overall calibration uncertainties at 7 per cent in all SPIRE bands (Griffin & Lim 2011).

## 2.4 2MASS, WISE and Planck catalogue data

Global 2MASS photometric measurements for NGC 4449 were taken from columns `j_m_ext`, `h_m_ext` and `k_m_ext` of the 2MASS All-Sky Extended Source Catalog, and were subsequently converted to flux densities using the zero-points tabulated by Cohen, Wheaton & Megeath (2003). WISE photometry was taken from columns `gmag` and `gerr` of the WISE All-Sky Source

<footnote><sup>a</sup>Given as diameters for circular apertures or major × minor axes for elliptical apertures. <sup>b</sup>Major axis of the ‘total’ aperture. <sup>c</sup>FWHM of a Gaussian fit.</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 1.** False-colour (R, G, B) broad-band images of NGC 4449: (a) Swift/UVOT (2486, 2221 and 1991 Å), (b) SDSS (r, g, u), (c) Spitzer/IRAC (8, 5.8 and 3.6 µm), (d) Herschel/PACS (160, 100 and 70 µm) and (e) Herschel/SPIRE (500, 350 and 250 µm). North is up, east is to the left. The images are centred at $12^{h}28^{m}11^{s}1$, $+44^{\circ}05'37''$ (J2000). The bar in the top-left image is 2 arcmin in length (2.2 kpc).</figcaption></figure>

Catalog. Colour correction was applied using interpolated correction factors suitable for a source emitting as $v^{-1.18}$ across the four WISE bands. Planck photometry was taken from column GAUFLUX(_ERR) of the Early Release Compact Source Catalogue. Similarly, the interpolated colour correction factors used were suitable for a source emitting as $v^{2.80}$ across the Planck bands centred at 350, 550 and 850 µm.

# 3 MODELLING METHOD

## 3.1 Overview

The numerical code MOCASSIN (Ercolano et al. 2003, 2005, 2008) was originally intended as a tool to construct realistic models of photoionized nebulae. It allows arbitrary three-dimensional geometries, separate for gas and dust, multiple ionizing sources emitting with a given input spectrum, variable gas chemistry, multiple dust species and arbitrary dust grain size distributions. Given these input parameters, the code self-consistently solves the radiative transfer in the coexisting gas and dust phases and calculates the ionization degree, electron temperature and dust grain temperature at every grid cell, and the overall emergent SED of the gas and dust.

In the sections that follow, we will describe how MOCASSIN (v2.02.70) can be used to model large physical systems, such as galaxies. Our model is constructed ‘bottom-up’, where as many parameters as possible are fixed a priori based on observations. The simulations are set up with the empirical elemental abundances and the observed radial gas distribution fixed. The input stellar spectrum corresponding to a mixture of stellar populations, scaled by the stellar luminosity $L_*$ and the dust to gas mass ratio ($M_{\text{dust}}/M_{\text{gas}}$) are free parameters. Based on these input parameters, which form a theoretical model of a galaxy, MOCASSIN produces a predicted low-resolution SED and a full set of predicted emission line intensities. If these predictions agree with observations, the original parameter set

<figure><img src="image_5.png" /><figcaption>**Figure 2.** Summary of the modelling method. The shaded regions group together the input variables (left) the outputs (middle) and the observational constraints (right).</figcaption></figure>

can be regarded as a true representation of the galaxy under the assumptions made. Otherwise, individual parameters can be adjusted and simulations repeated iteratively, as summarized schematically in Fig. 2.

In constructing our MOCASSIN model of NGC 4449 we have made simplifying assumptions about the geometry, gas density distribution, elemental abundances, dust composition and dust grain size distribution.

We have also made assumptions about the SFH. In Section 3.4, we discuss a general picture of the stellar populations in NGC 4449, obtained using the spectral decomposition code STARLIGHT. In Section 3.7.1, we characterize the youngest stellar population based on observational constraints. Finally, in Section 3.7.3, we fit the older populations taking into account the constraints from earlier sections.

The fitting of the SFH and dust is performed almost simultaneously for self-consistency. Therefore, the young stellar population,
the episodes of star formation and dust are never modelled on their
own. The youngest stellar population is fitted together with a repre-
sentative sample of older populations, and the final dust content is
determined relatively early in the process to ensure self-consistency.

The assumptions listed below form basis of the iterative scheme
described in detail in Section 3.7.

## 3.2 Gas density distribution

Spherical symmetry is assumed. Although MOCASSIN is a fully three-
dimensional code, this capability is not exploited at this time due
to technical limitations (see Section 3.6). However, the method
presented in this paper can be very easily adapted for fully three-
dimensional modelling and the first step towards such modelling is
discussed in Section 4.4.

We used the H I radial profile of Swaters et al. (2002) to approx-
imate n<sub>H</sub>(r)/cm<sup>-3</sup> as four second-order polynomials in the form
ar<sup>2</sup> + br + c for radii up to 3.3 kpc (3.0 arcmin), corresponding to
the extent of the galaxy in the FIR (Fig. 1). Although NGC 4449 is
highly irregular, for simplicity we treat the star formation regions
cumulatively, and based on the optical image (Fig. 1) we assume
that all gas is ionized within a distance of ~0.4 kpc from a sin-
gle ionizing source. The total observed H I mass for NGC 4449
is 2–2.5 × 10<sup>9</sup> M<sub>☉</sub> (Bajaja, Huchtmeier & Klein 1994; Swaters
et al. 2002), which is distributed between a system of streamers of
mass 9 × 10<sup>8</sup> M<sub>☉</sub> and the central region with diameter 24–35 kpc
and mass 1.1–1.25 × 10<sup>9</sup> M<sub>☉</sub> (Hunter et al. 1998, 1999). The H I
mass within the radius of 3.3 kpc derived using our formulation
is 5.5 × 10<sup>8</sup> M<sub>☉</sub>, which is consistent with the estimated total of
1.1–1.25 × 10<sup>9</sup> M<sub>☉</sub> for the inner system given the radial surface
density profile of Swaters et al. (2002).

The mass of H₂, based on $L_{CO(1-0)} = 8.4 \times 10^6$ K km s⁻¹ pc² (corrected for D = 3.8 Mpc; Böttner et al. 2003) and assuming $\alpha_{CO} \sim 10-20$ M☉ (K km s⁻¹ pc²)⁻¹ (cf. Sandstrom et al. 2012; Schruba et al. 2012 and references therein) can be estimated at $\sim 0.8-1.7 \times 10^8$ M☉. However, we note that CO-to-H₂ conversion factors for low metallicity objects are very uncertain, and therefore we include no estimates of H₂ mass in our gas mass totals.

We assumed that all gas exists in small clumps described by
the filling factor ϵ, which is the same at all radii. Considering a
Strömgren sphere of radius 0.4 kpc with L<sub>Hα</sub> given by Hunter et al.
(1999) and assuming that the average electron density is equal to
the central n<sub>e</sub> given by Martin (1997), we derive an initial estimate
of ϵ ~ 0.003. This is revisited in Section 3.7.1. Using similar as-
sumptions and the formalism of Barlow (1987) we estimate that the
mass of ionized gas is ~1.5 × 10<sup>6</sup> M<sub>☉</sub>.

## 3.3 Elemental abundances

Elemental abundances are assumed to be constant throughout the
galaxy. The abundances were taken from Vigroux, Stasińska &
Comte (1987), except for carbon, neon and sulphur, where the Large
Magellanic Cloud (LMC) ratios for C/O, Ne/O and S/O were as-
sumed (Russell & Dopita 1992). Although the sophisticated frame-
work offered by MOCASSIN can also be used to determine elemen-
tal abundances (Ercolano, Wesson & Bastian 2010), in this paper
MOCASSIN is not used in this context.

## 3.4 Star formation history

SFHs are difficult to constrain because of observational limitations
(e.g. instrumental sensitivity), but also because of the more funda-

mental age–metallicity degeneracy (e.g. Worthey 1994; Ferreras,
Charlot & Silk 1999; Ferreras & Silk 2003). However, a galac-
tic spectrum contains age-sensitive features (e.g. the Hβ absorp-
tion line; Worthey & Ottaviani 1997) and metallicity-sensitive fea-
tures (e.g. the [MgFe] index; González 1993; Thomas, Maraston &
Bender 2003), which can be useful in breaking this degeneracy.
At lower spectral resolutions the absorption features produced by
massive stars, for example Balmer lines, can be weakened by co-
incident emission within a starburst galaxy. However, higher order
Balmer lines are less prone to this effect because the strengths of cor-
responding nebular Balmer emission lines strongly decrease with
decreasing wavelength (González Delgado et al. 1998). Therefore,
a combined analysis of higher order Balmer lines and other tracers
may give an indication of the ages and metallicities of the dominant
underlying stellar populations in a galaxy. To obtain a range of pos-
sible ages and metallicities in NGC 4449 we used the spectral fitting
code STARLIGHT (v04; Cid Fernandes et al. 2005; Asari et al. 2007)
and the global optical spectrum acquired by Kennicutt (1992). This
analysis will provide a first indication of the SFH, which will form
basis for more detailed stellar population fitting in Section 3.7.3.

From the observed optical spectrum, sampled every 2 Å, we
selected the range 3700–4900 Å for fitting, with the contami-
nating Balmer emission lines Hβ–Hη, as well as the unresolved
[O II]λ3726,9 doublet and [Ne III]λ3869, masked out. We adopted 280
km s⁻¹ for the global recession velocity by direct measurement
from the observed spectrum, but the recession velocities are gen-
erally lower and vary as a function of position (Valdez-Gutiérrez
et al. 2002). The observed internal velocity dispersion, on the other
hand, is ~30 km s⁻¹ (Fuentes-Masip, Castañeda & Muñoz-Tuñón
2000). Since the grating selected by Kennicutt (1992) offered a
resolving power of only R ~ 900, the velocity dispersion cannot
be resolved by the observations and the code cannot take advan-
tage of the more suitable high-resolution stellar spectra for fitting.
Nevertheless, the distribution of stellar populations resulting from
spectral fitting performed here is not sensitive to the assumed inter-
nal velocity dispersion.

We used the synthetic SSP spectra of Bruzual & Charlot (2003) with R ~ 2000, which were generated with the stellar initial mass function (IMF) of Chabrier (2003) and provided as part of the STARLIGHT package. This assumed IMF is very similar to and consistent with the IMF of Kroupa, Gilmore & Tout (1991) and Kroupa, Aarseth & Hurley (2001) adopted in this work. STARLIGHT was run with 12 sets of 5–99 stellar templates, or ‘bases’, to probe the SFH at different sampling resolutions of age and metallicity. The lower and upper limits of the stellar ages were 1 Myr and 12 Gyr and 1–3 metallicities, corresponding to 0.2–1.0 Z☉, were tested at a time. One of the fits is presented in Fig. 3.

Our fits showed a significant spread in age–metallicity pairs at
high resolutions of age and metallicity. Therefore, in Table 3 we
present a summary of the fitted populations, grouped into ‘young’,
‘intermediate’ and ‘old’ and based on results from all fits. Table 3
shows that the most metal-poor populations are relatively well con-
strained and represent almost all of the stellar mass in NGC 4449.
We find no evidence for a significant mass of stars with metallicities
above 0.2 Z☉ amongst the old (≥1 Gyr old) population. Based on
our results, we estimate the first onset of star formation at approx-
imately 12 Gyr ago. Because STARLIGHT tended to select the oldest
available template for the oldest population, we emphasize that this
result is approximate and onset epochs ranging from 10 Gyr ago to
the age of the Universe should be considered equally plausible. Our
results also consistently point to an age of 400 ± 100 Myr as being
representative for the intermediate population.
<figure><img src="image_6.png" /><figcaption>**Figure 3.** One of 12 STARLIGHT fits, in the range 3700–4900 Å, to the global spectrum of NGC 4449 acquired by Kennicutt (1992). The input spectrum (black line) was decomposed using 33 different stellar ages at $Z/Z_⊙ = 0.2$. Contaminating emission lines in the observed spectrum, shown in blue, were masked out. See the text for more details.</figcaption></figure>

<table><caption>**Table 3.** Summary of results from STARLIGHT fits to the global spectrum of NGC 4449 acquired by Kennicutt (1992). The figures show the fractions of the total mass of stars contributed by populations in the given age and metallicity bin.</caption><thead><tr><th rowspan="2">Metallicity</th><th colspan="3">Mass fraction by age of stellar population</th></tr><tr><th>Young<br>(≤10 Myr)</th><th>Intermediate<br>(~10<sup>2</sup> Myr)</th><th>Old<br>(≥10<sup>3</sup> Myr)</th></tr></thead><tbody><tr><td>0.2 Z<sub>☉</sub></td><td>0.01</td><td>0.2–0.25</td><td>0.6–0.75</td></tr><tr><td>0.4 Z<sub>☉</sub></td><td>&lt;0.02</td><td>0.02</td><td></td></tr><tr><td>Z<sub>☉</sub></td><td>&lt;0.01</td><td>&lt;0.07</td><td></td></tr></tbody></table>

These results are supported by the resolved stellar population study of Annibali et al. (2008), suggesting continuous star formation activity for at least 1 Gyr, and by the results of Bothun (1986), who estimated that the underlying old stellar component has a mean age of 3–5 Gyr. A very young population of stars, on the other hand, is inferred from the presence of Wolf–Rayet stars (Martin & Kennicutt 1997).

Therefore, we assume that star formation started not later than ~4 Gyr ago and continued until very recently or continues to the present day. We also assume three continuous star formation episodes. This assumption reduces the complexity of the model, and, at the same time, hints at a possible SFH. The continuous-episode approach can be viewed as a generalization of a starburst approach, which is discussed in Section 4.1.

## 3.5 Interstellar dust

### 3.5.1 Distribution and composition

The dust distribution is assumed to follow that of the gas, with the dust to gas mass ratio (abbreviated here as DGR or $M_{\text{dust}}/M_{\text{gas}}$) initially kept constant throughout the galaxy (however, see, e.g. Bianchi 2008; Baes et al. 2010; MacLachlan et al. 2011; de Looze et al. 2012; Schechtman-Rook, Bershady & Wood 2012). We assumed two dominant dust species: amorphous carbon (Hanner 1988) and silicates (Laor & Draine 1993). To limit the number of free variables in the model, we further assume a constant mass ratio of carbonaceous dust to silicaceous dust (hereafter referred to as the carbon-to-silicate ratio) as a function of position. In what follows, we intend to illustrate explicitly the evolution of the carbon-to-silicate ratio for a galaxy continuously forming stars to allow us to fix this ratio for NGC 4449. We refer the reader to Dwek (1998), Morgan & Edmunds (2003), Galliano et al. (2008a) and Dwek &

Cherchneff (2011) for more complete studies of dust evolution in the ISM.

The refractory dust in the ISM is believed to form in the outflows from asymptotic giant branch stars (AGB stars; see, e.g. Matsuura et al. 2009 for the LMC), in type II core-collapse supernova (CCSN) ejecta (e.g. a multi-epoch study of SN 2004et by Fabbri et al. 2011), and has also been suggested to form in molecular clouds (Draine 2009). Dwek (1998) has suggested that low-mass AGB stars are the main source of carbon dust, while type II SNe are the main source of silicate dust. Thus, the relative abundance of carbon and silicate dust depends on the dynamics of dust formation by AGB stars and SNe, and on dust destruction processes and thus on the lifetimes of individual dust species.

Our calculations assume the dust yields of Ventura et al. (2012) for AGB stars, a constant and continuous star formation activity (Section 3.4, above), a constant metallicity, the IMF of Kroupa et al. (1991, 2001) for stellar masses $1.5 M_☉ < M < 40 M_☉$, and ignore dust destruction processes including consumption by ongoing star formation.

It is informative to first consider a flat distribution of yields across all SN masses and the carbon-to-silicate ratio inferred for SN 1987A by Matsuura et al. (2011, Model 1). Fig. 4 shows the carbon-to-silicate ratio for a galaxy actively forming stars for a continuous period of 4 Gyr. Since the data of Ventura et al. (2012) do not include AGB stars of mass lower than $1.5 M_☉$, our models cannot be reliably extrapolated beyond the first 4 Gyr. At early times, the carbon-to-silicate ratio is fixed by the assumed dust composition of SN ejecta. When the SN dust dominates the global dust budget, the ratio is not expected to change significantly with time (Fig. 4, left, red solid line). If AGB dust dominates dust production, the carbon-to-silicate ratio initially decreases after ~40 Myr due to the predicted injection of silicates by super-AGB stars. After ~400 Myr, less massive AGB stars begin injecting carbonaceous dust leading to carbon-dominated ISM dust after ~1 Gyr (Fig. 4, left, blue dotted line).

A more realistic model includes variations in dust mass production as a function of CCSN progenitor mass, based on CCSN element yields of Woosley & Weaver (1995). Since the metallicity assumed by Ventura et al. (2012) is $Z/Z_☉ = 0.05$, we chose the closest $Z/Z_☉ = 0.1$ $^{56}$Ni-producing models from Woosley & Weaver (1995) for all stellar progenitor masses; thus, ignoring fall-back of heavy elements (Moriya et al. 2010). We further scaled the CCSN elemental yields by a constant factor of 0.2 to bring the predicted dust masses into approximate agreement with the dust masses inferred from a range of SNe (Temim et al. 2006; Barlow et al. 2010; Wesson et al. 2010; Matsuura et al. 2011; Gomez et al. 2012). This scaling factor may be physically interpreted as the adopted condensation efficiency of dust. As amorphous olivine is observed to dominate the silicate mass in the diffuse ISM (Kemper, Vriend & Tielens 2004), we assumed silicate stoichiometry of $\text{MgFeSiO}_4$ for calculations of silicate dust masses.

The results in Fig. 4 (right) show that the carbon-to-silicate ratio increases with time, owing to more carbon-rich dust being produced by lower mass SNe. This is strengthened at later times by carbon-rich dust contribution from the less massive AGB stars.

In this model a single metallicity of $Z/Z_☉ \sim 0.1$ was assumed for the dust. However, the yields of Woosley & Weaver (1995) suggest that as metallicity increases, the relative production of carbon by CCSNe also increases. The opposite trend is suggested for AGB stars. In low metallicity environments AGB stars are expected to produce less silicates, while maintaining the same or an enhanced production of carbon dust (Sloan et al. 2006, 2012) than in higher
<figure><img src="image_7.png" /><figcaption>**Figure 4.** A simple model of the evolution of the global relative carbon to silicate dust content of NGC 4449, assuming continuous star formation and a constant metallicity $Z/Z_\odot \sim 0.1$. The left-hand panel shows the evolution of the carbon-to-silicate ratio for SN dust yields identical to those derived by Matsuura et al. (2011) for SN 1987A (solid line), adopting their lower limit of $M_{\text{dust,CCSN}} = 0.4 M_\odot$ and a carbon-to-silicate ratio of 0.21 (their model 1). Similar models with dust masses scaled to 0.01 $M_\odot$ (dashed line) and 0.001 $M_\odot$ (dotted line) are also presented, where the 0.001 $M_\odot$ model may be viewed as one where CCSNe are not important in global dust production. The right-hand panel shows the evolution of the carbon-to-silicate ratio for SN dust yields based on the elemental yields of Woosley & Weaver (1995, solid line). The thin solid and dashed lines show, respectively, the relative masses of carbon and silicate dust produced.</figcaption></figure>

metallicity environments. Since these trends operate on different time-scales, their combined effect is likely to result in a carbon-to-silicate profile which is more steeply increasing than the one shown in Fig. 4.

Therefore, given our assumptions about metallicity, the limited period of modelled star formation activity, and the uncertainty of the degree of dust enrichment by SNe, our final carbon-to-silicate mass ratio of 0.23 (Fig. 4, right) is likely to represent the lower limit of relative carbon content in a continuous star formation scenario. NGC 4449, with $\log(O/H) + 12 = 8.23$, is similar in terms of metallicity to the LMC (Russell & Dopita 1992), whose dust content suggests a slightly higher carbon-to-silicate mass ratio of 0.33 (Weingartner & Draine 2001). If our model is allowed to evolve beyond the first 4 Gyr, the relative carbon content in NGC 4449 would be expected to match or exceed the LMC value, depending on the age of NGC 4449. This observation may warrant the use of metallicity-based assumptions about the carbon-to-silicate ratio in future studies. Therefore, given that our spectral fitting results discussed in Section 3.4 point to an early onset of star formation, we adopt the LMC carbon-to-silicate mass ratio of 0.33.

### 3.5.2 Grain size distribution

While the absence of distinct silicate emission features in the *Spitzer* Infrared Spectrograph (IRS) spectra (Paper II) may indicate a small abundance of Si, our chemical evolution models (Section 3.5.1, above) predict a non-negligible mass of silicate dust. If all dust is assumed to reside in a thin layer around the central ionizing source, this would impose an upper limit of 20 per cent on the mass fraction of Si in ‘ultrasmall’ silicate grains ($< 15 \AA$; Li & Draine 2001) and suggest that the contribution of hot silicate dust to the global dust budget is relatively small (e.g. Smith et al. 2010). However, NGC 4449 is a face-on irregular galaxy and it is reasonable to assume, to a first-order approximation, that the radial distribution of dust follows the radial distribution of gas. Therefore, the absence of distinct silicate emission features in the global IRS spectrum may indicate effects resulting from geometry, for example extinguishing the emission from small grains by self-shielding.

We assume that the dust grain size distribution follows $n(a) \propto a^{-3.5}$ (Mathis, Rumpl & Nordsieck 1977) and we adopt the same range of grain sizes $a$ between 0.005 and 0.25 µm for both amorphous carbon and silicates.

### 3.5.3 Polycyclic aromatic hydrocarbons

Prominent emission features in the mid-infrared (IR) also known as the unidentified IR bands, are visible in the spectra of many types of objects, and are also detected in the *Spitzer* spectrum of NGC 4449 (Paper II). They are believed to originate from the C–H and C–C vibrations of large organic molecules, or PAHs (e.g. Allamandola, Tielens & Barker 1985; Peeters et al. 2002; Bauschlicher, Peeters & Allamandola 2009). The relative strengths of these emission features are believed to reflect the particle size distribution and the local physical environment of the PAH molecules (e.g. Galliano et al. 2008b; Bauschlicher et al. 2009).

We assumed the mathematical description of pure ionized PAHs from Draine & Li (2007) with the grain size distribution of Weingartner & Draine (2001) for grain sizes of 3.5–30 Å. The PAH masses in this approach are not constrained by abundances (cf. Zubko, Dwek & Arendt 2004). Since MOCASSIN currently allows only one grain size distribution per simulation, the PAHs have to be modelled separately.

## 3.6 Numerical setup

A series of preliminary simulations showed that the best compromise between high spatial resolution and reasonable computing time is obtained for a grid of 80 × 80 × 80 cells, corresponding to 0.5 million grid cells. In a spherically symmetric setup, MOCASSIN uses a Cartesian $x, y, z$ grid, with cells populating one octant of a sphere resulting in an eight-fold reduction of computing time for a given spatial resolution. The grid cells were not equal in size, and the resolutions used were 15, 50 and 185 pc for the inner 0.5 kpc region, intermediate radii up to 2.5 kpc, and for outer radii, respectively. The number of grid cells and the resolutions ensure adequate sampling of the gas density profile (Section 3.2) as a function of radius. A single ionizing source, describing the stellar content of
the entire galaxy, was placed at the origin of the coordinate system. The physical inner and outer radii of the galaxy were defined as 1 and 3300 pc.

The computation time for a single gas and dust simulation resulting from the setup described in this paper was ~9 h using eight CPUs with 8 GB of RAM per CPU.

## 3.7 Variables and convergence criteria

As described in Section 3.1 and shown in Fig. 2, the three variables, the global bolometric luminosity $L_*$, the DGR and the assumed SFH, are iteratively adjusted until all observables are well matched. Of these, $L_*$ is the most independent and can be fixed at early stages of the modelling.

The DGR and the SFH are not independent, as a higher dust content increases the extinction at the UV and optical wavelengths and therefore requires more ionizing photons to maintain the observed continuum levels. The SFH in itself is a multidimensional function, which can be explored to reach the most probable solution. We use the general picture of the SFH emerging from Section 3.4 as a constraint for constructing trial three-episode star formation scenarios. The steps described below may be viewed as an attempt to minimize a multidimensional mathematical expression by examining its partial derivatives locally.

### 3.7.1 The bolometric stellar luminosity $L_*$

A large fraction of all the energy emitted by a starburst galaxy comes from the youngest and most massive stars. These stars contribute to the observed UV and optical continuum and govern the observed global ionization states of species present in the ISM. The same radiation field is also partially absorbed by the interstellar dust and reradiated thermally in the IR. The degree of ionization of nebular species and the degree of reprocessing of radiation are also dependent on the clumpiness of the ISM.

Therefore, by focusing initially on the youngest stellar component and considering a representative set of star formation scenarios it is possible not only (i) to further simplify the SFH by constraining this youngest component, but also (ii) to fix the filling factor $\epsilon$, (iii) to fix $L_*$ and (iv) to provide an initial estimate for the DGR. Two main constraints, namely, the absolute level of the UV continuum and the ‘recent’ star formation rate (SFR) as traced by $H\alpha$ or a similar tracer, must both be matched by the total stellar luminosity and the average age of the youngest population given the adopted filling factor.

The ‘recent’ SFR for NGC 4449 is estimated at ~0.5 $M_{\odot}$ yr$^{-1}$ (Hunter et al. 1986, 1999) for a Salpeter IMF (Salpeter 1955; see discussion in Section 3.7.2). Since Dopita et al. (2006) showed that for a SSP all ionizing photons are emitted within 10 Myr, we adopt 10 Myr, the lifetime of O- and B-type stars with masses $\gtrsim 13 M_{\odot}$, as a representative period of continuous star formation for the youngest stars.

The trial SFHs (here also referred to as ‘scenarios’) were generated using STARBURST99 (Leitherer et al. 1999; v6.0.2) with the IMF of Kroupa (Kroupa et al. 1991, 2001), along with the Pauldrach/Hillier model atmospheres (Smith, Norris & Crowther 2002), the updated Padova AGB tracks (Vázquez & Leitherer 2005) and an assumed single constant metallicity of $Z/Z_{\odot} = 0.4$ (cf. Table 3).

For a representative set of scenarios consisting of three star formation episodes, we first varied the degree of clumpiness, represented by $\epsilon$, to allow enough gas to be ionized and produce the observed

nebular emission line intensities (Kobulnicky, Kennicutt & Pizagno 1999) to within a factor of a few. We found that $\epsilon = 0.033$ is most representative, which corresponds to an average $n_e$ of $41 \text{ cm}^{-3}$ within the modelled ionized region. This choice is discussed further in Section 4.3.

We constructed a grid in parameter space formed by a representative set of scenarios consisting of three star formation episodes, a range of SFRs of the youngest population which is continuously forming stars for a period of 10 Myr, and a small range of initial values for $L_*$ and the DGR. The predictions were then compared with observations. In particular,

(i) the predicted UV continuum level was compared with the observed SED (Table 2),

(ii) the estimated total number of ionizing photons, $Q(H^0)$, was compared with that derived from the integrated $H\alpha$ luminosity (Hunter et al. 1999) and

(iii) the predicted nebular emission line fluxes were compared with the measurements of Kobulnicky et al. (1999).

The best-fitting models were found by minimizing residuals between the model and the observations in a way similar to $\chi^2$ minimization. At this stage of modelling more weight was given to line intensities and to the level and shape of the continuum in the UV. To reduce the extent of the parameter space and the overall computation time, the iterative scheme in Fig. 2 was not allowed to advance automatically. Instead, after each iteration, candidate parameter sets were inspected to identify trends in each variable. The most promising parameter combinations were then manually expanded into a higher resolution parameter space and iterated according to Fig. 2.

Our results suggest an ongoing star formation activity with recent SFR $\approx 0.28 M_{\odot}$ yr$^{-1}$, $L_* \approx 5.7 \times 10^9 L_{\odot}$ and an initial DGR $\sim 1/250$; thus, allowing us to fix the youngest population and $L_*$. The empirical $L_*$, computed by integrating measurements in Table 2 between 0.15 and 3.5 $\mu$m, is $3.1 \times 10^9 L_{\odot}$. This suggests that a significant fraction of the total luminosity originates from the recently formed massive stars.

### 3.7.2 The initial mass function

The recent SFR of $0.28 M_{\odot}$ yr$^{-1}$, suggested by our initial models, is in good agreement with previous estimates (see below), when the revised distance of 3.8 Mpc (Annibali et al. 2008) and differences in the assumed IMFs are taken into account.

The Salpeter IMF, often used in the formulations of the SFR (e.g. Kennicutt 1998), is defined as $N(M)dM \propto M^{-2.35} dM$ (Salpeter 1955). However, the Kroupa IMF, which is assumed in this work, is defined in two intervals, with a flatter distribution of the number of stars at lower stellar masses. Therefore, the Salpeter IMF will overestimate the number of lower mass stars compared to the Kroupa IMF. This will in turn overestimate the total mass of stars resulting in a higher prediction of the SFR.

We generated distributions of stellar masses normalized to give the same number of stars at a stellar mass of $25 M_{\odot}$, which we assumed to be representative of the population generating the ionizing photons and giving rise to the observed $H\alpha$ emission. We computed the total mass of stars from $N(M)dM$ for $0.1 \le M/M_{\odot} \le 100$, which yielded a factor of 1.54 difference between the two IMFs. Assuming that the integrated luminosity of a galaxy scales linearly with the total mass, the SFRs based on the Salpeter IMF are therefore likely to be a factor of $\approx 1.5$ higher than if the Kroupa IMF were used. This result is also evident from the results of Dwek et al.
(2011), who tabulate the masses of all stars born per one SN event for a range of stellar IMFs.

Overall, our recent SFR of 0.28 M$_{\odot}$ yr$^{-1}$ is in good agreement with the distance-corrected and IMF-corrected estimates of 0.22–0.33 M$_{\odot}$ yr$^{-1}$ (Hunter et al. 1986), 0.30 M$_{\odot}$ yr$^{-1}$ (Hunter et al. 1999) and 0.21 M$_{\odot}$ yr$^{-1}$ (based on high ionization potential neon lines; Paper II).

### 3.7.3 Three episodes of star formation

To constrain $L_*$ and the youngest stellar population, we assumed a representative older stellar component. Since the contribution of stars older than ~10 Myr to the ionization structure, to the number of ionizing photons, or to the far-UV (FUV) continuum, is very small, the older populations can be studied in more detail without affecting the validity of our findings from Section 3.7.1.

The shape of the observed SED from the UV to the near-infrared (NIR) can be fitted by varying the lengths and characteristic epochs determining the remaining episodes of continuous star formation. Initially, we tested a set of 36 scenarios, which consisted of all combinations of:

(i) three epochs corresponding to the onset of star formation (4, 8 and 12 Gyr ago),

(ii) four epochs corresponding to the transition between the first (old) episode and the second (intermediate) episode of star formation (100, 200, 400 and 1000 Myr ago) and

(iii) three SFRs ranging from 0.05 to 0.15 M$_{\odot}$ yr$^{-1}$ further defining the second episode.

The star formation activity during the third (young) episode remained fixed in both duration and magnitude, and $L_*$ represented the luminosity of all stellar components, as described in Section 3.7.1. Therefore, the SFR of the first (oldest) episode was not an independent variable, but depended on $L_*$ and the durations and SFRs in items 1–3, for each scenario.

The resulting low-resolution SEDs were compared with the observed SED to select the best matches. Afterwards, following the scheme shown in Fig. 2, the transition epochs and the SFRs inferred from the best matches were expanded into a narrower, but higher resolution parameter space, analogous to conditions in items 2 and 3.

The criteria for finding the best-fitting combinations of stellar populations were similar to those used to find $L_*$, where we minimized residuals between the model and the observations. At this stage, more weight was given to the number of ionizing photons $Q(H^0)$, the line intensities and the level and shape of the entire stellar continuum from the UV to the NIR. As before, to reduce the extent of the parameter space and the overall computation time, the iterative scheme in Fig. 2 was not allowed to advance automatically. Instead, after each iteration, candidate parameter sets were inspected (i) to identify trends in each variable and (ii) to verify if the fit could be improved while still satisfying the observed number of ionizing photons $Q(H^0)$. The most promising parameter combinations were then manually expanded into a higher resolution parameter space.

### 3.7.4 A two-zone solution for $M_{dust}/M_{gas}$ and modelling PAHs

NGC 4449 is modelled as a centrally concentrated, spherically symmetric galaxy, even though its morphology is highly irregular (cf. Fig. 1). Therefore, spherical symmetry inevitably misrepresents the

interstellar radiation field (ISRF) within the galaxy. The single central source is much more luminous than any individual source in the real galaxy, as it must deliver enough energy to the correct mass of gas and dust to reproduce the observed strengths of the emission lines and the observed SED. Consequently, the conditions in the immediate proximity of the single central source result in unrealistically high ionization of the gaseous species and vaporization of any existing dust grains. This effect is discussed in more detail in Section 4.3.

In order to fit the dust emission in the FIR, we divided the galaxy into two zones: an inner zone with a lower DGR and an outer zone with a higher DGR, which ensures realistic absolute masses of hot dust near the centre. For technical reasons this was preferred to an alternative approach, in which the radial distribution of dust density is defined directly through the MOCASSIN input parameter $N_{dust}$. We found that for larger systems, such as galaxies, defining dust distribution through $M_{dust}$ gives best performance.

The best-fitting models were repeated with pure ionized PAH dust and the PAH grain size distribution (Section 3.5.3), and compared with the observed global SED (Table 2) and the global Spitzer/IRS spectrum (Paper II) to estimate the total mass of PAH molecules. In the final step, the resulting spectra were combined with the corresponding amorphous carbon and silicate models.

## 4 RESULTS AND DISCUSSION

### 4.1 Preliminary models

Our preliminary models included both starburst-only scenarios and continuous star formation scenarios. In the starburst-only scenario, a fit satisfying all observational constraints can be obtained for three representative starbursts occurring 3 Myr ago, 100 Myr ago and 4 Gyr ago, contributing to the total stellar mass in a ratio 1:200:1000. These models were later generalized by replacing starbursts with longer star formation episodes. For a scenario assuming two continuous star formation episodes the inferred SFRs were ~0.25 and ~0.09 M$_{\odot}$ yr$^{-1}$ between ~6 Gyr ago and ~120 Myr ago, ~120 Myr ago and the present day, respectively. The SFH in these models was constrained only by the observed SED and by the observed emission line intensities. We note that the representative populations and SFRs are not too dissimilar to the assumptions described in Section 3.4 and the recent SFR of 0.28 M$_{\odot}$ yr$^{-1}$ (Section 3.7.2). Therefore, it may be argued that combining the three convergence criteria, namely, (i) matching the observed SED, (ii) matching the total observed rate of ionizing photons, $Q(H^0)$ and (iii) matching the observed nebular emission line fluxes, may give informative results, even for a relatively simplified model.

### 4.2 Final model

Fig. 5 shows the best-fitting MOCASSIN model of NGC 4449 satisfying the observational constraints over the entire wavelength range from the UV to submillimetre. The individual SEDs arising from the assumed three continuous star formation episodes are also shown for comparison. The input parameters and the results are summarized in Table 4 and the details of the best-fitting parameters describing the SFH are given in Table 5.

The three episodes of star formation shown in Fig. 5 suggest that most of the UV emission, and most of the processed radiation
<figure><img src="image_8.png" /><figcaption>**Figure 5.** SED of the best-fitting MOCASSIN model of NGC 4449. The three star formation episodes were modelled simultaneously, but their individual SEDs are also shown for comparison. The photometric measurements have been corrected for foreground extinction, and the mid-IR spectrum acquired with Spitzer/IRS is discussed in detail in Paper II.</figcaption></figure>

<table><caption>**Table 4.** Summary of the input parameters and the results from the best-fitting MOCASSIN model of NGC 4449.</caption><thead><tr><th>Parameter</th><th>Value</th><th>Comments</th></tr></thead><tbody><tr><td colspan="3"><em>Numerical setup</em></td></tr><tr><td>Geometry</td><td>3D, spherically symmetric</td><td></td></tr><tr><td>Grid size</td><td>80 × 80 × 80</td><td>Ionizing source at (x, y, z) = (0, 0, 0)</td></tr><tr><td>Number of photons</td><td>2 × 10<sup>8</sup></td><td>Over all wavelengths</td></tr><tr><td colspan="3"><em>Observationally constrained parameters</em></td></tr><tr><td>Physical radius</td><td>3.3 kpc</td><td></td></tr><tr><td>Gas distribution</td><td>Exponential-like</td><td>Derived from the H<sub>I</sub> profile of Swaters et al. (2002)</td></tr><tr><td>M<sub>gas</sub></td><td>0.55 × 10<sup>9</sup> M<sub>☉</sub></td><td>Computed for 0 ≤ r ≤ 3.3 kpc; molecular gas not included</td></tr><tr><td>Filling factor (ε)</td><td>0.033</td><td></td></tr><tr><td>Chemical elements</td><td>H, He, C, N, O, Ne, S</td><td>Abundances constant throughout galaxy</td></tr><tr><td>Q(H<sup>0</sup>)</td><td>4.84 × 10<sup>52</sup> s<sup>-1</sup></td><td>Based on L<sub>Hα</sub> (Hunter et al. 1999)</td></tr><tr><td colspan="3"><em>Adopted dust characteristics</em></td></tr><tr><td>Dust composition</td><td>Amorphous carbon and silicates (1:3)</td><td>Hanner (1988), Laor & Draine (1993)</td></tr><tr><td>Dust grain sizes</td><td>0.005–0.25 µm (20 sizes)</td><td>Mathis et al. (1977)</td></tr><tr><td>PAH grain sizes</td><td>3.5–30 Å (10 sizes)</td><td>PAHs modelled separately, Weingartner & Draine (2001)</td></tr><tr><td colspan="3"><em>Best-fitting parameters</em></td></tr><tr><td>L<sub>*</sub></td><td>5.7 × 10<sup>9</sup> L<sub>☉</sub></td><td></td></tr><tr><td>M<sub>*</sub></td><td>1.05 ± 0.15 × 10<sup>9</sup> M<sub>☉</sub></td><td></td></tr><tr><td>M<sub>dust</sub></td><td>2.9 ± 0.5 × 10<sup>6</sup> M<sub>☉</sub></td><td>Including PAH grains</td></tr><tr><td>M<sub>PAH</sub></td><td>0.058 ± 0.005 × 10<sup>6</sup> M<sub>☉</sub></td><td>Ionized PAH grains; Draine & Li (2007)</td></tr><tr><td>M<sub>dust</sub>/M<sub>gas</sub></td><td>1/190</td><td>Effective, 1/1000 in centre, 1/182 otherwise</td></tr></tbody></table>

emitted in the FIR, originate from the intermediate-age stars of ages between 10 and ~400 Myr. The stars produced in the first episode, i.e. those older than ~400 Myr, dominate in the optical and in the NIR. The youngest stars of ages less than 10 Myr contribute mainly in the FUV.

Additionally, Fig. 6 shows the effects of varying the input parameters for the same best-fitting model (red solid line). For example,

the effect of varying $L_*$ is a vertical shift of the predicted SED in the stellar part (Fig. 6, top panel) and a corresponding change in $Q(H^0)$, as the IMF-weighted distribution of stellar luminosities is scaled by a constant.

The predicted and observed global emission line intensities, relative to $H\beta$, are shown in Table 6. The discrepancies between the predicted and the observed nebular line intensities of some species,
<figure><img src="image_9.png" /><figcaption>**Figure 6.** The effects of varying $L_*$ (top panel), the DGR (middle panel) and the relative dust composition (bottom panel) on the predicted SED. While $L_*$, the DGR and the distribution of the DGR were free variables in the models, the dust composition was fixed (Section 3.5) and its effect is presented here for illustration purposes only. The red solid lines represent the best-fitting model. See the text for more details.</figcaption></figure>

most notably [S II], result from the simplifying assumptions about the geometry and the degree of clumpiness, and are discussed in Section 4.3. Because the current version of MOCASSIN offers no treatment of PDRs, the PDR line [O I] in Table 6 is underpredicted. However, MOCASSIN is in the process of being expanded to accommodate PDRs via the recently developed three-dimensional PDR code 3D-PDR (Bisbas et al. 2012).

### 4.2.1 Stellar populations

In our three-episode model a stellar mass of $\approx 1 \times 10^9 M_{\odot}$ is produced at rates 0.25–0.10, 0.14 and 0.28 $M_{\odot}$ yr$^{-1}$ over three periods, as shown in Table 5. The best-fitting SED in Fig. 5 shows that this model does not fully account for the emission between 7000 Å and 2 µm, which may be a result of the assumption of only three
<table><caption>**Table 5.** The best-fitting three-episode SFH of NGC 4449 assuming Kroupa IMF. The assumed star formation activity is continuous, and the onsets of episode 2 and 3 coincide with the end of the preceding episodes.</caption><thead><tr><th></th><th>Episode 1<br>(Old)</th><th>Episode 2<br>(Intermediate)</th><th>Episode 3<br>(Young)</th></tr></thead><tbody><tr><td>Onset / yr ago</td><td>4–12 × 10<sup>9</sup></td><td>300–400 × 10<sup>6</sup></td><td>10 × 10<sup>6</sup></td></tr><tr><td>M<sub>*</sub> / M<sub>☉</sub></td><td>~10<sup>9</sup></td><td>~50 × 10<sup>6</sup></td><td>~3 × 10<sup>6</sup></td></tr><tr><td>SFR / M<sub>☉</sub> yr<sup>-1</sup></td><td>0.25–0.10</td><td>0.14</td><td>0.28</td></tr></tbody></table>

<table><caption>**Table 6.** Predicted and observed line strengths relative to Hβ and the predicted-to-observed ratio for the best-fitting MOCASSIN model of NGC 4449. References: (1) Kobulnicky et al. (1999) and (2) Paper II.</caption><thead><tr><th>Line</th><th>Predicted<br>[I<sub>Hβ</sub> = 1]</th><th>Observed<br>[I<sub>Hβ</sub> = 1]</th><th>Ratio</th><th>References</th></tr></thead><tbody><tr><td>[O II] λ3727</td><td>3.198</td><td>3.891</td><td>0.822</td><td>(1)</td></tr><tr><td>[Ne III] λ3868</td><td>0.173</td><td>0.189</td><td>0.916</td><td>(1)</td></tr><tr><td>Hα λ4101</td><td>0.260</td><td>0.279</td><td>0.932</td><td>(1)</td></tr><tr><td>Hγ λ4340</td><td>0.469</td><td>0.492</td><td>0.953</td><td>(1)</td></tr><tr><td>Hβ λ4861</td><td>1.000</td><td>1.000</td><td>1.000</td><td>(1)</td></tr><tr><td>[O III] λ4959</td><td>0.699</td><td>0.689</td><td>1.015</td><td>(1)</td></tr><tr><td>[O III] λ5007</td><td>2.086</td><td>2.069</td><td>1.008</td><td>(1)</td></tr><tr><td>He I λ5876</td><td>0.079</td><td>0.079</td><td>1.000</td><td>(1)</td></tr><tr><td>[S III] λ6312</td><td>0.008</td><td>0.019</td><td>0.423</td><td>(1)</td></tr><tr><td>[N II] λ6548</td><td>0.102</td><td>0.115</td><td>0.892</td><td>(1)</td></tr><tr><td>Hα λ6563</td><td>2.922</td><td>2.866</td><td>1.020</td><td>(1)</td></tr><tr><td>[N II] λ6584</td><td>0.313</td><td>0.338</td><td>0.927</td><td>(1)</td></tr><tr><td>He I λ6678</td><td>0.022</td><td>0.028</td><td>0.798</td><td>(1)</td></tr><tr><td>[S II] λ6716</td><td>1.117</td><td>0.476</td><td>2.347</td><td>(1)</td></tr><tr><td>[S II] λ6731</td><td>0.784</td><td>0.334</td><td>2.349</td><td>(1)</td></tr><tr><td>[S IV] 10.5 µm</td><td>0.264</td><td>0.219</td><td>1.205</td><td>(2)</td></tr><tr><td>[Ne III] 15.6 µm</td><td>0.228</td><td>0.387</td><td>0.588</td><td>(2)</td></tr><tr><td>[O I] 63 µm</td><td>0.012</td><td>0.413</td><td>0.030</td><td>(2)</td></tr><tr><td>[O III] 88 µm</td><td>0.654</td><td>0.647</td><td>1.011</td><td>(2)</td></tr><tr><td>[C II] 158 µm</td><td>1.112</td><td>1.307</td><td>0.851</td><td>(2)</td></tr></tbody></table>

episodes of star formation with a constant SFR for the duration of each episode. However, we note that a non-trivial three-dimensional dust distribution, where the degree of obscuration does not follow the distribution of gas or stars, may have a significant effect on the observed attenuation (e.g. Witt, Thronson & Capuano 1992; Baes & Dejonghe 2000).

We found that the results are degenerate for the oldest stars and therefore the first onset of star formation 4, 8 or 12 Gyr ago is equally plausible, resulting in a range of predicted SFRs. The spectral fits performed with STARLIGHT and discussed in Section 3.4 suggest an onset 12 Gyr ago or earlier. For an onset of star formation 12 Gyr ago, the average SFR for the duration of the first episode is 0.09 M☉ yr⁻¹. Because the mass of recently formed stars is small in comparison (cf. Tables 3 and 5), it is possible that up to 95 per cent of the total stellar mass M* in NGC 4449 was produced at higher instantaneous rates in an initial starburst.

The predicted onset of the second episode, 300–400 Myr ago, coincides with a possible encounter with the galaxy DDO 125 400–600 Myr ago (Theis & Kohle 2001). Given the simplifications of our model it would be difficult to quantify the effects of this postulated encounter on SFRs. However, we note that our STARLIGHT fits (Section 3.4) select a population ~400 Myr old as being representative for the intermediate-age populations. Similarly, the three-episode

MOCASSIN fits to the observed stellar SED converge to 300–400 Myr ago as the transition epoch between the first episode and the second.

Our model suggests that in the last 10 Myr star formation has been taking place continuously at an average rate of 0.28 M☉ yr⁻¹ (0.42 M☉ yr⁻¹ assuming Salpeter IMF; Section 3.7.2). This recent SFR, obtained by fitting both the UV continuum and the emission lines (Section 3.7.1), is consistent with other estimates based on the UV continuum emission (0.44 M☉ yr⁻¹; using data in Table 2 and following the formalism of Kennicutt 1998) or emission line intensities (Hunter et al. 1986, 1999; Paper II).

The FIR emission arises from dust heated by stars and can also be used to estimate the recent SFR. Integrating the best-fitting SED in Fig. 5 between 8 and 1000 µm yields I(FIR) ≈ 3.6 × 10⁻¹² W m⁻². The corresponding recent SFR is 0.28 M☉ yr⁻¹ assuming a Salpeter IMF (Kennicutt 1998), which is significantly lower than the estimate of 0.44 M☉ yr⁻¹ based on the UV emission. Indeed, the ratio SFR(FIR)/SFR(UV) for NGC 4449 based on the data in Table 2, and on the fit presented in Fig. 5, is 0.28/0.44 = 0.64, suggesting that approximately 35 per cent of the UV radiation does not contribute to the heating of dust. This UV ‘leakage’ (e.g. Relaño et al. 2012) is observed also in Haro 11 and NGC 4214 (Cormier et al. 2012; Hermelo et al. 2013), and may result from porosity of the ISM or different distributions for the stars and dust.

The low ratio SFR(FIR)/SFR(UV) may also suggest that, globally, in NGC 4449 the dust is heated by younger stellar populations. Similarly, a spatially resolved study by Galametz et al. (2010) suggests that the distribution of cooler dust within a different dwarf galaxy, NGC 6822, may be correlated with star formation activity. In spiral galaxies, on the other hand, several studies suggest that the dust is heated by both the younger and the evolved stellar populations (e.g. Bendo et al. 2010, 2012a; Boquien et al. 2011).

### 4.2.2 Interstellar dust

We infer a dust mass of 2.9 ± 0.5 × 10⁶ M☉ by making assumptions about the distribution of dust and its characteristics: the carbon-to-silicate ratio and the grain size distribution (cf. Section 3.5 and Table 4).

As explained in Section 3.7.4, in modelling the dust emission within NGC 4449 we adopted two zones with differing DGRs. Fig. 6 (middle panel) shows the two components, modelled simultaneously, in green and blue. The combined effect of the best-fitting DGRs and the adopted distribution of gas (Section 3.2) is a relatively constant distribution of dust as a function of radius, varying by less than a factor of 3 between the centre and r = 3.3 kpc, and yielding M<sub>gas</sub>/190 within the modelled region. The need for two zones demonstrates that the actual distribution of dust is significantly different from the distribution of gas, even for one-dimensional azimuthally averaged profiles. Fig. 6 also shows results for a smaller and a larger overall DGR, suggesting that the uncertainty in the best-fitting dust mass is ~0.5 × 10⁶ M☉.

The composition of dust has a major influence on the amount of attenuation produced at UV and optical wavelengths, and therefore on the total amount of dust required to reproduce the observed FIR emission. Fig. 5 shows that the photometric measurements agree with the model near the prominent absorption feature at 2175 Å, which is attributed to carbon dust (e.g. Fitzpatrick & Massa 1986). The Spitzer/IRS spectrum, on the other hand, is dominated by the PAH emission features and shows no discernible 10 µm silicate feature (Paper II), and no additional underlying component at 10 µm was needed to fit the PAH emission. Therefore, while
appreciating that the assumed carbon-to-silicate ratio of 1:3 may be underestimated (also see discussion in Section 3.5), we conclude that the adopted ratio is consistent with observations. We find a contribution from PAHs to the total mass of dust of 2 per cent, which is consistent with the $M_{PAH}/M_{dust}$ ratios expected for a galaxy with $1/3 Z_{\odot}$ (Galliano et al. 2008a).

For comparison, in Fig. 6 (bottom panel, dotted line) we show the SEDs corresponding to different compositions of dust for the same total mass of dust. Additional carbon-rich dust (dotted line and dashed line) produces more emission in the FIR. If such carbon-rich compositions were used in fitting the SED, the overall required mass of dust would be lower. More silicate-rich dust (dot-dashed line), on the other hand, would increase the required mass of dust. Dust mass estimates in the literature range from $2.0 \times 10^6 M_{\odot}$ (Engelbracht et al. 2008) to $3.8 \times 10^6 M_{\odot}$ (Böttner et al. 2003), reflecting different properties of dust adopted in interpreting the FIR emission.

The metallicity of NGC 4449 is similar to that of NGC 1705 and NGC 6822 studied previously (Galametz et al. 2010; O'Halloran et al. 2010). However, the DGRs, in the range 1/80–1/186, derived for these galaxies are higher. The higher ratios may result from the assumption of pure graphite or pure amorphous carbon dust, and from the adopted gas and dust budget.

We note that the DGR of 1/190 was computed using the total derived dust mass for NGC 4449, while the model covered the inner radius of 3.3 kpc, the extent of the galaxy in the FIR, and as such took into account only about 25 per cent of the total observed gas (Section 3.2). If we considered only the gas present within the modelled region, then that would require all metals in NGC 4449 to be locked into dust. However, NGC 4449 is a highly dynamic system, postulated to have been involved in a merger (Theis & Kohle 2001) and possessing a recently discovered companion (Rich et al. 2012). Since its optical body is not likely to be an isolated system, it is informative to include the extended H I envelope in the calculations, which yields a lower DGR of ~1/760. We also note that our gas masses do not include the mass of H$_2$ due to large uncertainties in the CO-to-H$_2$ conversion factors for low metallicity galaxies. Inclusion of an H$_2$ mass of ~1.7 $\times 10^8 M_{\odot}$ could decrease the global DGR even further to ~1/820.

On the other hand, an elevated flux at 850 µm, visible in Fig. 5 and not matched by the model, may suggest the presence of a significant mass of cold dust (e.g. Galliano et al. 2003). Correctly accounting for the submillimetre emission is likely to increase the total mass of dust, and consequently the DGR, and can be crucial for accurate dust mass determinations (e.g. Galametz et al. 2011).

Finally, the dust mass predicted by our model can be used to infer the fraction of metals locked into dust in NGC 4449. Adopting the total H I mass of $2.2 \times 10^9 M_{\odot}$ (cf. Bajaja et al. 1994; Swaters et al. 2002) and given the total dust mass of $2.9 \times 10^6 M_{\odot}$, the carbon-to-silicate ratio of 1:3, the LMC gas-phase C/O ratio of 0.5, the LMC gas-phase Si/O ratio of 0.29 (Russell & Dopita 1992) and the observed mean oxygen abundance O/H of $1.95 \times 10^{-4}$ (Vigroux et al. 1987), we estimate that approximately 30 per cent of carbon, 16 per cent of oxygen and 14 per cent of silicon atoms are locked in the solid phase as carbon and silicate dust. These estimates are consistent with the efficiency of 0.12 implied for SN 2003gd (Sugerman et al. 2006) and the efficiencies expected in early formed galaxies (Morgan & Edmunds 2003). They are also broadly consistent with the dust condensation efficiency of 0.2 adopted independently for the theoretical SN yields tabulated by Woosley & Weaver (1995) in order to match the dust masses inferred for a number of SNe (Section 3.5.1).

## 4.3 Spherical symmetry limitations

To achieve representative heating conditions in a spherical model, the numerical setup must allow for enough ionizing photons to interact with the correct mass of gas for a given strength of the averaged ISRF. The images in Fig. 1, as well as the *Hubble Space Telescope* observations of NGC 4449 presented by Annibali et al. (2008) and the study of its young stellar clusters by Reines, Johnson & Goss (2008) show a variety of distinct actively star-forming regions. A spherically symmetric model must reproduce not only the combined strength of the numerous individual sources, but also the average ionization structure resulting from the real three-dimensional distribution of these ionizing sources. As a consequence, the required strength of the central ionizing source may be slightly overestimated.

In general, Table 6 shows good agreement between the predicted optical and IR emission line fluxes and the observations. Since the ISRF in our models falls with distance from the central source as $1/r^2$, good agreement with observations can only be achieved by modifying the clumpiness or physical distribution of individual atomic or molecular species.

For example, our results do not correctly predict the ionization structure of sulphur: lines from the higher ionization potential $S^{2+}$ ion are systematically underpredicted, whereas lines from the lower ionization potential $S^+$ ion are overpredicted. This discrepancy cannot be explained by the assumed global abundances. However, the predicted relative sulphur line fluxes shown in Table 6 can be brought into very good agreement with observations by using a filling factor of $\epsilon = 0.055$ for sulphur in place of the lower global value of 0.033 adopted in this work. Physically, a smoother ISM shortens the mean free path of ionizing photons, increasing the degree of ionization within the modelled H II region. This result may suggest that the nebular emission in NGC 4449 originates from two or more gas phases with differing porosity, as well as from the transition zones between the phases. We conclude that a multiphase gas model has to be invoked in our spherically symmetric approach in order to reproduce the details of the observed ionization structure (cf. Cormier et al. 2012).

## 4.4 Distributed ionizing sources

A more realistic ISRF can be obtained by replacing the single central ionizing source with a uniform distribution of sources forming a central ‘ionizing sphere’. For this, we assumed a sphere with $r = 0.2$ kpc containing 100 identical sources, each contributing $L_*/100$ to the total stellar luminosity.

Fig. 6 (top panel, black solid line) presents the predicted global SED resulting from the model with distributed ionizing sources. The SED is similar to the best-fitting model (red solid line), except for the range ~10–30 µm, where the emission is weak in comparison. This discrepancy in the strength of the thermal emission by warm dust can be explained by the lower average photon energy available to interact with the dust in the case of distributed sources: dust in close proximity to one of the ionizing sources emitting with $L_*/100$ is heated to much lower temperatures than is a relatively smaller mass of dust placed near a single source emitting with $L_*$.

The plots in Fig. 7 illustrate the corresponding differences in the ionization structure as a function of radius. For example, the overall ionization fraction of $O^{2+}$ is significantly higher in the model with distributed sources, whereas the overall ionization fraction of $O^{3+}$ is significantly lower. Since in this case the resulting integrated mass of $O^{3+}$ is smaller, transitions of $O^{3+}$ are likely to produce significantly
<figure><img src="image_10.png" /><figcaption>**Figure 7.** The ionization structure of oxygen as a function of radius for the inner 1 kpc in the best-fitting MOCASSIN model of NGC 4449 (dashed lines) and the corresponding structure for the single ionizing source replaced with 100 ionizing sources distributed uniformly within a sphere of $r = 0.2$ kpc (solid lines).</figcaption></figure>

lower line intensities. The opposite would be expected for $O^{2+}$. Indeed, in the distributed sources model shown in Fig. 7, the line fluxes predicted for the [O III] lines were approximately 5 per cent higher, while the flux of the [O II]$_{λ3727}$ line was 10 per cent lower than the corresponding fluxes in the single source model (Table 6). As discussed above, the key physical factors determining individual line intensities are the mean free paths to the ionizing source and the strength of that ionizing source. Therefore, in more realistic models with distributed sources, the input luminosity $L_*$ would be expected to be lower by a few per cent to maintain good agreement between the predicted and observed line fluxes.

The physical distribution of ionizing sources affects the global ionization structure and, consequently, also affects the predicted line intensities, derived elemental abundances (cf. Ercolano et al. 2010) and may affect the overall stellar luminosity $L_*$. Therefore, we note that performing line fitting of global emission lines is a valuable constraint on the overall properties of a galaxy, but may not be meaningful without taking the morphology into consideration.

## 5 CONCLUSIONS

In this paper, we used results from previous studies, the available multiwavelength data, the spectral fitting code STARLIGHT and a simple chemical evolution code to construct a photoionization and radiative transfer MOCASSIN model of NGC 4449 through an iterative scheme.

We assumed a simplified SFH consisting of three continuous episodes. Using the shape of the observed SED, the observed rate of ionizing photons, $Q(H^0)$ derived from $L_{H\alpha}$, and the observed nebular emission line fluxes as criteria in our scheme, we infer $3 \times 10^6 M_{\odot}$ for the youngest stellar population (0–10 Myr old), which corresponds to an average recent SFR of 0.28 or 0.42 $M_{\odot}$ yr$^{-1}$, assuming a Kroupa IMF or a Salpeter IMF, respectively. We infer a bolometric stellar luminosity of $5.7 \times 10^9 L_{\odot}$ generated by a total stellar mass of $\approx 1 \times 10^9 M_{\odot}$. Although the ages of the oldest stellar populations in NGC 4449 are not well constrained by the model, we note that a very early onset, 12 Gyr ago or earlier, is selected by an independent analysis using the spectral fitting code STARLIGHT. A more in-depth study of stellar absorption line diagnostics at higher spectral resolutions would be necessary to reach more definite conclusions about these older populations.

We modelled the entire stellar component self-consistently along with a mixture of carbon and silicate dust. Our model yields a dust

mass of $2.9 \pm 0.5 \times 10^6 M_{\odot}$, which includes 2 per cent of PAHs. This estimate of $M_{\text{dust}}$ could be lower, if a higher proportion of carbon-rich dust is assumed. Overall, the results from our chemical evolution model for a galaxy continuously forming stars are consistent with the carbon to silicate dust mass ratio of 1:3 inferred for the LMC (Weingartner & Draine 2001). Interestingly, the dust condensation efficiency emerging from our MOCASSIN model and from the observed abundance of oxygen in NGC 4449 is comparable with the value of 0.2 adopted independently in our chemical evolution model to match the dust masses inferred for a range of SNe.

We note that a DGR of 1/190 was derived for the modelled region with a radius of 3.3 kpc. Including the extended H I envelope and the molecular gas is likely to lower the DGR to ~1/800. Similarly, our model does not account for cold dust, which may be responsible for an excess at 850 µm. Taking this additional component into account may have a significant effect on the derived DGR.

We conclude that our iterative scheme is a new tool, which can be used to model both the dominant stellar populations and the dust content in a self-consistent way. Although significant degeneracies in the derived parameters are expected, they can be reduced by supplementary spectroscopic data.

In the case of an irregular object, such as NGC 4449, we note that the assumption of spherical symmetry may lead to a misrepresentation of the ISRF giving rise to unphysical conditions near the central ionizing source. In particular, our results have shown that a single ISM phase does not fully reproduce the observed ionization structure, while the assumption of a single DGR overpredicts the mass of warm dust, suggesting that the radial density profiles of dust and gas are significantly different.

Our scheme is easily expandable to three dimensions, and in future work a realistic distribution of ionizing clusters could be used together with a diffuse evolved stellar population component and several gas and dust phases to construct more detailed representations of galaxies. In Paper II, we present a map of the recent SFR derived from [Ne II]$_{12.8}$ + [Ne III]$_{15.6}$, as well as spatially resolved measurements of the densities of ionized and neutral ISM phases, which may be useful in constructing a three-dimensional model of NGC 4449 in the future.

## ACKNOWLEDGEMENTS

OLK acknowledges the support of UCL's Institute of Origins. The authors acknowledge the use of the UCL Legion High Performance Computing Facility, and associated support services, in the completion of this work.

PACS has been developed by a consortium of institutes led by MPE (Germany) and including UVIE (Austria); KU Leuven, CSL, IMEC (Belgium); CEA, LAM (France); MPIA (Germany); INAF-IFSI/OAA/OAP/OAT, LENS, SISSA (Italy) and IAC (Spain). This development has been supported by the funding agencies BMVIT (Austria), ESA-PRODEX (Belgium), CEA/CNES (France), DLR (Germany), ASI/INAF (Italy) and CICYT/MCYT (Spain).

SPIRE has been developed by a consortium of institutes led by Cardiff University (UK) and including University of Lethbridge (Canada); NAOC (China); CEA, OAMP (France); IFSI, University of Padua (Italy); IAC (Spain); Stockholm Observatory (Sweden); Imperial College London, RAL, UCL-MSSL, UKATC, University of Sussex (UK) and Caltech/JPL, IPAC, University of Colorado (USA). This development has been supported by national funding agencies: CSA (Canada); NAOC (China); CEA, CNES, CNRS (France); ASI (Italy); MCINN (Spain); Stockholm Observatory (Sweden); STFC, UKSA (UK) and NASA (USA).
REFERENCES

Allamandola L. J., Tielens A. G. G. M., Barker J. R., 1985, ApJ, 290, L25

Amorín R., Pérez-Montero E., Vílchez J. M., Papaderos P., 2012, ApJ, 749, 185

Annibali F., Aloisi A., Mack J., Tosi M., van der Marel R. P., Angeretti L., Leitherer C., Sirianni M., 2008, AJ, 135, 1900

Asari N. V., Cid Fernandes R., Stasińska G., Torres-Papaqui J. P., Mateus A., Sodré L., Schoenell W., Gomes J. M., 2007, MNRAS, 381, 263

Baes M., Dejonghe H., 2000, MNRAS, 313, 153

Baes M. et al., 2010, A&A, 518, L39

Bajaja E., Huchtmeier W. K., Klein U., 1994, A&A, 285, 385

Barlow M. J., 1987, MNRAS, 227, 161

Barlow M. J. et al., 2010, A&A, 518, L138

Bauschlicher C. W., Jr, Peeters E., Allamandola L. J., 2009, ApJ, 697, 311

Bendo G. J. et al., 2010, A&A, 518, L65

Bendo G. J. et al., 2012a, MNRAS, 419, 1833

Bendo G. J., Galliano F., Madden S. C., 2012b, MNRAS, 423, 197

Bianchi S., 2008, A&A, 490, 461

Bisbas T. G., Bell T. A., Viti S., Yates J., Barlow M. J., 2012, MNRAS, 427, 2100

Boquien M. et al., 2011, AJ, 142, 111

Bothun G. D., 1986, AJ, 91, 507

Böttner C., Klein U., Heithausen A., 2003, A&A, 408, 493

Breeveld A. A., 2010, SWIFT UVOT CALDB RELEASE NOTE: AB Magnitude System (SWIFT-UVOT-CALDB-16-R01)

Breeveld A. A. et al., 2010, MNRAS, 406, 1687

Bruzual G., Charlot S., 2003, MNRAS, 344, 1000

Cardelli J. A., Clayton G. C., Mathis J. S., 1989, ApJ, 345, 245

Chabrier G., 2003, PASP, 115, 763

Cid Fernandes R., Mateus A., Sodré L., Stasińska G., Gomes J. M., 2005, MNRAS, 358, 363

Cohen M., Wheaton W. A., Megeath S. T., 2003, AJ, 126, 1090

Cormier D. et al., 2012, A&A, 548, A20

Corwin H. G., Jr, Buta R. J., de Vaucouleurs G., 1994, AJ, 108, 2128

de Looze I. et al., 2012, MNRAS, 427, 2797

de Vaucouleurs G., de Vaucouleurs A., Corwin H. G., Jr, Buta R. J., Paturel G., Fouque P., 1991, Third Reference Catalogue of Bright Galaxies. Springer-Verlag, New York

Dopita M. A. et al., 2006, ApJ, 647, 244

Draine B. T., 2009, in Henning T., Grün E., Steinacker J., eds, ASP Conf. Ser. Vol. 414, Interstellar Dust Models and Evolutionary Implications. Astron. Soc. Pac., San Francisco, p. 453

Draine B. T., Li A., 2007, ApJ, 657, 810

Dwek E., 1998, ApJ, 501, 643

Dwek E., Cherchneff I., 2011, ApJ, 727, 63

Dwek E. et al., 2011, ApJ, 738, 36

Engelbracht C. W., Rieke G. H., Gordon K. D., Smith J.-D. T., Werner M. W., Moustakas J., Willmer C. N. A., Vanzi L., 2008, ApJ, 678, 804

Ercolano B., Barlow M. J., Storey P. J., Liu X.-W., 2003, MNRAS, 340, 1136

Ercolano B., Barlow M. J., Storey P. J., 2005, MNRAS, 362, 1038

Ercolano B., Young P. R., Drake J. J., Raymond J. C., 2008, ApJS, 175, 534

Ercolano B., Wesson R., Bastian N., 2010, MNRAS, 401, 1375

Fabbri J. et al., 2011, MNRAS, 418, 1285

Ferreras I., Silk J., 2003, MNRAS, 344, 455

Ferreras I., Charlot S., Silk J., 1999, ApJ, 521, 81

Ferreras I., Cropper M., Kawata D., Page M., Hoversten E. A., 2012, MNRAS, 424, 1636

Fitzpatrick E. L., Massa D., 1986, ApJ, 307, 286

Fuentes-Masip O., Castañeda H. O., Muñoz-Tuñón C., 2000, AJ, 119, 2166

Galametz M. et al., 2010, A&A, 518, L55

Galametz M., Madden S. C., Galliano F., Hony S., Bendo G. J., Sauvage M., 2011, A&A, 532, A56

Galliano F., Madden S. C., Jones A. P., Wilson C. D., Bernard J.-P., Le Peintre F., 2003, A&A, 407, 159

Galliano F., Dwek E., Chanial P., 2008a, ApJ, 672, 214

Galliano F., Madden S. C., Tielens A. G. G. M., Peeters E., Jones A. P., 2008b, ApJ, 679, 310

Gehrels N. et al., 2004, ApJ, 611, 1005

Gomez H. L. et al., 2012, ApJ, 760, 96

González Delgado R. M., Heckman T., Leitherer C., Meurer G., Krolik J., Wilson A. S., Kinney A., Koratkar A., 1998, ApJ, 505, 174

González J. J., 1993, PhD thesis, Univ. California, Santa Cruz

Gordon K. D. et al., 2004, ApJS, 154, 215

Griffin M. J., Lim T. L., 2011, EAS Publications Series, Vol. 52, Herschel-SPIRE: An Update on Performance and Plans. Cambridge Univ. Press, Cambridge, p. 9

Griffin M. J. et al., 2010, A&A, 518, L3

Groves B., Dopita M. A., Sutherland R. S., Kewley L. J., Fischera J., Leitherer C., Brandl B., van Breugel W., 2008, ApJS, 176, 438

Guseva N. G., Papaderos P., Izotov Y. I., Noeske K. G., Fricke K. J., 2004, A&A, 421, 519

Hanner M., 1988, in Hanner M. S., ed., NASA Conf. Publ. Vol. 3004, Infrared Observations of Comets Halley and Wilson and Properties of the Grains. NASA, New York, p. 22

Hermelo I., Lisenfeld U., Relaño M., Tuffs R. J., Popescu C. C., Groves B., 2013, A&A, 549, A70

Hill R. S. et al., 1998, ApJ, 507, 179

Hoversten E. A. et al., 2011, AJ, 141, 205

Hunter D. A., Gillett F. C., Gallagher J. S., III, Rice W. L., Low F. J., 1986, ApJ, 303, 171

Hunter D. A., Wilcots E. M., van Woerden H., Gallagher J. S., Kohle S., 1998, ApJ, 495, L47

Hunter D. A., van Woerden H., Gallagher J. S., 1999, AJ, 118, 2184

James B. L., 2009, PhD thesis, University College London

Kemper F., Vriend W. J., Tielens A. G. G. M., 2004, ApJ, 609, 826

Kennicutt R. C., Jr, 1992, ApJS, 79, 255

Kennicutt R. C., Jr, 1998, ARA&A, 36, 189

Kobulnicky H. A., Kennicutt R. C., Jr, Pizagno J. L., 1999, ApJ, 514, 544

Kroupa P., Gilmore G., Tout C. A., 1991, MNRAS, 251, 293

Kroupa P., Aarseth S., Hurley J., 2001, MNRAS, 321, 699

Kuin N. P. M., Rosen S. R., 2008, MNRAS, 383, 383

Laor A., Draine B. T., 1993, ApJ, 402, 441

Leitherer C. et al., 1999, ApJS, 123, 3

Li A., Draine B. T., 2001, ApJ, 550, L213

MacLachlan J. M., Matthews L. D., Wood K., Gallagher J. S., 2011, ApJ, 741, 6

Madden S. C. et al., 2013, A&A, submitted

Martin C. L., 1997, ApJ, 491, 561

Martin C. L., Kennicutt R. C., Jr, 1997, ApJ, 483, 698

Martin D. C. et al., 2005, ApJ, 619, L1

Martínez-Galarza J. R., Groves B., Brandl B., de Messieres G. E., Indebetouw R., Dopita M. A., 2011, ApJ, 738, 176

Mathis J. S., Rumpl W., Nordsieck K. H., 1977, ApJ, 217, 425

Matsuura M. et al., 2009, MNRAS, 396, 918

Matsuura M. et al., 2011, Sci, 333, 1258

Morgan H. L., Edmunds M. G., 2003, MNRAS, 343, 427

Moriya T., Tominaga N., Tanaka M., Nomoto K., Sauer D. N., Mazzali P. A., Maeda K., Suzuki T., 2010, ApJ, 719, 1445

Morrissey P. et al., 2005, ApJ, 619, L7

Morrissey P. et al., 2007, ApJS, 173, 682

Müller T., Nielbock M., Balog Z., Klaas U., Vilenius E., 2011a, PACS Photometer - Point-Source Flux Calibration (PICC-ME-TN-037)

Müller T., Okumura K., Klaas U., 2011b, PACS Photometer Passbands and Colour Correction Factors for Various Source SEDs (PICC-ME-TN-038)

O'Halloran B. et al., 2010, A&A, 518, L58

Ott S., 2010, in Mizumoto Y., Morita K.-I., Ohishi M., eds, ASP Conf. Ser. Vol. 434, The Herschel Data Processing System - HIPE and Pipelines - Up and Running Since the Start of the Mission. Astron. Soc. Pac., San Francisco, p. 139

Peeters E., Hony S., Van Kerckhoven C., Tielens A. G. G. M., Allamandola L. J., Hudgins D. M., Bauschlicher C. W., 2002, A&A, 390, 1089

Pilgratt G. L. et al., 2010, A&A, 518, L1
Planck Collaboration, 2011, A&A, 536, A7

Poglitsch A. et al., 2010, A&A, 518, L2

Pohlen M. et al., 2010, A&A, 518, L72

Poole T. S. et al., 2008, MNRAS, 383, 627

Reines A. E., Johnson K. E., Goss W. M., 2008, AJ, 135, 2222

Relaño M., Kennicutt R. C., Jr, Eldridge J. J., Lee J. C., Verley S., 2012, MNRAS, 423, 2933

Rich R. M., Collins M. L. M., Black C. M., Longstaff F. A., Koch A., Benson A., Reitzel D. B., 2012, Nat, 482, 192

Roming P. W. A. et al., 2005, Space Sci. Rev., 120, 95

Roussel H., 2012, preprint (astro-ph/1205.2576)

Russell S. C., Dopita M. A., 1992, ApJ, 384, 508

Salpeter E. E., 1955, ApJ, 121, 161

Sandstrom K. M. et al., 2012, ApJ, preprint (astro-ph/1212.1208)

Sauvage M., 2011, Experiments in Photometric Measurements of Extended Sources (SAp-PACS-MS-0718-11)

Schechtman-Rook A., Bershady M. A., Wood K., 2012, ApJ, 746, 70

Schlegel D. J., Finkbeiner D. P., Davis M., 1998, ApJ, 500, 525

Schruba A. et al., 2012, AJ, 143, 138

Skrutskie M. F. et al., 2006, AJ, 131, 1163

Sloan G. C., Kraemer K. E., Matsuura M., Wood P. R., Price S. D., Egan M. P., 2006, ApJ, 645, 1118

Sloan G. C. et al., 2012, ApJ, 752, 140

Smith L. J., Norris R. P. F., Crowther P. A., 2002, MNRAS, 337, 1309

Smith H. A. et al., 2010, ApJ, 716, 490

Sugerman B. E. K. et al., 2006, Sci, 313, 196

Swaters R. A., van Albada T. S., van der Hulst J. M., Sancisi R., 2002, A&A, 390, 829

Temim T. et al., 2006, AJ, 132, 1610

Theis C., Kohle S., 2001, A&A, 370, 365

Thomas D., Maraston C., Bender R., 2003, MNRAS, 339, 897

Valdez-Gutiérrez M., Rosado M., Puerari I., Georgiev L., Borissova J., Ambrocio-Cruz P., 2002, AJ, 124, 3157

Valtchanov I., 2011, SPIRE Observers' Manual (HERSCHEL-DOC-0798, version 2.4)

Vasta M., Barlow M. J., Viti S., Yates J. A., Bell T. A., 2010, MNRAS, 404, 1910

Vázquez G. A., Leitherer C., 2005, ApJ, 621, 695

Ventura P. et al., 2012, MNRAS, 420, 1442

Vigroux L., Stasińska G., Comte G., 1987, A&A, 172, 15

Weingartner J. C., Draine B. T., 2001, ApJ, 548, 296

Werner M. W. et al., 2004, ApJS, 154, 1

Wesson R. et al., 2010, MNRAS, 403, 474

West A. A., Garcia-Appadoo D. A., Dalcanton J. J., Disney M. J., Rockosi C. M., Vezić Ž., Bentz M. C., Brinkmann J., 2010, AJ, 139, 315

Witt A. N., Thronson H. A., Jr, Capuano J. M., Jr, 1992, ApJ, 393, 611

Woosley S. E., Weaver T. A., 1995, ApJS, 101, 181

Worthey G., 1994, ApJS, 95, 107

Worthey G., Ottaviani D. L., 1997, ApJS, 111, 377

Wright E. L. et al., 2010, AJ, 140, 1868

Wright N. J., Barlow M. J., Ercolano B., Rauch T., 2011, MNRAS, 418, 370

York D. G. et al., 2000, AJ, 120, 1579

Zubko V., Dwek E., Arendt R. G., 2004, ApJS, 152, 211

## APPENDIX A: PHOTOMETRY OF EXTENDED SOURCES WITH Swift/UVOT

The UVOT (Roming et al. 2005) is one of three instruments on board Swift (Gehrels et al. 2004) designed to detect and observe gamma-ray bursts and their afterglows in seven optical and UV bands. To this end, the instrument is sensitive to single-photon events and the data reduction software has been tailored to point-source observations. However, with its FOV of 17 arcmin × 17 arcmin, the UVOT can also be of interest for studies of extended sources. Although this FOV is considerably smaller to that of GALEX (Martin et al. 2005), the UVOT offers a higher angular resolution, with PSFs of only 2.4–

<figure><img src="image_11.png" /><figcaption>**Figure A1.** Broad-band effective passbands for UV observations with GALEX and Swift.</figcaption></figure>

2.9 arcsec for the three UV bands (uvw2, uvm2 and uvw1) covering a comparable wavelength range to one of the two GALEX bands (Fig. A1; Morrissey et al. 2005; Breeveld et al. 2010).

Obtaining reliable photometric measurements for extended sources poses a technical challenge, which is inherent to the design of the UVOT. The signal from each incoming photon is electronically multiplied to generate a splash of photons at the CCD stage of the UVOT detector, and the centroid of the splash gives positional information to a sub-CCD-pixel accuracy. However, the photon-counting detector has a readout rate of ~90 s⁻¹, which leads to systematic undercounting ('coincidence loss' or 'pile-up') when two or more photons arrive in a similar location on the detector within one frame. In such case, the photons are not only counted as one, but the detection is also misplaced by the centroiding algorithm. At count rates of ~10 s⁻¹ this effect results in a 10 per cent loss in the number of counts (Poole et al. 2008).

High-background cases of point-source observations can be viewed as close analogues of extended sources. A range of background levels was investigated in the models of Breeveld et al. (2010), who showed that backgrounds higher than ~0.07 s⁻¹ per unbinned image pixel can no longer be fully corrected for coincidence loss without introducing an additional linear correction factor at each affected pixel. This limiting count rate was used by Hoversten et al. (2011) in their analysis of M81 to highlight regions for which coincidence loss is significant. Those regions in M81 coincided with point-like star-forming knots and could be corrected individually as point sources in the low background regime. However, NGC 4449 is a dwarf starburst galaxy, where the star formation activity, comparable to that of M81 (Gordon et al. 2004), is localized in a relatively small volume of space. The UV count rates in NGC 4449 are high throughout the galaxy, as illustrated in Fig. A2, making a similar approach not feasible.

Any method of obtaining photometry must take into account the variation of emission across the UVOT image, as areas with higher count rates suffer from a greater coincidence loss and require a greater correction. However, coincidence loss cannot be corrected for on a pixel-by-pixel basis because the counts in neighbouring pixels are not independent of each other. The position of each count
<figure><img src="image_12.png" /><figcaption>**Figure A2.** Regions affected by coincidence loss in the uvw2-band image of NGC 4449. White indicates no coincidence loss. North is up, east is to the left. The bar is 1 arcmin in length and the cross denotes the optical centre of the galaxy.</figcaption></figure>

in the image is calculated from a photon splash over five or more physical CCD pixels. Consequently, all counts within the same CCD pixel (64 image pixels) have been detected through similar photon splashes over the same area of the CCD. In the coincidence loss regime individual detections in a particular area on the detector are not independent of each other and should be considered collectively for coincidence loss correction.

## A1 Isophotal correction

We used NASA's HEASOFT (v6.7; released on 2009-08-19) to obtain corrected photometry for three UV and three optical bands in NGC 4449. The white band was excluded from this study. The task UVOTSOURCE was run repeatedly with two user-defined apertures: one enclosing a region of interest (or 'source') and one defining the background.

The aperture size recommended for UVOT photometry is 5 arcsec (Poole et al. 2008). By default, if the source aperture is greater than 5 arcsec, the coincidence loss correction factor is determined from a 5 arcsec circular region at the centre of the user-defined aperture. Consequently, for larger apertures the default procedure is likely to give a significantly biased corrected count rate, depending on features present in the central area of the user-defined aperture. Therefore, it is expected that UVOTSOURCE is generally not applicable 'as-is' to sources of angular extent greater than the size of the standard 5 arcsec aperture.

In the method presented below, the image was divided into regions following suitably chosen isophotal contours. The coincidence loss correction factors were determined from representative 5 arcsec 'test' apertures and applied to correct photon count rates in the corresponding regions. Afterwards, the background was subtracted and the count rates were summed to give the total count rate. An illustration of this isophotal setup is given in Fig. A3. Region 1 in Fig. A3 encloses areas with count rates above the threshold of 0.028 s⁻¹ pixel⁻¹ for 2 × 2 binning, where the effects of coincidence loss become non-negligible (Breeveld et al. 2010; Hoversten et al. 2011). Regions 2–4 in Fig. A3 enclose areas with count rates of 0.1,

<figure><img src="image_13.png" /><figcaption>**Figure A3.** The isophotal setup for the uvw2-band image of NGC 4449. Four isophotal regions numbered 1–4 are shown in grey, red, green and blue, respectively. Region 0, enclosed by an ellipse, and a background aperture defined as an elliptical annulus, both extend beyond the image and are not shown for clarity. The test regions are indicated by 5 arcsec circles. North is up, east is to the left. The bar is 1 arcmin in length and the cross denotes the optical centre of the galaxy.</figcaption></figure>

0.2 and 0.3 s⁻¹ pixel⁻¹. Finally, region 0 is an ellipse enclosing the entire galaxy but not aligned with its optical centre.

Measurements in all individual regions were performed by running the HEASOFT command:

```
uvotsource image=image_sk_co-add.fits \
    expfile=image_ex_co-add.fits \
    srcreg=REGION.reg bkgreg=bkg.reg \
    apercorr=NONE centroid=NO clobber=NO \
    frametime=0.0110329 output=ALL sigma=5 \
    syserr=YES outfile=photometry.fits
```

where REGION.reg is a description of (i) a test region 0–4, or (ii) one of the regions comprising regions 0–4 in Fig. A3. It is useful to define A, C<sub>r, err</sub>, C<sub>test, err</sub>, f<sub>test</sub> and C<sub>bkg, err</sub> based on the following columns in the output FITS table:

<table><tr><td>SRC_AREA</td><td>area of the source aperture in arcsec<sup>2</sup>, A</td></tr><tr><td>RAW_TOT_RATE(_ERR)</td><td>raw count rate (error) in the source aperture, C<sub>r, err</sub></td></tr><tr><td>COI_TOT_RATE(_ERR)</td><td>corrected count rate (error) in the source aperture, C<sub>test, err</sub></td></tr><tr><td>COI_STD Factor</td><td>coincidence loss correction factor in the source aperture, f<sub>test</sub></td></tr><tr><td>COI_BKG_RATE(_ERR)</td><td>corrected count rate (error) arcsec<sup>-2</sup> in the background aperture, C<sub>bkg, err</sub></td></tr></table>

The values of C<sub>test, err</sub> and f<sub>test</sub> were measured from the five test regions. The measurements from individual regions in each of regions 0–4 were combined to obtain integrated measurements for regions 0–4, A<sub>n</sub>, C<sub>r, n</sub> and C<sub>r, err, n</sub>. Since the same background aperture bkg.reg, defined as an elliptical annulus, was used in all instances, C<sub>bkg, err</sub> represents a single measurement. The total
corrected and background-subtracted count rate $C_{src}$ was then obtained by summing individually corrected stripes, or ‘isophotes’:

$$C_{src} = \sum_{n=0}^{4} C_{src,n}, \qquad (A1)$$

where

$$C_{src,n} = C_{i,n} \times f_{test,n} - C_{bkg} \times A_{i,n}, \qquad (A2)$$

$$C_{i,n} = \begin{cases} C_{r,n} - C_{r,n+1} & 0 \le n \le 3 \\ C_{r,4} & n = 4, \end{cases} \qquad (A3)$$

and

$$A_{i,n} = \begin{cases} A_n - A_{n+1} & 0 \le n \le 3 \\ A_4 & n = 4. \end{cases} \qquad (A4)$$

The upper limit to the statistical error in $C_{src,n}$ can be estimated by assuming that the fractional error in the count rate $C_{i,n}$ in each isophotal region is equal to the fractional error in the corrected count rates in the corresponding test region given by $C_{test, err, n}/C_{test, n}$. The errors associated with the corrected count rates for higher correction factors are significantly higher than the Poisson error, as shown by Kuin & Rosen (2008). Both $C_{test, err}$ and $C_{bkg, err}$ are binomial errors associated with the corrected count rates (Kuin & Rosen 2008; Poole et al. 2008).

## A2 Results

Detailed measurements in band *uvw2* are presented in Table A1. These results show that even in region 1 the average count rate per pixel is 0.055 s⁻¹, which is significantly higher than the coincidence loss threshold of 0.028 s⁻¹ (Hoversten et al. 2011).

The total count rates in all six UVOT bands are given in Table A2. The corresponding fluxes were obtained using the calibrations of Breeveld (2010), and were subsequently corrected for foreground extinction using $E(B - V) = 0.019$ (Schlegel et al. 1998) and the extinction law of Cardelli et al. (1989) with $R_V = 3.1$.

The uncertainties in the count rates listed in Table A2 combine (i) the binomial errors in each of the test regions 0–4 (Kuin & Rosen 2008) scaled to the area of the corresponding region, (ii) the estimated uncertainties due to high background (Breeveld et al. 2010) and (iii) the uncertainty resulting from the choice of thresholds for the five regions. For the UV and optical bands, the uncertainties in (ii) were estimated at 6 and 8 per cent, respectively. The uncertainties in (iii) were found to be at most 2 per cent in all bands. Thus, the

<table><caption>**Table A1.** Detailed photometric measurements for the five isophotal regions in NGC 4449 in band *uvw2* shown in Fig. A3. $C_{bkg} = 9.08 \pm 0.03 \times 10^{-4} \text{ s}^{-1} \text{ arcsec}^{-2}$. See Section A1 for column definitions.</caption><thead><tr><th>n</th><th>A<sub>i,n</sub><br>(arcsec<sup>2</sup>)</th><th>f<sub>test, n</sub></th><th>C<sub>i, n</sub><br>(s<sup>-1</sup>)</th><th>C<sub>src, n</sub><br>(s<sup>-1</sup>)</th></tr></thead><tbody><tr><td>0</td><td>130000</td><td>1.003</td><td>525.6 ± 1.4</td><td>409 ± 21</td></tr><tr><td>1</td><td>8580</td><td>1.034</td><td>470.2 ± 1.1</td><td>479 ± 5</td></tr><tr><td>2</td><td>2870</td><td>1.086</td><td>396.8 ± 0.9</td><td>428 ± 3</td></tr><tr><td>3</td><td>747</td><td>1.134</td><td>181.5 ± 0.5</td><td>205 ± 1</td></tr><tr><td>4</td><td>426</td><td>1.494</td><td>206.8 ± 0.2</td><td>309 ± 1</td></tr></tbody></table>

<table><caption>**Table A2.** Final count rates $C_{src}$ and fluxes $F_{src}$ in six UVOT bands for NGC 4449. $f_{max}$ gives the maximum correction factor used for each band. For comparison, given on the right are the fluxes in the ‘as-is’ approach. All fluxes have been corrected for foreground extinction.</caption><thead><tr><th rowspan="2">Band</th><th rowspan="2">f<sub>max</sub></th><th colspan="2">This work</th><th rowspan="2">'as-is'<br>F<sub>src</sub><br>(mJy)</th></tr><tr><th>C<sub>src</sub><br>(s<sup>-1</sup>)</th><th>F<sub>src</sub><br>(mJy)</th></tr></thead><tbody><tr><td>uvw2/1991 Å</td><td>1.49</td><td>1830 ± 120</td><td>175 ± 12</td><td>172 ± 12</td></tr><tr><td>uvm2/2221 Å</td><td>1.32</td><td>1140 ± 80</td><td>189 ± 14</td><td>202 ± 14</td></tr><tr><td>uvw1/2486 Å</td><td>1.43</td><td>1870 ± 120</td><td>202 ± 15</td><td>207 ± 14</td></tr><tr><td>u/3442 Å</td><td>1.51</td><td>3480 ± 300</td><td>249 ± 22</td><td>468 ± 40</td></tr><tr><td>b/4321 Å</td><td>1.40</td><td>4670 ± 400</td><td>463 ± 41</td><td>1200 ± 100</td></tr><tr><td>v/5410 Å</td><td>1.18</td><td>2270 ± 210</td><td>611 ± 57</td><td>1250 ± 100</td></tr></tbody></table>

statistical and systematic uncertainties in the final fluxes (Table A2) amount to ~7 and ~9 per cent overall for the UV and the optical bands.

Table A2 also gives total fluxes calculated by UVOTSOURCE directly from region 0 in the ‘as-is’ approach, i.e. ignoring the spatial variations in the correction factors across UVOT images. In bands *uvw2*, *uvm2* and *uvw1* the fluxes agree to within 7 per cent, whereas in the optical bands the discrepancy is more significant. The agreement in the UV bands may suggest that the centre of region 0 probed the global average of the count rate distribution for the UV bands, weighted by the extent of coincidence loss. However, in sources that are generally fainter in the UV, better agreement in these bands may be expected as the required correction factors are smaller. Although the values of $f_{max}$ listed in Table A2 show that the maximum correction factors were high in all bands, the emission in the UV is sharply peaked compared to the more uniform emission in the optical. This suggests that applying corrections to the UV bands is likely to be affected by a smaller systematic uncertainty than applying similar corrections to the optical bands. In general, these results confirm that UVOTSOURCE is not applicable ‘as-is’ to photometric measurements of extended sources, as expected from Section A1.

The UV and optical fluxes in Table 2 are consistent with those obtained from *GALEX* and SDSS. A plot of the global SED of NGC 4449 in the range 1000–7000 Å is presented in Fig. A4.

In the UV, the fluxes derived from *GALEX* and *Swift/UVOT* agree very well, but the absolute near-UV (NUV) flux is slightly lower than the corresponding fluxes obtained from *Swift/UVOT*. It should be noted that *GALEX*, in a similar way to *Swift/UVOT*, also suffers from local non-linearities near bright sources. These are more difficult to quantify and are not automatically corrected for (Morrissey et al. 2007). The count rates measured for NGC 4449 in a 3 arcmin diameter aperture (~1000 s⁻¹ in FUV and ~4000 s⁻¹ in NUV) suggest that this effect is likely to lower the measured count rates and add a significant systematic uncertainty to the derived *GALEX* fluxes. Therefore, the true *GALEX* FUV and NUV fluxes may be higher and may be associated with larger uncertainties than those shown in Fig. A4.

## A3 Conclusions

In the method presented above, we used the existing understanding of coincidence loss for point-source observations and applied it to obtain photometric measurements for an extended source,
<figure><img src="image_14.png" /><figcaption>**Figure A4.** Global photometry of NGC 4449 with *GALEX*, *Swift*/UVOT and SDSS using data from Table 2.</figcaption></figure>

NGC 4449, in three UV and three optical *Swift*/UVOT bands. The derived fluxes are in good agreement with the available global *GALEX* and SDSS photometry, with the overall uncertainties estimated at 7 and 9 per cent for the UV and the optical bands.

Extended-source observations with the three narrow UV bands and the small PSFs of *Swift*/UVOT enable more detailed studies of the youngest stellar component in galaxies. Although the uncertainties in our method are comparable with those of *GALEX*, we note that at higher count rates *GALEX* also suffers from non-linearity, which is not corrected for or included in the published uncertainties.

This paper has been typeset from a *TEX/LATEX* file prepared by the author.