THE STAR FORMATION HISTORY OF THE VERY METAL-POOR BLUE COMPACT DWARF I Zw 18 FROM HST/ACS DATA^*

The blue compact dwarf (BCD) galaxy I Zw 18 is one of the most intriguing objects in the local universe and has fascinated generations of astronomers since its discovery (Zwicky 1966). With a metallicity between 1/30 and 1/50 Z_☉ (Searle & Sargent 1972; Lequeux et al. 1979; Davidson & Kinman 1985; Dufour et al. 1988; Pagel et al. 1992; Skillman & Kennicutt 1993; Kunth et al. 1994; Stasińska & Leitherer 1996; Garnett et al. 1997; Izotov & Thuan 1998) it holds the record of the second lowest metallicity and lowest helium content measured in a star-forming galaxy (Izotov et al. 2009). The dynamical mass of I Zw 18 measured at a radius of 10''–12'' is ∼2–3 × 10⁸ M_☉ (Petrosian et al. 1997; van Zee et al. 1998; Lelli et al. 2012). A large amount of gas, corresponding to ∼70% of the total mass, is detected all around the system, but only 10⁷ M_☉ of H i is associated with the optical part of the galaxy (Lequeux & Viallefond 1980; van Zee et al. 1998). The H i associated with the starburst region forms a compact rapidly rotating disk, and the steep rise of the rotation curve in the inner parts indicates that there is a strong central concentration of mass (Lelli et al. 2012). H i observations have revealed a neutral hydrogen bridge connecting I Zw 18 with a faint companion galaxy (the so-called C component or secondary body), which has also been resolved into stars with Hubble Space Telescope (HST) imaging (Dufour et al. 1996; Izotov & Thuan 2004; Aloisi et al. 2007), and a ∼13.5 kpc H i tail extending to the south of the I Zw 18 main body (Lelli et al. 2012).

I Zw 18 shows very blue colors, U − B = −0.88 and B − V = −0.03 (van Zee et al. 1998), suggesting the presence of a very young stellar population, with a current star formation rate (SFR) much higher than the past mean value (Searle & Sargent 1972). All these observational characteristics make I Zw 18 resemble a primeval galaxy in the nearby universe. In fact, soon after its discovery, the question arose whether I Zw 18 is so metal poor because (1) it started forming stars only recently, so that they have not had much time to pollute metals in the interstellar medium, (2) because its star formation (SF) activity, although occurring over a long period of time, has proceeded at a rate too low for an efficient chemical enrichment, or (3) because strong galactic winds have removed from the system most of the metals. The nature of I Zw 18 has important cosmological implications. If indeed some BCDs turned out to be young galaxies, their existence would support the view that SF in low-mass systems has been inhibited until the present epoch (e.g., Babul & Rees 1992). On the other hand, the lack of such primordial systems would provide strong constraints on chemical evolution and hydrodynamical models of metal-poor galaxies, e.g., on their SF regime (continuous or bursting?) or the onset of a galactic wind.

With the advent of the HST, it has been possible to resolve the individual stars in I Zw 18 and thus to characterize its evolutionary status. From HST/WFPC2 data, Hunter & Thronson (1995) and Dufour et al. (1996) argued for a continuous SF over the last 30–50 Myr. Thanks to a more sophisticated treatment of the same data set, Aloisi et al. (1999) were able to go deeper and to detect asymptotic giant branch (AGB) stars with ages of at least several hundreds of Myr. These results were subsequently confirmed by Östlin (2000) through deep HST/NICMOS imaging. Izotov & Thuan (2004) failed to detect red giant branch (RGB) stars with the Advanced Camera for Surveys (ACS) on board the HST and concluded that the galaxy is at most 500 Myr old. However, both Momany et al. (2005) and Tosi et al. (2007), from independent reanalyses of the same ACS data set, suggested that I Zw 18 should be older than at least 1–2 Gyr, since it does appear to contain RGB stars as well.

In order to shed light on the situation, we acquired in 2006 new time-series HST/ACS photometry to study Cepheid stars in I Zw 18 and pin down its distance (GO 10586; PI: Aloisi). By combining the new data with archival ones, we both identified the RGB tip (TRGB) at I₀ = 27.27 ± 0.14 mag (D = 18.2 ± 1.5 Mpc, Aloisi et al. 2007), and detected for the first time a few Cepheids whose light curves allowed us to independently derive the distance to ≈1 Mpc of accuracy (D = 19.0 ± 0.9 Mpc; Fiorentino et al. 2010; Marconi et al. 2010). The detection of RGB stars in I Zw 18 implies that it has started forming stars at least ∼1 Gyr ago, and possibly at epochs as old as a Hubble time, ruling out the possibility that it is a truly primordial galaxy formed recently in the local universe. In this paper we use synthetic color–magnitude diagrams (CMDs) to reconstruct the entire star formation history (SFH) within the reachable look-back time and to quantify the stellar mass formed at old, intermediate-age, and young epochs. A qualitative study of the spatial distribution and evolutionary properties of I Zw 18's resolved stars has been already presented by Contreras Ramos et al. (2011).

Observations and data reduction have been extensively described in previous papers (Aloisi et al. 2007; Fiorentino et al. 2010; Contreras Ramos et al. 2011). Here we only recall the main aspects. The data were collected between 2005 October and 2006 January with the ACS/WFC in 13 different epochs in F606W (∼broad V), and 12 different epochs in F814W (∼I), for total integration times of ≈27,700 s and ≈26,200 s. The single exposures per epoch per filter were reprocessed with the most up-to-date version of the ACS calibration pipeline (CALACS). Then, for each filter, we co-added all the exposures into a single image using the software package MULTIDRIZZLE (Koekemoer et al. 2002), in order to obtain two deep images in F606W and F814W, respectively. A portion of the F606W deep image showing I Zw 18's main and secondary bodies is shown in Figure 1. The pixel size of the final resampled drizzled images is 0035 (0.7 times the original ACS/WFC pixel size). We also retrieved and combined archival ACS/WFC data in F555W (∼V) and F814W (GO program 9400; PI: Thuan), taken in 2003 over a period of 11 days, into two deep master images with total integration times of ≈43,500 s and ≈24,300 s, respectively.

Zoom In Zoom Out Reset image size

Photometry was performed with DAOPHOT (Stetson 1987) in the IRAF environment⁵ for both the archival and proprietary data sets. The instrumental magnitudes were estimated via a point-spread-function (PSF)-fitting technique. The PSF was created choosing the most isolated and clean stars in the vicinity of I Zw 18's main and secondary bodies, and was modeled with an analytic Moffat function plus additive corrections derived from the residuals of the fit to the PSF stars. In order to push the photometry as deep as possible, stars were detected in a sum of the V and I images in each data set, choosing a detection threshold of 2.4 times the local background level. This procedure allowed us to recover faint objects with very blue or red colors. Aperture photometry with PHOT, and then PSF-fitting photometry with the DAOPHOT/ALLSTAR package (Stetson 1987), were performed separately on the V and I images at the position of the objects detected in the sum image. The finding procedure, the aperture photometry, and the PSF-fitting photometry were then re-iterated on the subtracted V and I images. This allowed us to recover additional faint stars showing up only after the subtraction of their brighter neighbors. The HST F606W and F555W magnitudes of the two data sets were transformed into V Johnson–Cousins, and the F814W magnitudes into I Johnson–Cousins using the transformations in Sirianni et al. (2005). The CMD presented in Aloisi et al. (2007) was obtained by cross-correlating the photometry in the two data sets requiring good matches in position (within one pixel) and magnitude (within 3.5σ), and a Daophot sharpness parameter ⩽0.5. Then, the final V and I magnitudes were obtained averaging the magnitudes in the two data sets. In this paper, we use this matched catalog with ≈2000 sources to derive the SFH of I Zw 18.

To account for the highly variable crowding within I Zw 18, we selected five regions (A, B, and C in the main body, and D and E in the secondary body). The corresponding I, V − I CMDs are shown in Figure 2. In the main body, region A is the most crowded one, and is embedded into a shell of ionized gas, while region C, at the outskirts, is the least crowded one. In the secondary body, region E is more crowded than region D. We detect 111, 481, 590, 548, and 265 sources in regions A, B, C, D, and E, over an area of ≈0.06, 0.66, 4.83, 0.17, and 0.70 kpc², respectively. However, in region C the majority of stars are concentrated into a smaller area of 0.76 kpc², and thus we will refer to this value hereafter to compute the SFR normalized to the region area (see Table 1). The 8.7 day period classical Cepheid of Fiorentino et al. (2010) is located in region C, the other two ultra–long-period Cepheids (P ∼ 130 days; see also Fiorentino et al. 2012) are found in regions C and A, while the two long period variables (P > 100 days) are in regions B and C.

Zoom In Zoom Out Reset image size

On the CMDs in Figure 2, we have overplotted the Padova 94 (Fagotto et al. 1994) stellar tracks for different masses to show the mass ranges (and thus the look-back times) sampled in the different regions. In region A, a concentration of bright stars between the 20 and 30 M_☉ tracks, with ages 7–10 Myr, and a few stars with mass as low as ∼12 M_☉ (∼20 Myr old) are detected. Lower-mass stars are lost because of the severe photometric incompleteness due to the extreme crowding. In region B, the look-back time sampled is longer than in region A, with stellar masses as low as 2 M_☉ (∼1 Gyr old) present on the CMD. However, only in the least crowded region C are RGB stars with masses <2 M_☉ and ages >1 Gyr detected. This also holds for region D in the secondary body, while in the more crowded region E only a few potential RGB stars are observed. From the comparison of the CMDs with the theoretical tracks, we also notice the clear presence of an age gradient moving from the most crowded regions to the least crowded ones. In fact, the bright part of the blue plume (I < 26 mag) is significantly more populated in regions A, B, and E than in regions C and D, and reaches brighter magnitudes. This implies more massive stars and younger SF. For a quantitative derivation of the SFH in the different regions we refer to Section 5.

To evaluate the role of incompleteness and blending in the data, we performed artificial star experiments on the original frames, following the same procedure described in Annibali et al. (2008). These tests serve to probe observational effects associated with the data reduction process, such as the accuracy of the photometric measurements, the crowding conditions, and the ability of the PSF-fitting procedure to resolve partially overlapped sources. We performed the tests for the individual F555W, F606W, and for the two F814W frames, according to the following procedure. We divided the frames into grids of cells of chosen width (50 pixels) and randomly added one artificial star per cell at each run. This procedure prevents the artificial stars from interfering with each other, and avoids biasing the experiments toward an artificial crowding not present in the original frames. The position of the grid is randomly changed at each run, and after a large number of experiments the stars are uniformly distributed over the frame. In each filter, we assign to the artificial star a random input magnitude between m1 and m2, with m1 ≈ 3 mag brighter than the brightest star in the CMD, and m2 ≈ 3 mag fainter than the faintest star in the CMD. At each run, the frame is re-reduced following the same procedure as for the real data. To account for the fact that in the real data the source detection was performed on a sum of the V+I frames, we adopted a smaller threshold than the value of 2.4 adopted for the real data. The new threshold is derived by performing tests on the original frames, and is such that the number of detections in a single frame equals the number of detections in the sum of the V+I frames with a threshold of 2.4. This procedure is somewhat approximate. However, it has the advantage of simplifying the artificial star experiments computationally, because the stellar color does not appear explicitly as an independent variable. Instead, the dependence of the completeness and photometric errors on the stellar color follows implicitly by combining the results individually obtained for different bands. This approach is sufficiently accurate for the purposes of the present paper, since we do in fact perform our photometry of the data separately in each band. Moreover, we have not found our SFH results discussed below to depend sensitively on the exact incompleteness at the faintest magnitudes.

The output photometric catalog is cross-correlated with a sum of the original photometric catalog of real stars and the list of the artificial stars added into the frame. This prevents cross-correlation of artificial stars in the input list with real stars recovered in the output photometric catalog. We simulated ∼500,000 stars for each filter. Stars with input–output magnitude >0.75 were considered lost, because such a difference implies that they fell on a real star of their same luminosity or brighter. For each magnitude bin, the completeness of our photometry was computed as the ratio of the number of recovered artificial stars over the number of added ones.

As an example, we show in Figure 3 the distribution of the artificial star output-minus-input magnitude (Δm) in F814W for the Aloisi data set. We show the different distributions relative to regions A, B, C, D, and E selected in I Zw 18, in order to characterize the photometric error as a function of magnitude and crowding. The systematic deviation from 0 of the mean Δm is due to the increasing effect of blending at fainter magnitudes: faint stars tend to be systematically recovered with brighter magnitudes because they happen to overlap with objects of comparable or brighter magnitude. The photometric errors systematically decrease from the most crowded region A to the least crowded regions C and D. At the same time, the completeness increases, as shown in Figure 4.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

5.1. Methodology

5.2. CMD Fit

5.2.1. Region A

Figure 5 shows the comparison between the observed CMD selected in the most crowded region A (left) and two synthetic CMD realizations drawn from the best-fit SFH with the Baltimore code (middle) and with the Bologna code (right). In both cases, when searching for the best-fit solution, we excluded from the χ² minimization the region below the 60% completeness limit, where the fractional uncertainty in the completeness becomes significant. In fact, due to the small spatial extension of region A (∼0.06 kpc²) and the relatively large size of the grid cells (∼0.02 kpc²) in the artificial star procedure, the total number of artificial stars injected into region A turns out to be quite low (∼3000), providing ∼60 stars per 0.25 mag bin, and implying an uncertainty of $\sqrt{N}/N\approx 13\%$ in the completeness curve. Furthermore, the observational error simulation, which is based on a small number of recovered stars (=number of injected stars × completeness), also is not well performed below this limit.

Zoom In Zoom Out Reset image size

When running the Baltimore code, we adopted a constant pixel size of 0.5 in both color and magnitude. Since the stellar tracks fail to reach the reddest supergiants in the observed CMD, we decided to treat the box at 0.8 < V − I < 2, 21.5 < I < 23 as a single pixel. On the other hand, in the Bologna procedure, the CMD was divided into a regular grid of color–magnitude cells with Δ(V − I) = 0.25 mag and ΔI = 1 mag. The choice of the pixel size, which is driven by the need to minimize statistical fluctuations without losing the information in the CMD, remains somewhat subjective. Thus, we purposely chose different grid parameters in the two codes to test if any of the results would depend sensitively on them. Both the Baltimore and the Bologna codes provided reduced χ² values less than 1 for the best-fit solutions. The absolute χ² of the fit gives a measure of how well the data are reproduced, but it is often difficult to interpret in practice. In fact, the value of χ² can only be interpreted with the help of the mathematical χ² distribution, if the data are distributed around the model with a symmetric Gaussian error distribution. This is often not the case here, because the distribution in each CMD bin is Poissonian. This approaches a Gaussian only in the limit of large numbers of stars. Furthermore, the effective number of degrees of freedom is difficult to estimate when many Hess diagram pixels are empty (such as in region A), and do not effectively contribute to constraining the SFH. We therefore prefer to assess the quality of fit in this paper through other metrics, including data-model comparisons of luminosity functions (LFs).

Figure 5 shows that both the Baltimore and the Bologna solutions well reproduce the overall morphology of the observed CMD; however, an inspection of the LFs in Figure 6 reveals some discrepancies between observations and simulations. More specifically, the Baltimore solution tends to underpredict the blue (V − I < 0.6) counts at 22 < I < 22.5, and the red (V − I ⩾ 0.6) counts at 22 < I < 23; these stars are blue and red supergiants with masses higher than ∼15 M_☉ and ages less than ∼10 Myr in their post-main-sequence (post-MS) phase (we will ignore the shortage of simulated stars at I > 25.5 mag, since this region was excluded from the fit and the completeness correction may be uncertain). Also, the Bologna solution tends to underpredict the number of red supergiants, but provides a better match to the bright blue supergiants at 22 < I < 22.5. However, this occurs at the expenses of an excess of fainter blue counts at 24 < I < 25. Concluding, we notice that both solutions tend to under-reproduce the correct ratio of brighter over fainter massive post-MS stars. As shown in Figure 7, where we have coded the mass in the simulated CMD, this problem cannot be solved with a flatter IMF, because the same masses populate at the same time both bright and faint regions in the CMD. Possible explanations are: (1) the completeness behavior as a function of magnitude is not properly taken into account; indeed the severe crowding in region A may prevent a robust determination of the completeness function through artificial star tests; (2) some of the brightest stars in the CMD are in fact blends of two or more stars or unresolved clusters rather than individual stars; however, as already discussed in Contreras Ramos et al. (2011), the colors of the reddest (V − I ≳ 1) bright objects in I Zw 18 are incompatible with the cluster hypothesis; we can also exclude the fact that the large number of bright red stars observed in the CMD is due to reddening within I Zw 18; in fact, the largest reddening value estimated by Cannon et al. (2002) in the northwest component from the Hα/Hβ flux ratio amounts to just E(V − I) ∼ 0.1; (3) the evolutionary timescales in the brightest post-MS phase for massive stars are underestimated. The same discrepancy persists if we run some simulations with the new Padova PARSEC stellar models (Bressan et al. 2012; A. Bressan et al. 2013, private communication) created with a thorough revision of the major input physics, a new treatment of the opacities consistent with the chemical composition, and a new treatment of the mass loss. It is possible that models including additional input physics, such as stellar rotation (Meynet & Maeder 2002), could provide a better match to the brightest objects in the CMD, but these models are not available at the moment for a large set of masses at I Zw 18's metallicity.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The best-fit SFHs obtained with the two different procedures are shown in Figure 8. For both the Baltimore and the Bologna code, we show the stellar mass formed and the SFR in the different age intervals. The results are also summarized in Tables 1 and 2. The values are normalized to a Salpeter's IMF in the mass interval 0.1 to 120 M_☉. The symbol > implies that the value is a lower limit. This is because the look-back time may not be large enough to probe the entire age range in some regions. The look-back time in region A is ∼20 Myr. Both solutions show that a strong SF event occurred 6–10 Myr ago with a rate as high as ≈1.5 M_☉ yr⁻¹, comparable to the high rate derived in NGC 1569 (Greggio et al. 1998). A mass of ≈6 × 10⁶ M_☉, similar to that of the largest globular clusters ever discovered (e.g., ω Cen in the Milky Way and G1 in M31) was formed during this episode. If the SFR is divided by the area sampled by region A, we obtain ∼2.5 × 10⁻⁵ M_☉ yr⁻¹ pc⁻², e.g., just a few factors lower than the rate derived in the strong star-forming region 30 Doradus in the LMC. From our best-fit SFH we cannot exclude the presence of an even more recent burst, or ongoing SF.

Zoom In Zoom Out Reset image size

Table 2. Stellar Mass (in 10⁶ M_☉) Formed in I Zw 18 at Different Epochs

Regiona	<10	10–100	0.1–1	>1	Total
	(Myr)	(Myr)	(Gyr)	(Gyr)
A	11.07	>2.22	...	...	>13.29
	[10.26]	>[1.34]	...	...	>[11.60]
B	0.20	1.14	>0.61	...	>1.95
	[0.23]	[1.46]	>[1.68]	...	>[3.37]
C	0.02	0.14	0.70	1.1	1.96
	[0.02]	[0.19]	[1.25]	[2.4]	[3.86]
MB	11.29	>3.50	>1.31	>1.1	>17.20
	[10.51]	>[2.99]	>[2.93]	>[2.4]	>[18.83]
D	0.003	0.12	0.98	0.37	1.47
	[0.004]	[0.13]	[0.87]	[0.01]	[1.01]
E	0.005	0.34	0.40	0.28	1.02
	[0.020]	[0.27]	[0.48]	[0.02]	[0.79]
SB	0.008	0.46	1.38	0.65	2.50
	[0.024]	[0.40]	[1.35]	[0.03]	[1.80]

Notes. Values are given for the Baltimore solution and, in brackets, for the Bologna solution. The symbol > is for lower limits. ^aRegion—MB: main body; SB: secondary body.

Download table as:  ASCIITypeset image

5.2.2. Region B

The comparison between the observed CMD in region B and two synthetic CMD realizations from the best-fit SFH with the Baltimore and the Bologna codes is shown in Figure 9. We excluded from the fit the area below the dotted line where the completeness is below 40%. Given the higher density of stars on the CMD than in region A, we adopted here a smaller pixel size of 0.25 in both color and magnitude with the Baltimore code, while we chose a pixel size of 0.25 in color and 0.5 in magnitude to run the Bologna code. When running the Baltimore code, we treated the box at 0.8 < V − I < 2, 21.5 < I < 23 as a single pixel (as for region A), and excluded the box at 0.8 < V − I < 2, 24 < I < 27, where thermally pulsing AGB (TP-AGB) stars are located, from the fit. The reason for this is that the TP-AGB phase is not included in the Padova 94 tracks that we use in this paper. Although more recent Padova models (see Girardi et al. 2010) provide a much improved treatment of the TP-AGB stars, this is still a complex and highly uncertain phase of stellar evolution. Moreover, a significant but unknown fraction of TP-AGB stars may be obscured by dust shells (Boyer et al. 2009) making it difficult to use such stars to constrain the SFH. On the other hand, the Bologna code was run considering the entire CMD down to the 40% completeness limit. We purposely adopted a different approach for the TP-AGB region in the two codes to be able to assess the impact of this region on the final inferred SFH. This is valid for all the other regions as well.

Zoom In Zoom Out Reset image size

As shown in Figure 9, both solutions reproduce quite well the overall morphology of the observed CMD. The agreement is reasonable also when comparing the observed and simulated LFs in Figure 10. The main discrepancies that we notice are the following: (1) both solutions tend to under-reproduce the number of bright blue stars at I ≲ 24 mag and bright red stars at 22 ≲ I ≲ 23 (we noticed an analogous problem in region A); and (2) the Baltimore solution tends to under-reproduce the faint red counts at I > 27.5 mag, while the Bologna solution over-predicts the blue counts at I > 27.5 mag; however, we recall that these regions fall below the 40% completeness limit and were excluded from the fit when searching for the best-fit solution.

Zoom In Zoom Out Reset image size

The best-fit SFHs obtained with the two different procedures are shown in Figure 11.⁶ In both solutions, the SFR decreases with increasing look-back time (∼300 Myr is the maximum look-back time reached in this region); the highest peak of SF occurred at a rate of ≈0.04 M_☉ yr⁻¹ during the last 10–20 Myr (depending on the solution; see Table 1). The total mass formed in region B over the last ∼300 Myr is ≈2–3 × 10⁶ M_☉, depending on the solution (see Table 2).

Zoom In Zoom Out Reset image size

5.2.3. Region C

As shown in Figure 2, the CMD of region C is deep enough to reach RGB stars and thus to sample the SF at epochs older than ∼1 Gyr. The comparison between the observed CMD and two synthetic CMD realizations from the best-fit SFH with the Baltimore and the Bologna codes is shown in Figure 12. Here we were less conservative than in regions A and B, and performed the χ² minimization using the CMD down to the 20% completeness limit, in order to reach faint stars in the RGB phase (I ≳ 27.5 mag). When running the Baltimore code, we adopted a pixel size of 0.25 in magnitude and color and, as for region B, we excluded from the fit the box at 0.8 < V − I < 2, 24 < I < 27 where TP-AGB stars are located. On the other hand, the Bologna code was run adopting a pixel size of 0.5 × 0.5 and 0.25 × 0.25 in the magnitude ranges 22 < I < 25 and 25 < I < 29, respectively, and without the exclusion of CMD regions down to the 20% completeness limit. Both solutions provide a good agreement with the overall CMD morphology in Figure 12 and with the LFs in Figure 13. The Bologna solution tends to exceed the red counts at 28 < I < 28.5, while the Baltimore solution tends to under-predict the red counts at the faintest (I > 28.5) magnitudes; we recall that this region, where the completeness is very uncertain, was excluded from the fit. In general, both solutions provide a good agreement with the observed counts at the RGB tip (I ∼ 27.5 mag); here the completeness is ≈70%, high enough to allow a reliable derivation of the old (>1 Gyr) SFH.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The best-fit SFHs obtained with the Baltimore and the Bologna procedures are shown in Figure 14. According to both solutions, the average SFR over the last 100 Myr was ≈10 times as large as the average rate over the 1−10 Gyr interval (see Table 1 for details). The total stellar mass formed at epochs older than 1 Gyr is ≈(0.9–1.4) × 10⁶ M_☉, depending on the solution. However, we notice that the error bars in both SFR and mass are very large at ages >1 Gyr, compatible with no SF at all. Indeed, at a metallicity as low as Z = 0.0004 the RGB color is highly degenerate with age, and our photometric errors are large. This implies that, despite the choice of a constant 0.25 dex age bin, we are not able to resolve the SF at ages >1 Gyr. Thus, we repeated the fit adopting just one age interval between 1 and 10 Gyr, while keeping the 0.25 dex bins at younger ages. The results for both the Bologna and the Baltimore methods are shown in Figure 14. The error bars associated with the mass and the SFR imply that the detection of old SF is significant at ∼2σ: with the Baltimore and the Bologna codes we find that the mass at ages older than 1 Gyr is ≈1.1 ± 0.4 × 10⁶ M_☉ and ≈2.4 ± 1.3 × 10⁶ M_☉, respectively.

Zoom In Zoom Out Reset image size

Given that the RGB tip is so close to the completeness limit, we also tested the impact of the distance modulus on the final SFH. We found that, within the distance uncertainties derived by Aloisi et al. (2007; 18.2 ± 1.5 Mpc), there are no significant effects on the final SFH.

5.2.4. Region D

The CMD of region D in the secondary body appears quite similar to that of region C in the main body (Figure 2). The photometry is deep enough to potentially reach RGB stars sampling the SF at epochs older than ∼1 Gyr. The comparison between the observed CMD and two synthetic CMD realizations from the best-fit SFH with the Baltimore and the Bologna codes is shown in Figure 15. As for region C, the χ² minimization was done using the CMD down to the 20% completeness limit, in order to reach faint stars in the RGB phase. When running the Baltimore code, a pixel size of 0.25 in both color and magnitude was adopted and the region of TP-AGB stars at 0.8 < V − I < 2, 24 < I < 27 was excluded from the fit. On the other hand, with the Bologna code we adopted a pixel size of 0.25 × 0.25 and 0.5 × 0.5 in the magnitude ranges 22 < I < 24 and 24 < I < 29, respectively, and used the entire CMD down to the 20% completeness limit for the χ² minimization. Both solutions reproduce quite well the overall CMD morphology. The agreement with the LFs in Figure 16 is quite good too, with the only exception being the region at 27.5 < I < 28.5, V − I ⩾ 0.6, where the Baltimore solution tends to over-predict the counts.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The best-fit SFHs obtained with the two different procedures are in very good agreement with each other (Figure 17). In both solutions, the highest peak of SF occurred ∼250 Myr ago at a rate of ∼5 × 10⁻³ M_☉ yr⁻¹, and the SFR has been gradually decreasing since then. The Baltimore solution is compatible with SF during the 1–10 Gyr age interval, with an upper limit of ≈10⁶ M_☉ for the stellar mass (see Figure 17), but the error bars are very large. On the other hand the Bologna solution seems to exclude the presence of significant SF at these epochs. As done for region C in Section 5.2.3, we tried to reduce the errors on the old SF by repeating the CMD fit with only one age interval between 1 and 10 Gyr, while keeping the 0.25 dex bins at younger ages. However, in this case we find, both with the Bologna and the Baltimore method, that the mass formed at ages older than 1 Gyr is still below the errors. For this reason, we cannot claim a significant detection of old SF in region D.

Zoom In Zoom Out Reset image size

5.2.5. Region E

The comparison between the observed CMD in the most crowded region E of the secondary body and two synthetic CMD realizations from the best-fit SFH with the Baltimore and the Bologna codes is shown in Figure 18. Due to the higher crowding than in region D, it is more difficult here to sample RGB stars. In fact, the completeness is as low as ∼30% at the RGB tip. The fit was done considering the CMD down to the 20% completeness limit. When running the Baltimore code, we excluded the region of the TP-AGB stars, and adopted a pixel size of 0.25 in both color and magnitude. On the other hand, when running the Bologna code, we considered the entire CMD down to the 20% completeness limit and adopted a pixel size of 0.5 × 0.5 and 0.25 × 0.25 in the magnitude ranges I < 26 mag and I ⩾ 26 mag, respectively.

Zoom In Zoom Out Reset image size

From Figures 18 and 19 we notice that the Baltimore solution tends to under-predict the blue (V − I < 0.6) counts at I > 28 mag while overproducing the red (V − I ⩾ 0.6) counts at 24.5 < I < 27.5; since this is the region of TP-AGB stars that was excluded from the fit, the poor agreement between simulations and observations should not be a surprise. On the other hand, the Bologna solution seems to better reproduce both the overall CMD morphology and the LFs. Despite these differences, the best-fit SFHs obtained with the two procedures (Figure 20) are in very good agreement with each other. The highest peak of SF occurred 30–10 Myr ago at a rate of ≈7 × 10⁻³ M_☉ yr⁻¹. This is much lower than what found for the most active regions of the main body. The Baltimore solution is compatible with the presence of SF at ages older than 1 Gyr (see Tables 1 and 2).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In previous papers (Aloisi et al. 2007; Fiorentino et al. 2010; Contreras Ramos et al. 2011) we have tried to qualitatively characterize the different stellar populations in I Zw 18 and have demonstrated the presence of RGB stars, excluding the possibility that this is a truly primordial system in the present-day universe. In this paper we use HST/ACS data in combination with the method of synthetic CMDs to quantitatively derive the SFH in I Zw 18. Previous similar studies have been presented by Aloisi et al. (1999) and by Jamet et al. (2010). However, the SFH derived by Aloisi et al. (1999) was based on shallower WFPC2 data and on the assumption of a distance of 10 Mpc, much closer than what is suggested by the deeper ACS CMD and by the Cepheid study. On the other hand, Jamet et al. (2010) used ACS data (only the Thuan data set) for the derivation of the SFH in the secondary body of I Zw 18, but inferred a distance of ∼27 Mpc from the CMD analysis, incompatible with the distance of ∼18 Mpc obtained from the Cepheid light curves. In this paper we present for the first time the SFH of I Zw 18 based on a robust distance assumption.

In order to account for the different crowding conditions in I Zw 18, we have selected three regions (A, B, and C) and two regions (E and D) in the main and secondary bodies, respectively. In the main body, region A is the most crowded one and corresponds to what has been referred to as the northwest component in previous papers (see, e.g., Figure 1 in Papaderos & Ostlin 2012), while region C is the least crowded one and is the most suitable to study the old stellar population. In the secondary body, region E follows the comma-shaped distribution of the brightest blue stars, while region D comprises the redder, less crowded "halo."

The first important result of this paper is an estimate of the stellar mass locked up in old (age >1 Gyr) stars. Our study indicates a mass of ≈1–2 × 10⁶ M_☉ in old stars in region C, which amounts to at least half of the total stellar mass in this region. This detection is significant at ≈2σ. On the other hand, the severe crowding in regions A and B prevents us from reaching RGB stars and deriving the SFH at ages older than 1 Gyr. From this result we can derive a lower limit for the total mass in old stars in the whole main body. Assuming in region A and B the same surface density of old stars as in region C, and using the areas reported in Table 1, we obtain a mass of ≈2 × 10⁶ M_☉. Since stars tend to be more concentrated toward the center, it is likely that the true mass in old stars in the main body is significantly higher than this. As for the secondary body, the detection of old SF is less significant, and we do not attempt to provide a lower limit.

The second important result of our study concerns the very recent SF in I Zw 18. The main body has been forming stars very actively in recent epochs, with an average SFR over the last 10 Myr as high as ≈1 M_☉ yr⁻¹ in the northwest component (region A). This corresponds to a rate normalized to the area of ≈2 × 10⁻⁵ M_☉ yr⁻¹ pc⁻², which is just a few factors lower than the rate derived in the strong star-forming region 30 Doradus in the LMC. We recall, however, that this result relies on the assumption that all the brightest stars in the CMD are in fact individual objects rather than blends of two or more stars or unresolved star clusters, although we have demonstrated that the colors of the reddest (V − I ≳ 1) bright objects are incompatible with the cluster hypothesis (Contreras Ramos et al. 2011). Indeed, our value can be compared with the rate inferred from the Hα luminosity. From Cannon et al. (2002) and assuming a distance of 18 Mpc, we derive a Hα luminosity of ∼5.9 × 10³⁹ erg s⁻¹ for the northwest component. Notice that this only includes the ionized gas around the peak of the continuum emission (≈region A) but not the ionized gas extending outside the stellar component. Using the relation from Kennicutt et al. (1994) for a Salpeter's IMF in the 0.1–100 M_☉ mass interval, we obtain an SFR of ∼0.05 M_☉ yr⁻¹, significantly lower than what was derived from the synthetic CMDs. We notice that the occurrence of SF at very recent epochs in the main body is in agreement with the presence of ultra–long-period variables (Fiorentino et al. 2010; Marconi et al. 2010). A potential rate of ≈1 M_☉ yr⁻¹ in region A implies a stellar mass of ≈10⁷ M_☉. This high concentration of mass is interesting in view of the results by Lelli et al. (2012), who analyzed archival Very Large Array data and found that the H i associated to the starburst region forms a compact fast-rotating disk. The rotation curve is flat with a steep rise in the inner parts, indicating the presence of a strong central concentration of mass.

In the main body, the average SFR over the last 10 Myr decreases from the most crowded region A to the least crowded region C down to ≈2 × 10⁻³ M_☉ yr⁻¹. However, this is still comparable to the average rate inferred in region C over the last ∼1 Gyr. On the other hand, the secondary body was much less active than the main body during the last ∼10 Myr, in agreement with the absence of nebular emission everywhere but in its central region (see also Papaderos & Ostlin 2012). The peak of SF was found to occur ∼15 Myr ago in the more crowded region E and ∼250 Myr ago in the less crowded region D.

• 1.  
  We confirm that I Zw 18 has started forming stars earlier than ∼1 Gyr ago, and possibly at epochs as old as a Hubble time. Thus it is not a truly young galaxy at its first bursts of SF, as argued in previous studies (e.g., Izotov & Thuan 2004; Jamet et al. 2010; Papaderos & Ostlin 2012).
• 2.  
  In the periphery of I Zw 18's main body, where crowding is low enough to potentially detect the old population, we estimate a mass of 1–2 ×  ≈ 10⁶ M_☉ in stars older than ∼1 Gyr, accounting for at least half of the total stellar mass in this region. On the other hand, crowding is too severe to allow the detection of such a population in the more central regions. Assuming there the same surface density of old stars as in the periphery, we obtain a lower limit of ≈2 × 10⁶ M_☉ for the total mass in old stars in I Zw 18's main body. The presence of a significant amount of old stars in the secondary body is more uncertain.
• 3.  
  The main body has been forming stars at a very high rate in recent epochs. In the most crowded northwest region (region A), the average rate over the last ∼10 Myr was as high as ≈1 M_☉ yr⁻¹, corresponding to a rate normalized to the area of ≈2 × 10⁻⁵ M_☉ yr⁻¹ pc⁻². However, this result relies on the assumption that all the brightest stars in the CMD are in fact individual objects rather than blends of two or more stars or unresolved star clusters.
• 4.  
  While in the main body, the peak of activity occurred during the last ∼10 Myr, the secondary body was much less active at these epochs, in agreement with the absence of significant nebular emission. In the secondary body, the peak of activity was found to occur ∼15 Myr ago in the more crowded central region and ∼250 Myr ago in the periphery.
• 5.  
  The high current SFR in I Zw 18 explains why this galaxy is so blue and has a high ionized gas content, resembling primeval galaxies in the early universe. Detailed chemical evolution models are required to quantitatively check whether the SFH from the synthetic CMDs can explain the low measured element abundances, both in the ionized (e.g., Garnett et al. 1997; Izotov & Thuan 1998) and in the neutral interstellar medium (Aloisi et al. 2003; Lecavelier des Etangs et al. 2004; James et al. 2013), or if galactic winds with loss of metals are needed, as originally suggested by Matteucci & Tosi (1985), Pilyugin (1993), and Marconi et al. (1994).

This paper is dedicated to the twins Astrid and Leonardo, born right after its submission. We thank the referee for comments and suggestions that helped us to improve the presentation of our results. F.A. has been supported by Cofis-ASI-INAF-07, ASI I009/10/0, and PRIN-INAF-2010. G.C., M.C., and M.T. have been also partially funded with the above grants. M.C., G.F., and M.T. are grateful to the International Space Science Institute (ISSI) in Bern for enlightening conversations and warm hospitality during Team meetings.