Disentangling the metallicity and star formation history of HII galaxies through tailor-made models

We present a self-consistent study of the stellar populations and the ionized gas in a sample of 10 HII galaxies with, at least, four measured electron temperatures and a precise determination of ionic abundances following the"direct method". We fitted the spectral energy distribution of the galaxies using the program STARLIGHT in order to quantify the contribution of the underlying stellar population to EW(Hbeta), which amounts to about 10% for most of the objects. We then studied the Wolf-Rayet stellar populations detected in seven of the galaxies. The presence of these populations and the corrected EW(Hbeta) values indicate that the ionizing stellar populations were created following a continuous star formation episode of 10 Myr duration, hence WR stars may be present in all of objects even if they are not detected in some of them. The derived stellar features, the number of ionizing photons and the relative intensities of the strongest emission lines were used as input parameters to compute tailored models with the photoionization code CLOUDY. Our models are able to y reproduce their thermal and ionization structure as deduced from their collisionally excited emission lines and, hence, no abundance discrepancy factors are implied for this kind of objects. Only the electron temperature of S+ is overestimated by the models, pointing to the possible presence of outer shells of diffuse gas in these objects. This kind of geometrical effects can affect the determination of the equivalent effective temperature of the ionizing cluster using calibrators which depend on low-excitation emission lines.


INTRODUCTION
The understanding of the physical properties of gas, stars, and dust and their interrelations and evolution in different scenarios of star formation constitute one of the most intriguing open issues on astrophysics. This type of phenomenon is well studied at different spatial scales, since massive young star clusters are powerful sources of radiation, detectable even in the most distant galaxies. For in-stance, from a statistical point of view, large surveys of starburst galaxies have been used to confirm, with a great confidence base, some of the observational relations between the physical properties of different types of galaxies at different cosmological epochs. This is the case of the luminositymetallicity relation or its equivalent, the mass-metallicity relation, which establishes that in the local Universe the most massive galaxies have, in average, higher metallicities than the less massive ones (SDSS: Tremonti et al., 2004; 2dFGRS: Lamareille et al., 2004). This observational fact is a consequence of higher stellar yields in the more massive galaxies, which could find less troubles to retain the metals ejected by stars and pushed out towards the interstellar and intergalactic medium by supernovae and stellar winds (Larson, 1974).
Nevertheless, since this kind of studies are mostly based on the relation between the stellar mass and the oxygen gasphase abundance, they present some drawbacks, involving several aspects in the determination of their properties, including the following: 1) The relation between the total mass of a galaxy and its stellar mass depends strongly on different environmental and morphological circumstances (Bell & de Jong, 2001). 2) The stellar mass is usually estimated by means of stellar synthesis population model fitting of the spectral energy distribution (SED). These models are often affected by degeneracies between age, metallicity, and dust-extinction which are never well solved (Renzini & Buzzoni, 1986). 3) The metallicity is only estimated in starburst galaxies through the oxygen gas phase abundance derived from the relative intensities of the most prominent emission lines. This excludes from these studies the most quiescent and early-type galaxies. 4) The spatial distribution of metallicity in massive, non-compact objects is often neglected, assigning a unique value to the whole distribution. However, it is known that the radial profile of the oxygen abundance in some close spiral galaxies can spread over more than an order of magnitude (Searle, 1971). 5) The strongline parameters used to derive metallicities are often calibrated employing sequences of models which do not follow the distribution of metallicity obtained with more confident observational methods like, for instance, the direct method (Pérez-Montero & Díaz, 2005), which is based on the previous determination of the electron temperature of the ionized gas by means of the fainter auroral emission lines.
HII galaxies usually represent the low-end of the massmetallicity relation in different surveys. In fact, the work of Lequeux et al. (1979), the first where the correlation between oxygen gas phase abundance and stellar mass was found, is based on this type of blue dwarf irregular galaxies. Their spectra are characterized by very bright emission lines while the total luminosity and color are dominated by the presence of intense bursts of recent star formation (Sargent & Searle, 1970;French, 1980). They are dwarf and compact so they seldom show chemical inhomogeneities. The metallicity distribution of HII galaxies has an average value (12+log(O/H) ∼ 8.0, Terlevich et al., 1991) sensibly lower than the rest of starburst galaxies and, in fact, the objects with the lowest metal content in the Local Universe belong to the subclass of HII galaxies (IZw18, Searle & Sargent, 1972).
Sloan Digital Sky Survey (SDSS) represents the largest sample of galaxies in the Local Universe and it has been explored for the study of the mass-metallicity relation (Tremonti et al., 2004). Unfortunately, this survey presents some drawbacks that prevent a precise determination of the metallicity in the emission lines-like objects: 1) Since the spectral coverage of the survey starts at 3800Å, the sample objects at redshift ∼ 0 have not any detection of the [Oii] λ 3727Å line and hence no direct measurement of the total abundance of oxygen is possible. Only an estimate can be given in terms of the weaker auroral lines at λλ 7319,7330 A but with a high degree of uncertainty (Kniazev et al., 2003). 2) Only a fraction of the objects show in their spectra the weak auroral emission lines necessary to estimate the electron temperature and, hence, to calculate chemical abundances following the direct method .
3) The only auroral lines seen in the most part of this subsample is [Oiii] λ 4363Å and, therefore, only the electron temperature of the high excitation zone can be derived. The total ionization structure of the nebula is then often obtained based on this temperature through a very simplified assumption (Pérez-Montero & Díaz, 2003). 4) For the rest of the objects with no direct electron temperature derivation, a strong-line method is used, so the determination of metallicity is hugely sensitive to the chosen calibration (Kewley & Ellison, 2008). 5) The 3 arc secs aperture of the SDSS fiber precludes a good determination of the physical properties of the starburst region of these galaxies.
Double-arm spectrographs are suitable tools to observe the entire spectral range between 3500Å and 1 micron in one single exposure at good spectral resolution the star forming regions of HII galaxies, avoiding the effects due to a different covering of the object in different spectral ranges. In the case of compact HII galaxies, their bursts of star formation can be covered entirely with a narrow slit. Assuming that different electron temperatures are measured in the optical spectrum, this implies that different excitation regions inside the ionized gas can be characterized, leading to a more precise determination of the ionic chemical abundances. There are some reported discrepancies between abundance determinations obtained from bright collisional lines and the fainter recombination ones (Peimbert et al., 2007), which do not depend on the determination of the electron temperature. Nevertheless, the existence of these abundance discrepancy factors (ADFs) is still not conclusively established. In fact, in the case of HII galaxies, no significant deviations have been found between the electron temperatures derived from collisional lines and from the Balmer (Hägele et al., 2006) or Paschen discontinuities . Besides, the high ADFs derived from optical recombination lines are sensibly lower than those calculated from mid-infra-red lines Lebouteiller et al., 2008) or from planetary nebulae (Williams et al., 2008 ).
In this work we model the stellar content and the abundances of the most representative ions of a sample of HII galaxies previously observed in the SDSS and re-observed with different double arm spectrographs in the spectral range between 3500Å and 1 micron (Hägele et al., 2006, hereafter H06 andHägele et al., 2008, hereafter H08). The study of the stellar population in the star forming knots of these galaxies includes underlying and ionizing populations. We show that the correct characterization of the ionizing stellar population, including Wolf-Rayet stars, and an appropriate geometry of the gas, dust content and the abundances of the main elements are enough to reproduce with precision the intensities of the observed emission lines. We demonstrate too that it is possible to draw a consistent picture of the thermal and ionization structure of the gas based on the measurement of as many electron temperatures as possible and, hence, to calculate ionic abundances with high precision using optical to far red collisional emission lines.
In the next section we describe the sample of objects that we have modeled. In section 3 we describe our results concerning the stellar content through model fitting of the SED and the analysis of the Wolf-Rayet features. We also describe the photoionization models of the gas, once the underlying stellar population has been removed. In section 4, we discuss the star formation history and the physical properties obtained from the models and we compare them with observations. In section 5, we summarize our results and present our conclusions. Finally, an appendix has been added to compare the total abundances and ionization correction factors obtained from models with those obtained in H06 and H08.

THE SAMPLE OF STUDIED OBJECTS
The modeled galaxies were selected from the SDSS using the implementation of the database in the INAOE Virtual Observatory superserver 1 (see H06 and H08 for more details) . We selected the brightest nearby narrow emission line galaxies with very strong lines and large equivalent widths of the Hα line from the whole SDSS data release available when we planned each observing run. Once AGN-like objects were removed by BPT diagnostic diagrams (Baldwin, Phillips and Terlevich, 1981), the sample was restricted to the largest Hα flux and highest signal-to-noise ratio objects (López, 2005). The HII galaxies studied in this work comprise 10 objects belonging to the final list and observed in two independent observing runs with different double-arm spectrographs in order to improve the spectral range and the signal-to-noise. Three out of them were observed with the ISIS spectrograph mounted on the William Herschel Telescope (H06) and the other 7 were observed with the TWIN instrument in the 3.5 m. CAHA telescope (H08). In Table 1 we list the galaxies, with the names adopted in this work, the spectroscopic redshift as listed in the SDSS catalogue and their extinctioncorrected Hα luminosities in erg·s −1 . Distances have been derived from the measured redshifts for cosmological parameters: H0 = 73 km s −1 Mpc −1 , Ωmatter = 0.27, Ωvaccum = 0.73. We also list the oxygen abundance as calculated in H06 and H08 and, finally, other names by which the galaxies are also known. Both observing runs had similar configurations. In WHT observations the spectral range goes from 3200Å up to 10550Å with a spectral dispersion of 2.5Å/px FWHM in the blue arm and 4.8 in the red. In CAHA observations, the spectral range covers from 3400 up to 10400Å (with 1 http://ov.inaoep.mx/ a gap between 5700 and 5800Å) with spectral resolutions of 3.2 and 7.0Å/px FWHM in the blue and the red part respectively. The observations were taken at parallactic angle. These configurations allow to detect simultaneously the bright emission lines from [Oii] at 3727Å up to [Siii] at 9532Å, allowing the derivation of five different electron tem- Siii]) and t([Nii]) ) in three objects of the sample and four temperatures (the same listed above with the exception of t ([Nii]) in the other seven). The derivation of all these electron temperatures allows a more precise determination of the thermal structure of the gas and, hence, of the ionic abundances of the detected species in the observed spectra. Our sample galaxies have oxygen abundances in the range 7.93 < 12+log(O/H) < 8.19 and equivalent widths of Hβ in the range 83Å < EW(Hβ) < 171Å.

Model fitting of the stellar population
Underlying stellar populations in starburst galaxies have several effects in their spectra that can affect the measurement of emission lines. Firstly, since this population has a weight in the galaxy continuum, most of the times the measurement of the equivalent widths of the emission lines is not a trustable estimator of the age of the ionizing stellar population . Secondly, the presence of a conspicuous underlying stellar population depresses the Balmer and Paschen emission lines and do not allow to measure their fluxes with an acceptable accuracy (Díaz, 1988). This can alter those properties dependent on the fluxes of these recombination lines, like the reddening or the ionic abundances. In the case of the helium recombination lines, this effect is present too, affecting the derivation of the helium abundance (Olive & Skillman, 2004).
We have subtracted from the observed spectra the spectral energy distribution of the underlying stellar population found by the spectral synthesis code STARLIGHT 2 (Cid Fernandes et al. 2004Fernandes et al. , 2005Mateus et al. 2005). STARLIGHT fits an observed continuum spectral energy distribution using a combination of multiple simple stellar    Figure 2. Upper panel, in black, optical spectrum of the galaxy J1624 in the spectral range 3500 -3912Å, around the Balmer jump and the high order Balmer series, together with the fitting (red line) performed using the STARLIGHT spectral synthesis code. Lower panel, subtraction of this fitting from the WHT spectrum in the same spectral range; the solid red line shows the zero level.
population synthetic spectra (SSPs; also known as instantaneous burst) using a χ 2 minimization procedure. In order to keep the consistency with the photoionization models used to model the gas emission, we have chosen for our analysis the SSP spectra from the Starburst99 libraries (Leitherer et al., 1999;Vázquez & Leitherer, 2005), based on stellar model atmospheres from Smith et al. (2002), Geneva evolutionary tracks with high stellar mass loss (Meynet et al., 2004), a Kroupa Initial Mass Function (IMF; Kroupa, 2002) in two intervals (0.1-0.5 and 0.5-100 M⊙) with different exponents (1.3 and 2.3 respectively), the theoretical wind model (Leitherer et al., 1992) and a supernova cut-off of 8 M⊙. We have fixed the metallicity of the stellar populations to Z = 0.004 (= 1/5 Z⊙) and Z = 0.008 (= 2/5 Z⊙) depending on the closest total oxygen abundance as measured using the direct method for each galaxy ( see H06 and H08). The STARLIGHT code solves simultaneously the ages and relative contributions of the different SSPs and the average reddening. The reddening law from Cardelli, Clayton & Matthews (1989) has been used. Prior to the fitting procedure, the spectra were shifted to the rest-frame, and re-sampled to a wavelength interval of 1Å in the entire wavelength range between 3500Å and 9000Å by interpolation conserving flux, as required by the program. Bad pixels and emission lines were excluded from the final fits.
To illustrate the fitting results, we show in the upper panel of Figure 1 the spectrum of J0021 with the best fit found by STARLIGHT overimposed and, in the lower panel of the same figure, the emission line spectrum once the underlying stellar continuum has been subtracted. The measured intensities of the most prominent Balmer emission lines in this subtracted spectrum are, in all cases, equal within the errors to those reported in H06 and H08, where a pseudo-continuum was adopted. It was also shown in those works that the measurement of the intensities of the Balmer lines is consistent with a multi-gaussian fitting of the lines. This method fits independently the emission line and the absorption feature, which is visible through their wings in some bright Balmer lines. In the upper panel of Figure   we plot in detail the WHT spectrum of J1624 in the spectral range 3500-3912Å, around the Balmer jump and the high order Balmer series, and in solid red line we show the STARLIGHT fitting. In the lower panel of the same figure we show the subtraction of the fit from the observed spectrum. We can appreciate that small errors in the fit of the underlying stellar population could become a great unquantifiable error in the emission line fluxes. Nevertheless, for the strongest emission lines, the differences between the measurements done after the subtraction of the STARLIGHT fit and those using the pseudo-continuum approximation are well below the observational errors, giving almost the same result. It must be highlighted that in this type of objects the absorption lines, used to fit the underlying stellar populations, are mostly affected by the presence of emission lines.
There are a few that are not, such as some calcium and magnesium absorption lines. To make the STARLIGHT fit it is necessary to mask the contaminated lines, thus the fit is based on the continuum shape and takes into account only a few absorption lines. Hence, it is not surprising that the STARLIGHT results are not so good for this kind of objects. It must be noted too that the STARLIGHT methodology is very advantageous for statistical and comparative studies. Those studies only consider the results derived from the strongest emission lines which have small proportional errors. However, we must be careful when the aim of our work is the detailed and precise study of a reduced sample of objects because, as we have shown, we can not quantify the error introduced by the fitting. Moreover, these errors are higher for the weakest lines, including the auroral ones. Regarding helium absorption features the multigaussian approximation is not possible in any line, because their wings are not visible. Nevertheless, we have checked that the contribution of the depressed components of the emission lines predicted by the STARLIGHT fit are not larger than taking into account a slightly lower position of the continuum, as it is explained in both H06 and H08.
¿From this analysis we only take as a valuable input information for our tailor-made photoionization models, the ratios predicted by STARLIGHT for the ionizing and the underlying stellar populations in order to correct the equivalent width of Hβ and to take it as an estimate of the properties of the ionizing stellar population.
Nevertheless, we can do a comparative analysis of the star formation histories of every galaxy as predicted by the population synthesis fitting technique. In Figure 3, we show the histogram of the age distribution of the contribution to the visual light of each stellar population used Table 2. Properties of the stellar populations as obtained with STARLIGHT for each studied object. In the last two columns we list the measured EW(Hβ) and its value corrected for contribution to the continuum from the young stellar component only. by STARLIGHT to fit the optical SED for the ten studied galaxies As we can see, the blue light in these objects seems to be largely dominated by a very young stellar population (younger than 10 Myr) which is responsible for the ionization of the gas. At same time, there are different underlying stellar populations, with a component older than 1 Gyr in all cases. This "old" component dominates the total stellar mass in all the objects, providing more than 99 % of it, with the exception of J1509, for which this percentage is about 92 %. At any rate, we find very similar age distributions for the stellar populations in all the objects pointing to very similar star formation histories with the presence of recursive starburst episodes resembling the one observed at present. In table 2, we summarize the properties of the best fitting to the optical spectrum of each galaxy as found by STARLIGHT. We list in columns 2 to 5 the visual extinction, A(V), the total stellar mass, M * , the mass of the stellar population younger than 10 Myr, Mion * , and the percentage contribution of this to the total stellar mass, % Mion * . We also list in columns 6 and 7 the measured equivalent width of Hβ, EW(Hβ), and the corrected value, once the contribution to the continuum by the underlying stellar population has been subtracted. As we can see, our procedure yields a negative value for the extinction in J1616. Nevertheless, this is still consistent with the very low extinction found using the Balmer decrements.
We can check the accuracy of these fittings by comparing the stellar mass of the ionizing population as derived from STARLIGHT, and the values of the same masses as derived from the extinction-corrected Hα fluxes, using the expression from Díaz (1998): This comparison is shown in Figure 4. It can be seen that the agreement between both derived values is very good for most of the studied objects. Only J0021 shows a higher young stellar mass when derived from its Hα flux. Horizontal error bars in the plot correspond to the underestimate of the Hα flux due to the UV dust absorption. We have quantified this effect taking the absorption factors derived from the photoionization models described in section 3.3 below. Nevertheless, taking into account these corrected values of L(Hα) in the derivation of Mion * from (1) does not increase the ionizing stellar masses by more than 0.3 dex, except in the case of J1455, with 0.42 dex. . Comparison between the ionizing stellar mass as derived from the extinction-corrected Hα fluxes and from SED fitting for an age younger than 10 million years. The red solid line represents the 1:1 relation. Error bars to the right take into account the dust absorption factors obtained from the photoionization models described below.

Detection and measurement of Wolf-Rayet bumps
Wolf-Rayet (WR) stars appear in the first stages after the main sequence phase of the evolution of massive stars, so they are already visible soon after the beggining of a burst of star formation (approximately 2 Myr) and for a relatively short period of time (about 3 Myr, on average). The existence of this WR stellar population can be very useful for the study and characterization of the ionizing stellar population in starburst galaxies (e.g. Pérez-Montero & . The strength of the stellar winds produced by these stars can be sometimes measured in the integrated spectra of the starburst galaxies, which are then identified as WR galaxies (Conti, 1991). The presence of the emission broad features ("bumps") produced by WR stars are therefore associated to the existence of processes of intense star formation. These WR features are the blue bump, centered Table 3. Observed extinction corrected Hβ intensities and properties of the WR bumps for the objects of the sample. We list equivalent widths and relative intensities, normalized to 100 times the Hβ value, measured at 4650Å (blue bump, bb) and at 5808Å (red bump, rb).
The corrected values take into account the contribution of the underlying stellar population and the nebular emission to the continuum and the dust absorption of the UV predicted by the photoionization models in each object.
We have detected the blue bump in seven of the ten galaxies of our sample, while the red bump has been detected in only three of these seven galaxies. All the three galaxies observed with WHT had already been identified as WR galaxies in H06 and have also been included in the catalogue of WR galaxies of the SDSS (Brinchmann et al., 2008). The other four galaxies with a detection of WR features, observed with the 3.5 CAHA telescope are not listed in that catalogue. An example of the WR emission can be seen in Figure 6 of H06 for the J1624 galaxy.
We have measured the bumps in the following way. First we have adopted a continuum shape and we have measured the broad emission over this continuum at the wavelengths of the blue bump, at 4650Å, and the red bump, at 5808 A. Finally, we have subtracted the emission of the narrower lines not emitted by the WR stars winds, which are [Feiii] 4658Å, the narrow component of Heii 4686Å, Hei + [Ariv] 4711Å and [Ariv] 4740Å in the blue bump. The equivalent widths and the deredenned relative intensities of the blue and red bumps are shown in Table 3. Both have been corrected in order to study the properties of the ionizing stellar population. In the case of EWs, we have subtracted the contribution of the underlying stellar population as derived from STARLIGHT and from the nebular continuum as derived from the photoionization models described below. In the case of the relative intensities, we have taken into account the dust absorbed component of Hβ derived from the same photoionization models for each object.

Photoionization models of the ionized gas
Photoionization models have been calculated using the code CLOUDY v. 06.02c (Ferland et al., 1998) with the same atomic coefficients used to derive ionic abundances in H06 and H08 and the recent dielectronic recombination rate coefficients from Badnell (2006). We have taken the same Star-burst99 stellar libraries used by STARLIGHT to fit the observed SEDs. We have assumed different initial input param-eters and then, we have followed an iterative method to fit the relative intensities of the strongest emission lines in the integrated spectra, along with other observable properties.
We have assumed a constant star formation history, with a star formation rate calculated according to the measured and corrected number of ionizing photons. The age of this burst of star formation has been limited to the range 1-10 Myr in order to reproduce the corrected equivalent width of Hβ once the older stellar populations have been removed, except in the case of J1540 for which a burst duration of 15.85 Myr has been considered in order to reproduce the observed EW(Hβ). Total abundances have been set to match the values derived in H06 and H08, including He, O, N, S, Ne, Ar and Fe. The abundances of the rest of the elements have been scaled to the measured oxygen abundance of each galaxy and taking as reference the photosphere solar values given by Asplund et al. (2005). The inner radii of the ionized regions have been chosen to reproduce a thick shell in all cases. This kind of geometry reproduces best the ionization structure of oxygen and sulphur simultaneously . We have assumed a constant density equal to that derived from the measured [Sii] emission line ratio (∼ 100 particles per cm 3 in all the objects). We have kept as free parameters the filling factor and the amount of dust. The inclusion of dust allows to fit correctly the measured electron temperatures because it is a fundamental component to reproduce the thermal balance in the gas. The heating of the dust can affect the equilibrium between the cooling and heating of the gas and, thus, the electron temperature inside the nebula. We have assumed the default grain properties given by Cloudy 06, which has essentially the properties of the interstellar medium and follows a MRN (Mathis, Rumpl & Norsieck 1997) grains size distribution.
After some first estimative guesses of the UV absorption factor due to the internal extinction, f d , we have recalculated the age of the cluster and the number of ionizing photons to match the observed values. Finally, we have adjusted the chemical abundances to fit the relative intensities of the emission lines.
The average number of models needed to fit the observable properties of the galaxies with a reasonable accuracy is about 50 for each galaxy. In Table 4 we show the intensities of the most representative emission lines as measured in H06 and H08 compared to the values obtained with the best model for each object. The agreement in most cases is excellent, even for the faint auroral lines.
In Table 5, we show the different line temperatures obtained from the photoionization models compared to those derived in H06 and H08. In the same table, we show the elemental ionic abundances derived using the corresponding electron temperature representative of the region where the ion is formed: ) in the other cases. These abundances are compared with the ionic abundances obtained from the corresponding models. We also show the temperature from the oxygen recombination lines predicted by our models in the high excitation region and the temperature fluctuations (t 2 , Peimbert, 1967) obtained from the difference between the temperatures derived from recombination and collisional emission lines. It can be seen that the predicted fluctuations are practically zero and, therefore, no ADFs are expected according to our models. The comparison between the model predictions and the temperature fluctuations inferred from the observed values of t([Oiii]) and the Balmer jump temperature for the three objects observed with the WHT gives compatible results, with the exception of J0032, whose t 2 is higher as obtained from the observations.

Properties of the ionizing population
In Table 6 we list the main observed properties of the young stellar population responsible for the ionization of the surrounding gas and their comparison with the same values as predicted by the corresponding tailor-made models. Column 2 lists the age of the burst of continuous star formation adopted for each model in Myr. Column 3 gives the metallicity (Z) as derived from the relative oxygen abundance measured in the gas-phase and column 4 the adopted Z of the stellar ionizing cluster, which is in each case the value closest to the oxygen abundance. Although the hypothesis of equal metallicity for the gas and the ionizing stars is not well justified and the available metallicities in the stellar libraries are still poor to some extent, the differences in metallicity among models are not crucial to improve their accuracy. Column 5 in the table lists the number of ionizing photons, Q obs (H), as derived from the observed luminosity of Hα using the following expression (see, for example, Osterbrock, 1989): with L(Hα) in erg · s −1 . This Q abs (H) should be compared to the the number of photons produced by the model cluster once divided by the dust absorption factor, f d , predicted by the models, which is listed in column 6. The total number of ionizing photons produced by the ionizing cluster is then: where both Q(H) are expressed in photons · s −1 . Therefore, in order to be compared with the observations we list in column 7 the ratio Q(H)/f d predicted by the models. Finally, columns 7 and 8 give the equivalent width of Hβ, EW(Hβ), measured on the spectra once the underlying stellar populations have been removed. This is to be compared with the corresponding EW(Hβ) value given by the models, which already takes into account the dust absorption factor and the nebular continuum emission. In the lower panel of the table we show in columns 2 and 3 the logarithm of the ionization parameter, log U , calculated using the following expression derived by Díaz et al. (1991)  which offers a reliable estimation of the ionization parameter in nebulae with a non-plane parallel geometry. This value is compared with the corresponding log U obtained from our models. The logarithm of the filling factor (log ǫ) and the adopted inner radius (R0) of the ionized region are listed in columns 4 and 5 while columns 6 and 7 give the external radius (RS) from the models as compared with the radius derived from the spatial profile of Hα in the long slit. Finally, the last three columns give the dust-to-gas ratio in each model and the visual extinction, A(V), as derived from the observed Balmer decrements and compared with the values obtained from the models. As in the case of the emission lines, the agreement between the derived properties and those obtained from the models is good with the exception of the outer radii which are much larger in the observations by factors ranging from 2.2 (J1616) to 15.5 (J1540). Possible explanations for these discrepancies could be either a very low density layer of diffuse gas in the outer regions of the ionized gas or lower filling factors as compared with the values found in our models. We find as well some discrepancies between the extinction values found in our models and those derived using Balmer decrements. This could be symptomatic of very complex structures of dust in the ionized gas region, not easy to reproduce in the models.
Regarding the age of the burst, the equivalent width of an intense Balmer emission line, like Hβ, is generally associated to the age and properties of the ionizing population. It has long been accepted that the low values of EW(Hβ) found in starburst galaxies as compared with the predictions from Population Synthesis Models are due to the presence of an underlying stellar population that enhances the continuum light without significantly contributing to the ionization of the gas. Nevertheless, as can be seen in Table 2 the subtraction of this population from the integrated spectra does not imply significant corrections to EW(Hβ) (between 1.5 % in the case of J1540 and J1729 and 24 % in the case of J0032).
However, the excitation conditions found in these objects, are only reachable with the energy supplied by very young and massive O stars, whose associated EW(Hβ) is four or five times larger than the observed values. This difference can be partially explained by a certain dust absorption of UV photons. The dust absorption factors obtained in our models range between 1.5 and 2, as we can see in Table 6. These help to put the modeled EW(Hβ) values in agreement with the observed ones.  . Left panels: relation between the equivalent width of Hβ and the age of the corresponding synthesis population model for metallicities of 0.004 and 0.008 as predicted by Starburst99. Crosses represent the observed values measured on the young stellar continua and multiplied by the corresponding absorption factor obtained in each tailor-made model. The assigned age to each object is that predicted by the same models. Therefore, the J1540 point does not appear as the fitted age is older than 10 Myr. Right panels: relation between the ratio of helium to hydrogen ionizing photons and the age of the ionizing cluster. In all panels, the black solid line represents an instantaneous star formation and the red dashed line, a continuous star formation as predicted by Starburst99. As in the previous case, crosses represent the values obtained from the models for each object, with the exception of J1540.
In the left panels of Figure 5 we show the evolution of EW(Hβ) with time for the two metallicities closest to those in our sample for both instantaneous (black solid line) and continuous (dashed red line) star formation laws as obtained from Starburst99 models. In the right panels of the same figure we show the evolution of the ratio of helium to hydrogen ionizing photons, Q(He)/Q(H) according to the same Star-burst99 models. This ratio, that can be considered a good estimate of the ionizing power of the cluster, falls dramatically for an instantaneous star formation law once the WR phase has finished (∼ 5 Myr), while it maintains a constant high value in the case of a continuous star formation. We Table 5. Comparison between the derived electron temperatures and ionic abundances derived in H06 and H08, and the values obtained from the best tailor-made photoionization model for each object.  Table 6. Comparison between some properties obtained from the best tailor-made model and derived from the observations for each object, including the age of the ionizing cluster, metallicity of the gas and the adopted metallicity of the cluster, dust absorption factor, number of ionizing photons, equivalent width of Hβ, ionization parameter, filling factor, inner and outer radii, dust-to-gas ratio and visual extinction. can compare these diagrams with the observed EW(Hβ), once the underlying stellar population has been removed and the dust absorption factor has been taken into account. They appear as crosses in Figure 5 and lie in a range which corresponds to the asymptotic value of the continuous star formation line. In the case of the evolution of Q(He)/Q(H), crosses represent the values derived from the tailor-made models.

ID
Although it is true that both EW(Hβ) and Q(He)/Q(H) could be consistent with an instantaneous burst at present in the WR phase, statistically it is not probable that all the objects of the sample present such similar properties due to their being in the same evolutionary stage. At the same time, all the calculated photoionization models in our work require a great amount of O and B stars younger than the WR phase in order to reach the ionization conditions responsible for the observed emission lines. We have therefore assumed for all the studied objects a continuous star formation history for the ionizing source, with a duration longer in all cases than 5 Myr.

Derived properties of the WR population
Seven out of the ten observed galaxies have detectable WR features in their spectra. Five of them were already included in the catalog by Kniazev et al. (2004;SHOC) with a flag indicating the detection of both the blue and red bumps from WR stars. In our spectra, the red bump is observed in only three galaxies. The SHOC catalog, however, does not provide actual measurements for the features. In a more recent paper, Brinchmann et al. (2008) list only the "blue bump" luminosities and equivalent widths of three galaxies from the present work: J002101 (SHOC 11), J003218 (SHOC 22) and J162410 (SHOC 536). The direct comparison between the equivalent widths of the blue bump measured by these authors inside the 3-arcsec fiber spectra of the SDSS and our long-slit observations yields discrepancy factors for J0021, J0032 and J1624 of 1.7, 2.0 and 1.1, respectively, larger in our WHT observations. This could be a consequence of the smaller aperture of our long-slit observations which includes a lower contribution to the continua from extended ionized gas and hence produce larger equivalent widths. Table 7 shows the extinction corrected luminosities of the measured WR bumps. We show in the same table the number of O and WR stars as derived from the extinction corrected Hα and WR blue bump luminosities, respectively. In the case of O stars, we have taken into account as well the dust absorption factor, f d . The total number of stars have been estimated by comparison with the properties of the Starburst99 models at the age derived from the photoionization models described above. For those galaxies with a red bump detection we have used the same libraries to derive the ratio between WC and WN stars, which is also listed in Table 7.
We can compare the relative intensities and equivalent widths of the WR bumps detected in our objects with the predictions of Starburst99 populations synthesis models as a function of the metallicity, age and star formation law of the cluster. To make this comparison we have taken the corrected values listed in Table 3: in the case of EWs from the effect of the underlying stellar populations and the contribution of the nebular continuum and in the case of the relative intensities from the effect of the dust absorption factor to the hydrogen recombination lines. In Figure 6 we show the predicted equivalent widths and intensities relative to Hβ of both the blue and red bumps as a function of the age of the cluster for metallicities Z = 0.004 (0.2·Z⊙), left panels, and Z = 0.008 (0.4·Z⊙), right panels. The black solid lines represent the evolution of an instantaneous burst, while the red dashed lines do for a continuous star formation history with a constant star formation rate. Crosses represent the corrected values as listed in Table 3 with their corresponding error bars. As we can see, the models predict the brightest and most prominent features for the highest metallicity, agreeing with the stellar model atmospheres for WR stars (Crowther, 2007 and references therein). At the same time, we can appreciate that in the instantaneous star formation case the WR features appear during a time interval between 2 and 5 Myr reaching higher intensities than in the case of continuous star formation. In this latter scenario, the WR features appear at the same age and, despite reaching lower intensities, they converge to a non-zero value at older ages. For an average metallicity between the two cases of study these values are EW ∼ 2Å and 0.013·I(Hβ) for the blue bump and EW ∼ 1Å and 0.004·I(Hβ) for the red one. By comparing with the measured relative intensities and equivalent widths in our sample objects we can see that all of them are larger than the Starburst99 model predictions for a constant star formation law and, in some of them, like J1509 and J1540, even larger than the values predicted for an instantaneous burst.
A great deal of work has been done in order to reconcile this disagreement between the observed luminosities of the WR bumps detected in the integrated spectra of galaxies and the model stellar atmospheres of WR stars and evolutionary synthesis models of clusters. Suggested hypotheses include the existence of a binary channel (Eldridge et al., 2007) or the consideration of rotation in the stellar evolutionary tracks, which reduces the mass required to form WR stars (Meynet & Maeder, 2005). Nevertheless, some other reasons can be explored such as geometrical effects. If we consider the whole burst of star formation and we assume a continuous star formation law, we must compare the luminosities of the WR bumps with the intensities of the Balmer lines and the continuum coming from the whole burst. This approach was already used by  to match the observed fluxes of the WR features in Mrk 209. On the contrary, if only the ionizing stellar population producing the WR stars is considered, only the burst region where the WR population concentrates must be analyzed. This could explain as well the lack of WR emission in the other three galaxies with an assumed continuous star formation history. In the case of our sample, no aperture correction factors have been taken into account either for the Hα luminosities or the continuum contribution at the WR bumps wavelength, due to the compact nature of the objects. Indeed, the discrepancy factors between our Hα luminosities and those measured in the SDSS catalogue using a 3 arcsec fiber is not larger than 1.5 in any case and even smaller than 1 in two of them (J1528 and J1729). Therefore, the aperture factor does not seem to be the cause of the abnormally high values found for the EW and the relative intensities of the WR features.
The analysis of the relative number of WR to O stars for the ratio of WC to WN stars. In both cases, results for two metallicities, Z= 0.004 and Z=0.008, are shown. Solid lines represent instantaneous star formation models and dashed line continuous star formation ones. Crosses represent the derived numbers for each galaxy taking into account the measured luminosities of the bumps as compared with those predicted by the Starburst99 models at the predicted age. The number of O stars has been derived using the same procedure for L(Hα). Table 7. Luminosities of the Wolf-Rayet features measured for the sample of HII galaxies at 4650Å (blue bump, bb) and at 5808Å (red bump, rb), derived number of WR stars, ratios of WR to O stars and, for those galaxies with a detection of the red bump, of WC to WN stars.
Object ID log L(bb) log L(rb)  . Comparison between the electron temperatures (in units of 10 4 K) as derived from the observations and from the tailormade photoionization models. Black circles represent a comparison between a temperature directly derived in the observations, while white squares represent a temperature of another ion, assumed to be similar in the thermal structure. and of the number of WC to WN stars leads to very similar conclusions. In all cases the derived ratio of the number of WR to O stars is much higher than the asymptotic value reached with a continuous star formation, which is 0.022 for a metallicity in between the two studied cases. For J1540, the obtained value, 0.44, is even much higher than the value expected for an instantaneous burst. As in the case of the luminosities of the WR features, we can reconcile the derived values enhancing the number of O stars on the basis that a non-negligible flux of the radiation from these stars is obscured by internal extinction. As in the previous case, probably geometrical effects can be determinant. Regarding the ratio between WC and WN stars, the derived values are lower in all cases than expected. This is due to the fact that the luminosities of the red bump predicted by the models are in much better agreement with the observed values than those of the blue bump and, hence, the relative number of WC stars is lower.

Physical conditions and ionic abundances of the gas through tailor-made models
The correct determination of ionic abundances in ionized gaseous nebulae using collisional emission lines versus recombination lines has become a matter of debate in the last years due to the disagreements predicted by photoionization models in the high metallicity regime (e.g. Charlot & Longhetti, 2001) or the discrepancies found in local nebulae (e.g. Peimbert et al., 2007 and references therein). Many times, these discrepancies are due to the simple assumptions made about the ionization structure of the nebulae. Indeed, the only electron temperature measured in many cases is t([Oiii]) because the emission line of [Oiii] at 4363Å is the brightest auroral line. This temperature is assumed to be representative of the high excitation zone and the low excitation zone temperature is obtained under very simple assumptions. This can produce large deviations in high metallicity nebulae, where the low excitation ions have a high relative weight in the total abundance of a specific species (Pérez-Montero & Díaz, 2005). One way to overcome this issue is to try to trace as precisely as possible the inner thermal structure of the gas.
In Figure 8 we show a comparison between derived and modeled temperatures for the objects of our sample, for which at least four electron temperatures have been measured. In the case of ions whose temperatures can not be derived due to the lack of observable lines in the optical, another line temperature, assumed to be representative of the zone where these iones are formed, has been used: t ( [Nii] in the objects where this was not available. As we can see, the agreement for t([Oiii]) is excellent, mostly due to the input conditions in the models to fit the auroral line of [Oiii]. Regarding t([Oii]) and t([Siii]), the agreement is still good, although slightly worse that in the previous case. Finally, the largest deviations are found for the electron temperature of [Sii]. Except in the cases of J1455 and J1616, models predict much higher temperatures than observed. This is not surprising if we take into account that the emission of [Sii] is often found as well in the diffuse medium (e.g. Hunter & Gallagher, 1990) and therefore we are probably tracing in our observations a non-negligible portion of colder and less dense gas, which has not been considered in our models. This is compatible with the significatively lower radii found in our models in comparison with the values obtained from the spatial profile of Hα. Perhaps, this could be overcome if we considered in our models a lowdensity region in the outer parts of the nebula where lines from low excitation ions like [Sii] were produced. In the case of the three objects with a measurement of the electron temperature of [Nii], the agreement is rather good in the case of J0021, while it is sensibly underestimated by the models in    Figure 10. Comparison between ionic abundances as derived from the observations using the direct method and the predictions of tailor-made photoionization models.
the case of J1624 and J1729. This temperature is assumed to be equal to t([Oii]) in objects without a measurement of the auroral line of [Nii] at 5755Å. As we can see in the corresponding panel, this assumption seems to work for most of the objects, except in the cases of J0032 and J1455. For these objects, we can not rule out some geometrical effect not considered in our photoionization models. In Figure 9 we illustrate the good agreement found between models and observations for the galaxy J1616. In each panel, in red, green and blue solid lines, the radial distribution of the ionic fractions of the different elements as computed by the model (from left to right and from up to down: O, S, N, Ne, Ar and Fe) are plotted. The magenta thick solid lines show the radial distribution of the electron temperature predicted by the same model. Finally, black squares show the electron temperatures obtained from the observations. We have associated these points to the average zone where the corresponding ion is formed. In the cases of Ne 2+ and Fe 2+ , t([Oiii]) is shown; for Ar 2+ ,t([Siii]) is plotted and in the cases in which no measurements of t([Nii]) are available, this is represented by t([Oii]).
In Figure 10, we show the comparison between the relative ionic abundances, derived using the corresponding measured electron temperatures and the observed emission-line intensities, and the same abundances as predicted by our photoionization models. We can see that the agreement is somewhat better for high excitation ions, like O 2+ , than for low excitation ones, like O + and S + . At the same time these latter ions present abundances spanning a wider range than the former. This might indicate that geometrical effects (e.g. density variations, matter bounded nebulae, etc...) affect the diagnostics of the gas in the low excitation region, hence a simple scheme can not be adopted for its analysis. In the case of our sample of HII galaxies, these differences are not appreciable due to their low metallicities and high excitation conditions, so the inner parts have a large weight on their average conditions. On the contrary, the analysis of disk HII regions and high metallicity objects are tied to a much larger uncertainty if these effects are not taken into account using as many indicators of the physical conditions of the gas as possible.

The hardening of the radiation
The stellar effective temperature has been traditionally the functional parameter most difficult to obtain. Vílchez & Pagel (1988) defined the so-called softness parameter, η, which allows to estimate the hardening of the radiation field by means of the ratio of the ionic abundances of oxygen and sulphur, derived from the optical spectrum.
The same ratio, defined using the intensities of the corresponding emission lines, can also be used to figure out the scale of equivalent effective temperature: Both expressions agree in finding that HII galaxies present very high equivalent effective temperatures as compared to other HII regions of higher metallicity. Nevertheless, there is a disagreement between the position of these objects in diagnostic diagrams and the scale of temperatures of the stellar atmospheres of O stars. HII galaxies usually lie in a region of the diagrams which corresponds to temperatures even higher than the maximum value of model atmospheres, which is 50 kK. The cause has been sometimes attributed to some extra source of heating in HII galaxies (e.g. Stasińska & Schaerer, 1999). The left panel of Figure 11 shows the logarithmic ionic fractions O + /O 2+ vs. S + /S 2+ . The different lines of constant slope 1 in this diagram correspond to different values of the η parameter and hence of ionizing radiation temperature. The values derived from actual measurements of our observed objects are plotted as black solid circles while blue squares correspond to the values derived from the computed models. Red solid lines represent polynomial fits to different photoionization models computed using as ionizing source WM-Basic stellar atmospheres (Pauldrach et al., 2001) for different values of the effective temperature. These model sequences can help to establish the scale of the equivalent effective temperature. We can see that the models are not able to reproduce the high values of log(S + /S 2+ ), and in some cases log(O + /O 2+ ) found in our objects that, in fact, would correspond to extremely high values of equivalent effective temperatures.
This disagreement disappears when the emission lines instead of the actual ionic ratios are considered. This is due to the fact that emission-line intensities are well fitted in all the tailor-made models while this is not the case for electron temperatures (see Figure 8), specially for T([SII]) which is, in general overestimated by the models leading to an underestimate of the S + /H + abundance ratio. Again, the disagreement could be caused by the geometry of the diffuse gas in these HII regions, which is not adequately modeled, thus avoiding the need for an extra source of heating to explain the apparently high stellar effective temperatures in these objects. At the same time, this probably means that the calibration of the effective temperature is not unique but in fact sensitive to the object class involved, through the different geometries of the ionized gas.

SUMMARY AND CONCLUSIONS
We have studied the underlying and ionizing stellar populations and have modeled the properties of the emitting gas in a sample of HII galaxies. This sample was observed (H06 and H08) using double-arm spectrographs, which allow the simultaneous analysis of the entire spectral range from [OII] 3727Å to [SIII] 9532Å providing the derivation of the electronic temperatures of, at least, four ions in all the cases: O + , O 2+ , S + , and S 2+ and hence a precise determination of the ionization structure of the gas.
The fitting of the blue part of the continuum was made using the program STARLIGHT and the SEDs computed from the Starburst99 stellar population synthesis models. This approach was used to study in a consistent manner the differences between the stellar populations in our sample of galaxies. The resulting distribution of ages and masses from this fitting revealed the presence of various stellar populations. An underlying stellar population of about several Gyr dominates the stellar mass, while the blue to visual light is in all cases dominated by young stellar populations with ages younger than 10 Myr. The presence of the older population can affect the reddening determination but its contribution to the measured EW(Hβ) is only about 10% for most of the sample objects, being around 25% in the case of J0032 and almost zero for J1616 and J1729. In fact, regarding the subtraction of Balmer absorption features from the underlying stellar populations, we did not find that the spectral fitting technique constitute a substantial improvement over the simultaneous fitting of absorption and emission components as performed in previous works (H06 and H08).
We analysed the Wolf-Rayet clusters detected in seven of the ten studied galaxies through the WR broad features seen in their spectra. The fitting of the stellar population and the modeling of the ionized gas allows to correct them for the contribution of the underlying population and the nebular emission to the continuum and for the dust absorption factor in the case of the emission of the Balmer lines. The comparison between the corrected fluxes and equivalent widths of the WR blue bump with the predictions from Starburst99 models at the corresponding metallicity leads to the estimate of the number of WR stars which ranges from 9±6 to 2500±600 depending on the galaxy. In three of them, also the red bump was measured allowing to estimate the WC to WN number ratio, Nevertheless, the fluxes of both the blue and red bumps relative to Hβ are larger than the values predicted by the models in almost all the cases.
A continuous star formation history of the stellar ionizing population has been assumed in all the objects of the sample. This is supported by the simultaneous presence of very massive O and WR stars, along with the relatively low measured values of EW(Hβ), once corrected for underlying contribution and dust absorption. The ages of the star formation bursts, the number of ionizing photons as derived from the Hα fluxes and the intensities of the most prominent emission lines were used to constrain photoionization models for each galaxy. The Starburst99 libraries were used to calculate the SED of the ionizing cluster in order to be consistent with the previous fitting of the stellar populations. Some dust is assumed to be present associated with the gas. The amount of dust in each model was derived in order to reproduce the observed EW(Hβ). These derived values agree well with the main observed quantities in each object, with the exception of the sizes of the ionized regions which result 5-10 times smaller in the models as compared to measurements from Hα in the spatial profiles of the slits.
The observed electron temperatures and ionic abundances are, in most cases, well reproduced by the models which implies that the direct method using collisional emission lines provides robust answers if the thermal structure of the nebula is well traced. Therefore, our models were able to reproduce the ionization structure of the nebula without appealing to temperature fluctuations or any other geometrical effects that can affect the abundances derived from collisionally excited emission lines. This result is consistent with the temperature measured from the Balmer jump in some of the objects in our sample (H06). Only the electron temperature and the ionic abundances of S + /H + is not well reproduced by our models, which could be consequence of the presence of diffuse gas in these galaxies, not taken into account in the models. This would also be consistent with the small sizes predicted by the models as compared with Hα spatial profiles in the slits.
Finally, regarding the hardness of the radiation field, model predictions agree with observations when the softness parameter η, which parametrizes the temperature of the ionizing radiation, is expressed in terms of emission line intensity ratios. However, there is a disagreement when ionic abundance ratios are used instead. This is probably caused by the overestimate of the electron temperature of S + by the models and the corresponding underestimate of the S + /S 2+ abundance ratio rather than by existence of an additional heating source in HII galaxies. Table A1. Ionization correction factors and total abundances as derived from photoionization models of the main elements with emission lines in the optical spectral range of the studied HII galaxies. ICF(S + +S 2+ ) 1.11 1.14 1.28  Figure A1. Comparison between different ionization correction factors (ICFs) derived for different elements from the computed photoionization models (black circles) and some of the most used relations.

T([Ariii]) ≈ T([Siii]
) (Garnett, 1992 In the same work, another ICF is proposed for the case in which [Ariv] is also measured: These two ICFs can be expressed as well as a function of the ionic abundances of sulphur, S + and S 2+ , when only red spectroscopic observations are available: ICF (Ar 2+ ) = 0.596 + 0.967 · 1 − S 2+ S + + S 2+ + + 0.077 · 1 − S 2+ S + + S 2+ −1 (A9) In the bottom panels of Figure A1, we can see that the agreement between these expressions and the ICFs found in our models is good.
Finally, iron abundances in the gas-phase can be calculated by correcting the ionic abundances of Fe 2+ , derived from the intensity of the emission line of [Feiii] at 4658Å, using the following expression from Rodríguez & Rubin (2004, RR04 in the plot): with a value of α = 0.09. Nevertheless, the fit to the ICFs found in our models gives a best value for α = 0.24, which leads to lower total abundances of iron.