A&A 447, 709-717 (2006)
DOI: 10.1051/0004-6361:20052979

Spectral energy distribution for GJ406

Ya. V. Pavlenko1,2 - H. R. A. Jones1 - Yu. Lyubchik2 - J. Tennyson3 - D. J. Pinfield1

1 - Centre for Astrophysics Research, University of Hertfordshire, College Lane, Hatfield, AL10 9AB, UK
2 - Main Astronomical Observatory, Academy of Sciences of the Ukraine, Golosiiv Woods, Kyiv-127, 03680 Ukraine
3 - Department of Physics and Astronomy, University College London, Gower Street, London WC1E 6BT, UK

Received 3 March 2005 / Accepted 18 October 2005

We present results of modelling the bulk of the spectral energy distribution (0.35-5 $\mu$m) for GJ406 (M6V). Synthetic spectra were calculated using the NextGen, Dusty and Cond model atmospheres and incorporate line lists for H2O, TiO, CrH, FeH, CO, MgH molecules as well as the VALD line list of atomic lines. A comparison of synthetic and observed spectra gives $T_{\rm eff}$ =  $2800\pm 100$ K. We determine $M_{\rm bol} = 12.13 \pm 0.10$for which evolutionary models by Baraffe et al. (2003, A&A, 402, 711) suggest an age of around 0.1-0.35 Gyr consistent with its high activity. The age and luminosity of GJ406 correspond to a wide range of plausible masses (0.07-0.1  $M_{\odot}$).

Key words: stars: individual: GJ406 - molecular data - stars: fundamental parameters - stars: late-type - stars: atmospheres - stars: evolution

1 Introduction

Studies of M dwarf spectra are of interest to many branches of modern astrophysics. Indeed, perhaps 70% of stars within 10 parsecs are M dwarfs and it is very probable that this number density prevails throughout our Galaxy. The population of these numerous low-mass stars (0.08  $M_{\odot}$ $\la$ M $\la$ 0.6  $M_{\odot}$), together with substellar objects (brown dwarfs; $M \leq 0.075$  $M_{\odot}$) would contain an appreciable amount of the baryonic matter in the Galaxy. Estimates of brown dwarf number densities currently suggest the same order as for stars ($\sim$0.1 per pc-3), therefore their contribution to the total mass should not exceed 15% (Reid et al. 1999). Nonetheless, the large errors associated with age and mass determinations for brown dwarfs make such estimates very uncertain.

The verification of the theory of stellar evolution and structure of stars, the detection among M dwarfs of a subset of young brown dwarfs, and the physical state of plasma in their low temperature atmospheres are among a few of the interesting problems that may be addressed through the detailed study of M-dwarfs.

Table 1: Telescope and Instrument configurations used to collect our GJ 406 dataset.

Some authors reference GJ406 (other names are V*CN Leo, EUVE J1056+07.0, [GKL99] 228, GSC 00261-00377, LFT 750, LHS 36, 2MASS J10562886+0700527, 2RE J1056+070, 1RXS J105630.3+070118) as an "archetype dwarf of spectral type M6V'', or one of the "well known spectral standards for its type'' (Mohanty et al. 2004). GJ 406 is located at 2.39 pc from the Sun (Henry et al. 2004). Altena et al. (1995) determined a proper motion of $\mu =4.696$. Radial velocities are of order $19 \pm 0.1$ km s-1 (Martín et al. 1997; Mohanty & Basri 2003; Fuhrmeister et al. 2005). Leggett (1992) found that this nearby dwarf has typical old disk properties. Deflosse et al. (1998) reported a rather low $v*{\rm sin} i~$ < 3 km s-1 (see also Mohanty & Basri 2003). Guetter et al. (2003) even use GJ 406 as one of the JHK standard stars on the CIT photometrical system. Indeed, optical and IR spectra of the dwarf do not contain any unusual features. They are governed by absorption of diatomic molecular band systems, such as TiO and VO (Kirkpatrick et al. 1991), as well as rotational-vibrational bands of H2O and CO (Jones et al. 1994).

On the other hand, GJ 406 has some very unusual properties. SIMBAD (http://simbad.u-strasbg.fr/sim-fid.pl) labels GJ 406 as a flare star. The dwarf has strong H$_{\alpha}$ of EW = 6.7 Å ( $L_{\rm H_{\alpha}}/L_{\rm bol} = -3.9$, see Mohanty & Basri 2003). An ultraviolet spectrum of GJ406 contains some emission lines (Fuhrmeister et al. 2004). Furthermore, Schmitt & Wichmann (2001) detected the Fe XIII forbidden coronal line at 0.33881 $\mu$m. Fuhrmeister et al. (2004) reported a high level of variability of this line on a timescale of hours which they ascribe to microflare heating. Recently an X-ray luminosity of log Lx = 26.97was detected by Schmitt & Liefke (2004). GJ 406 is the only known M 6 star yet observed with a strong chromospheric and coronal activity. Only a few stars in the solar vicinity are known with such a menagerie of activity phenomena.

In this paper we compute the synthetic energy distribution of several model atmospheres with a range of effective temperature and compared them with the observed fluxes of GJ406. Section 2 presents the spectral data used in our paper. Section 3 describes our procedure for computation. We psesent our results in Sect. 4. In Sect. 5 we discuss the implications of our results.

\end{figure} Figure 1: Observational data used for this paper. The wavelength coverage of the different instruments is shown.
Open with DEXTER

2 Observations

Table 1 lists the data and instruments used to obtain the observational spectra in this paper. The spectra are shown in Fig. 1 and come from measurements taken with a variety of different instruments on different telescopes. All the data are assessed to be of good quality, most have already been used for other papers. We refer to the spectra in wavelength order. The reduction procedures for the Integral Spectrometer (ISIS) data taken on the William Herschel Telescope (WHT) are described in Dobbie et al. (2004). The reduction procedures for the CGS4 data taken on the United Kingdom Infrared Telescope (UKIRT) here have been reported elsewhere in Jones et al. (1994, 1995, 1996). It should be noted that CGS4 observations in the spectral regions labelled in Fig. 1 by "m'' and "n'' are not continuous. Fluxes between them were filled by NextGen synthetic spectra. The reduction procedures for the Short Wavelength spectrometer (SWS) data taken with the Infrared Space Observatory (ISO) are reported elsewhere in Jones et al. (2002). The Near Infrared Camera Mosiac Spectrograph (NICMOS) data taken on the Hubble Space Telescope (HST) was reduced using the data processing software Calnicc (Version 2.5.7). The data has been compared with Jones et al. (1994) and is preferred due to its excellent flux calibration (5-10%, Pirzkal & Freudling 1998) and broad wavelength coverage.

In general we have renormalised the fluxes for the different spectral regions so as to ensure that flux levels are the same where the regions overlap. For the CGS4 spectra centered on 4.7 $\mu$m (region n), this was not possible since there is no overlap. For this spectral region, we simply used the instrument calibration to determine the flux level, and filled in the gap with a normalised section of NextGen model spectra.

2.1 Absolute flux calibration

In order to provide an absolute normalisation for the full GJ406 spectrum, we used the available near infrared photometry from Leggett et al. (2000), which we transformed onto the Mauna Kea Observatory (MKO) photometric system using Hawarden et al. (2001). We used the measured MKO J, H and K filter pass bands convolved with atmospheric transmission to estimate the ground based photometrically measured flux components in our spectra, by integrating over each band. We then performed the same task on a flux calibrated spectrum of Vega. Vega was assumed to be zero magnitude at all wavelengths, and the GJ406 flux could thus be scaled to match the photometry. The normalisation value that we found was $\sim$20% higher when normalising in the K-band than in the J-band, and was an intermediate value in the H-band. These differences presumably result primarily from the relative normalisations that we used to join the individual spectral regions, and provide an accuracy gauge on this procedure. We chose to make our final normalisation in the H-band, and duly estimate a likely $\pm$10% uncertainty in our absolute flux levels.

In order to derive a bolometric flux, we added a 4.8-20 $\mu$m spectral tail to our calibrated spectrum (using a NextGen 2800 K, $\log{g}=5.0$, [M/H] = 0 model spectra, shown to be appropriate in Sect. 4.2), normalised in its overlap region, and integrated out to 20 $\mu$m. Note however, that this synthetic spectral tail only contributes $\sim$3% to the bolometric flux, which we found to be $6.35\times 10$-12 W m-2. We then derived the bolometric magnitude (m $_{\rm bol}$) using the Sun as a standard (adopting $L_{\odot} =3.86\times 10^{26}$ W and $M_{\rm bol\odot}=4.75$), which yields $m_{\rm bol}=9.02$. This is consistent with the value derived by Leggett et al. (2000) of 9.07. Assuming a distance modulus of $m = -3.11 \pm 0.01$ (van Altena et al. 1995) and 10% uncertainty in our flux calibration, we thus determine that GJ406 has $M_{\rm bol} = 12.13 \pm 0.10$ and $L =\log{L_*/L_{\odot}}=-2.95 \pm 0.05$.

3 Theoretical spectra computation procedure

Theoretical spectral energy distributions[*] were computed for model atmospheres of dwarfs with effective temperatures $T_{\rm eff} = 2500$-3200 K from the NextGen grid of Hauschildt et al. (1999) for solar metallicity (Anders & Grevesse 1989). Hereafter we use the syntax "effective temperature/gravity/metallicity'', e.g. 2800/5.0/0 to signify the model atmosphere. Unless otherwise mentioned all models are for $\log~ g = 5.0$. Computations of synthetic spectra were carried out by the program WITA6 (Pavlenko 2000) assuming LTE, hydrostatic equilibrium for a one-dimensional model atmosphere and without sources and sinks of energy. The equations of ionisation-dissociation equilibrium were solved for media consisting of atoms, ions and molecules. We took into account $\sim$100 components (Pavlenko 2000). The constants for equations of chemical balance were taken from Tsuji (1973).

Molecular line data were taken from different sources. Lines of 1H216O were computed using the AMES database (Partrige & Schwenke 1998). The partition functions of H2O  were also computed from these data (see Sect. 3.1). 12C16O and 13C16O  line lists were computed by Goorvitch (1994). The CO partition functions were taken from Gurvitz et al. (1989). TiO line lists were taken from Plez (1998) and Schwenke (1998). CN lines came from CDROM 18 (Kurucz 1993); CrH and FeH lines were taken from Burrows et al. (2002) and Dulick et al. (2003), respectively. Atomic line list was taken from VALD (Kupka et al. 1999).

The profiles of molecular and atomic lines were determined using the Voigt function H(a,v). Parameters of their natural broadening C2 and van der Waals broadening C4 were taken from Kupka et al. (1999) or in their absence computed following Unsöld (1955). Owing to the low temperatures in M dwarf atmospheres and consequently, electron densities, Stark broadening could be neglected. As a whole the effects of pressure broadening prevail. Computations for synthetic spectra were carried out with a step 0.5 Å  for microturbulent velocity vt = 1 - 4 km s-1. The instrumental broadening was modelled by gaussian profiles set to approximate the resolution of the observed spectra. The relative importance of the different opacities contributing to our synthetic spectra is shown in Fig. 2.

\end{figure} Figure 2: The contribution of different molecules to the formation of the synthetic spectrum of in $T_{\rm eff}$/log g = 2800/5.0 model atmosphere.
Open with DEXTER

3.1 Partition functions of water

We recomputed the constants of chemical equilibrium following Kurucz (1970) taking into account weights si = 1/4 and 3/4 for levels of water of different symmetry. We followed the scheme described by Pavlenko (2002). Let us write an equation of ionisation-dissociation equilibrium for the molecule consisting of x, y,..., z atoms as

$\displaystyle n_x*...*n_z/n_{x...z}= \exp(-E_{xy...z}/T_{\rm ev}+b-c*$      
$\displaystyle +3/2*(m-k-1)*{\rm ln} T)$     (1)

where E and T are dissociation energy and temperature (in eV), nzis the number density of z-species, k and m are ionisation degree (0 for neutrals) and number of atoms per molecule, respectively (see Kurucz 1970, for more details). Computed constants a,b,c,d,e,fare given in Table 2.

\end{figure} Figure 3: Partition functions of molecules H2O and HDO. Our data for H2O are compared with Vidler & Tennyson (2000), data for HDO are compared with Hewitt et al. (2005). Our partition function for HDO is computed for s1 =1. The differences in the HDO partition functions at higher temperatures (T > 4000 K) are due to the use of more complete sets of deuterated water levels in the UCL model compared with the AMES model.
Open with DEXTER

\end{figure} Figure 4: Molecular densities of H2O computed for NextGen model atmospheres 2800/5.0 and 3800/5.0 for the cases of partition functions computed with si = 1, 4 and si =1/4, 3/4.
Open with DEXTER

Table 2: Constants of the dissociation equilibrium of H2O in formats of the ATLASxx, where xx labels a version of ASTLAS, the superscripts correspond to different fit temperature ranges: 1 - $300~{\rm K} < T < 1000$ K, 2 - $300~{\rm K} < T < 6000$ K, 3 - $60~{\rm K} < T < 6000$ K, 4 -  $50~{\rm K} < T < 10~000$.

The temperature dependence of computed partition functions of H2O are given in Fig. 3. In general, our new partition functions agree well with data of Vidler & Tennyson (2000) computed using a mixture of experimental data and a UCL model of the water vapour molecule. Some differences occur at T > 5000 K as the UCL model has more levels of high excitation energy.

It is worth noting that

the molecular densities of given molecules obtained from a solution of the system of equations of molecular equilibrium response to changes of absolute values of U, especially in high temperature regimes (T > 3000 K, see Fig. 4).

U(H2O ) depends strongly on temperature. Table 2 provides fitting constants obtained for different temperature regions. From a general point of view it would be reasonable to restrict our fitting to temperatures T > 300 K. Though for some astrophysical objects it would be interesting to have the partition function of water vapour for even lower temperatures. Thus we have provided these as well.

4 Results

4.1 Dependence of theoretical spectra on input parameters

First of all, we computed model SED's using different input parameters: effective temperatures, gravities, metallicities, microturbulent velocities. Ratios of fluxes computed with different sets of input parameters are shown in Fig. 5. It can be seen that a temperature change of 200 K is roughly equivalent to a change in metallicity of 0.5 dex or a gravity change of $\Delta$log g = 1.

\end{figure} Figure 5: Responses of computed spectra to variations of input parameters. A model atmosphere of 3000/5.0/0.0 was used as the reference model.
Open with DEXTER

In Fig. 5 we see the differential effects of the dependence of our model spectra on different parameters. Some of these effects depend more on changes in the opacities for different parameters. However, some effects are explained by changes in the structures of the model atmosphere. Indeed, stellar photospheres of different $T_{\rm eff}$, log g, [M/H] lie in different pressure regions.

\end{figure} Figure 6: Top: temperature structure of model NextGen, DUSTY and COND model atmospheres of 2800/5.0/0. Bottom: ratio of fluxes computed for these model atmospheres.
Open with DEXTER

\end{figure} Figure 7: Ratio of fluxes computed for 2800/5.0 NextGen model atmosphere with TiO line lists by Plez (1998) and Schwenke (1998).
Open with DEXTER

\end{figure} Figure 8: Min $F(f_{\rm s}, f_{\rm n}$) found for different model atmospheres.
Open with DEXTER

4.2 Dependence of theoretical spectra on different model atmospheres

We also examined the dependence of our results on the choice of model atmospheres: NextGen, Dusty, Cond (see Allard 2005 for references). These model atmospheres have different T=f(P)structures (Fig. 6) due to differences in the physical treatment of the dust formation which cause some changes in the opacities and the molecular equilibrium. In general, the DUSTY and COND models have more hot inner layers and cooler outer layers (see first plot in Fig. 6).

We computed model theoretical fluxes for the NextGen, Dusty and COND model atmospheres with parameters 2800/5.0/0 in the spectral region of interest and compared them. Theoretical fluxes were convolved with a spectral resolution element of FHWM = 4 Å. In Fig. 6 we show the ratio of the convolved fluxes.

It is worth noting a few results:

Infrared spectra containing strong water bands agree rather well.

We see rather big differences in the optical and blue parts of the spectra. Here the dependence of computed spectra to changing temperature is much stronger in comparison to the IR spectral region. In general, both the COND and DUSTY model atmospheres are a bit hotter in the inner regions, and we therefore see an increased flux in the continuum flux levels (or increased flux for lower opacity wavelength ranges). Some molecular lines here become stronger due to the lower temperatures in DUSTY and COND model atmospheres.

CO bands strength responds to the change in temperature in the photospheric layers, due to the high sensitivity of dissociation equilibrium of this molecule to temperature.

4.3 Dependence of theoretical spectra on the use if different TiO line lists

There are two TiO line lists of widespread use by Schwenke (1998) and Plez (1998). They are based on the improved Langhoff (1997) model of the TiO molecule but differ in details. Plez (1998) added an a-f system at 0.5 $\mu$m to the line list. Schwenke subsequently computed a corresponding list of transitions complete to the higher excitation energies. Plez (1998) provided a line list for a solar mixture of Ti isotopes. Schwenke's (1998) database provides lists for each TiO isotope.

We compare synthetic spectra computed with these two line lists for a NextGen model atmosphere 2800/5.0/0 (Fig. 7). Synthetic spectra were computed with a 0.1 Å  step, then convolved with a spectral resolution element of ${\it FWHM} = 1$ Å.

We find the largest differences in the blue part of spectrum, which is more affected by incompleteness of molecular line lists for other molecules, chromospheric effects, veiling, and strong atomic absorption. Therefore from inspection of Fig. 7 we conclude that differences between synthetic spectra computed with Plez (1998) and Schwenke (1998) line lists do not affect our main results (see also Lyubchik & Pavlenko 2001).

4.4 Fits to GJ406 spectra

In this paper we are interested in the dependence of computed spectral energy distributions to $T_{\rm eff}$. To find the best fit of computed spectra to observed fluxes we use a minimisation procedure described in Jones et al. (2002) and Pavlenko & Jones (2003). Namely, the best fits are found for the ${\rm min}~F(f_{\rm s},f_{\rm n},
f_{\rm g})$, where $f_{\rm s}$, $f_{\rm n}$ and $f_{\rm g}$ are relative shift of the spectra, normalisation constant for the computed spectra, and broadening parameter, respectively. We found a rather weak dependence of the min F on the $f_{\rm g}$ broadening parameter and set $f_{\rm g} = 6$ Å. Previous studies have considered GJ406 as a typical M6 dwarf, and we thus assumed $\log~ g = 5.0$ in its atmosphere.

To determine a self-consistent solution we fit the theoretical spectra to the observed fluxes in all spectral regions and estimate the quality of the fit by computing $F(f_{\rm s}, f_{\rm n}$). Two spectral regions were excluded from our analysis: 0.35-0.4 $\mu$m due to incompleteness of our opacity sources, and 4.3-0.461 $\mu$m due to the gap in the observed data.

In Fig. 8 we show computed $F(f_{\rm s}, f_{\rm n}$) for a grid of our theoretical spectra of different $T_{\rm eff}$. We find a weak dependence of $F(f_{\rm s}, f_{\rm n}$) on $V_{\rm t}$. The best fit can be found for the min $F(f_{\rm s}, f_{\rm n}$) at $2800\pm 100$ K (Fig. 9 shows these fits in linear and logarithmic scales).

Some problems with fitting spectral features are seen at $1.3< \lambda < 1.7$ $\mu$m. Partially these discrepancies can be exlained by problems with modelling of strong H2O bands located here (Jones et al. 2005). Fits to TiO and VO in optical spectral regions and H2O  (beyond 1.7 $\mu$m) are of particularly good quality (see Pavlenko 1998; Jones et al. 2002 for more details).

We also made comparison of the observed spectrum with synthetic spectra based on the Cond model atmospheres and obtained the fits practically of the the same quality. The "best fit'' synthetical spectra computed with NextGen and Cond model atmospheres coincide over the whole spectral region. The largest differences do not exceed 1-2%. Consequently, the function $f(x_{\rm s}, x_{\rm n})$ for fits to the Cond synthetic spectra has the minimum at 2800 K, as for NextGen (see Fig. 10).

\par\includegraphics[width=14cm,height=6.4cm]{2979_10.eps} %
\end{figure} Figure 9: Fits with linear and logarithmic flux scales to our GJ406 spectrum with a theoretical spectrum computed for a solar composition NextGen model atmosphere 2800/5.0.
Open with DEXTER

To confirm this explanation of the discrepancy between theory and observation, it would be desirable to carry out a similar analysis of late-type dwarfs like GJ 406, but employing a more complete water vapor line list for the theoretical atmosphere calculations. Such a line list is expected to be available in the near future (Jones et al. 2005).

4.5 Evolutionary model fits

Using evolutionary models by Baraffe et al. (2003), we have estimated a mass and age for GJ 406 (Fig. 11). Our measured effective temperature $T_{\rm eff}$ =  $2800\pm 100$ K, and luminosity $L =\log{L_*/L_{\odot}}=-2.95 \pm 0.05$ agree well with estimates of other authors (Leggett et al. 2002; Jones et al. 2002). We know from the lack of lithium in GJ 406 that its age must be >108 yr (Magazzu et al. 1993). Indeed, within the error bars we find that the measured $T_{\rm eff}$ and luminosity are consistent with the Baraffe et al. (2003) evolutionary models for stars of $T_{\rm eff}$ = 2700-2900 K and L = -3.0, which complies an age of 0.1-0.3 Gyr, and, respectively, M = 0.07-0.94 $M_{\odot}$. For L = -2.9 we obtain the age 0.2-0.35 Gyr, and M = 0.1-0.09 $M_{\odot}$.

5 Discussion

A comparison of observed and computed spectral energy distributions provides a unique tool with which to assess the completeness and quality of our knowledge about the structure and properties of late type dwarf atmospheres; the physical state of their matter; opacity sources and line lists, atmospheric temperature and pressure structure, effective temperature scales.

Here we have modelled the spectral energy distribution of the M6 dwarf GJ 406 from 0.4-4.9 $\mu$m. The optical spectrum is formed primarily by absorption of the saturated bands of VO and TiO. In general, the response of optical fluxes to the variations of input physical parameters is stronger than at infrared wavelengths (see Fig. 5). However, the infrared spectrum is more sensitive to the H2O absorption bands, with fluxes coming from deeper atmospheric layers. Therefore, if one is to understand both the outer and inner atmospheric structure of cool dwarfs, one must simultaneously account for both the optical and infrared spectra. In general, we achieved good agreement between our theoretical spectra and observation, although we note some problems with fits at certain wavelengths (0.35-0.4 $\mu$m, 1.3-1.7 $\mu$m, these are presumably a consequence of missing molecular opacities). Our fits support the idea that the TiO/VO and H2O  line lists covering these wavelength ranges are of good quality. Most probably, our problems in the blue part of the spectrum can be solved with proper fits to strong atomic lines located there. Then, the blue part of the spectrum should be more affected by chromospheric like phenomena. The detailed analysis of these and related problems is beyond the scope of this paper. We plan to consider those in forthcoming papers.

\end{figure} Figure 10: A comparison of $f(x_{\rm s}, x_{\rm n}$) computed for the fitting ofobserved spectra by NextGen and Cond synthetic spectra.
Open with DEXTER

Our best fit age range is consistent with the age constraints from depleted lithium. However, GJ406 has been kinematically classified as an old disk star, and as such, one might expect its age to be greater than $\sim$600 Myr (the age of the Hyades, which traditionally represents an upper age limit for the young disk population; e.g. Leggett 1992). Our model fit age range thus suggest that GJ406 is more youthful than this, and that its old disk classification should be interpreted solely as a kinematic description, and not as an age constraint. This is not contradictory of course, since the dispersion in the kinematics of the young disk population naturally places some young stars outside of the canonical young disk UVW kinematic region.

Although the spectrum of GJ406 shows strong H$_{\alpha}$ emission and activity (as summarised in Sect. 1), it should be noted that one cannot use this to place any strong constraints on age. As Gizis et al. (2002) explain, one expects a significant spread in the H$_{\alpha}$ emission strength of a population if the age is less than some value that depends on the stellar colour or spectral type (see their Fig. 11). For GJ406 ( V-I = 4.06), this age upper limit is greater than the age of the disk, and a high level of H$_{\alpha}$ emission may thus be expected for some M 6 dwarfs spanning the full age of the disk. It is at least clear that the emission properties of GJ406 are not inconsistent with our age estimate.

\end{figure} Figure 11: Top: location of GJ 406 in respect to the Baraffe et al. (1998) evolutionary tracks and $T_{\rm eff}$ = 2800, 2600 K and $L = L_*/L_\odot =-2.90$, -3.00 planes. The planes shown on plot determine the cuboid of our errors of $T_{\rm eff}$ and L. Bottom: location of GJ 406 in respect to evolutionary tracks in coordinates L, t, M. The thick lines shown on plot determine the cuboid of our errors of age, luminosity $L_*/L_\odot $ and mass M*.
Open with DEXTER

The William Herschel Telescope and United Kingdom Infrared Telescope are operated for the Particle Physics and Astronomy Research Council (PPARC). ISO is an ESA project with the participation of ISAS and NASA funded from member states. HST is a NASA project funded in partly by ESA. This work was partially supported by the Royal Society and the Leverhulme Trust. The computations were in part performed using the resources of HiPerSPACE computing facility at UCL which is partly funded by the UK Partical Physics and Astronomy Research Council (PPARC). We thank Maria Rosa Zapatero Osorio for her helpful comments. We thank the anonymous Referee for useful remarks. This research has made use of the SIMBAD database, operated at CDS, Strasbourg, France.



Copyright ESO 2006