E Moraux

J Bouvier

J R Stauffer

D Barrado Y Navascués

J-C Cuillandre

The lower mass function of the young open cluster Blanco 1: from 30M Jup to 3M ⊙

Keywords: Stars, low-mass, brown dwarfs -Stars, mass function -Open clusters and associations, individual, Blanco 1

Aims. We performed a deep wide field optical survey of the young (∼ 100 -150 Myr) open cluster Blanco 1 to study its low mass population well down into the brown dwarf regime and estimate its mass function over the whole cluster mass range. Methods. The survey covers 2.3 square degrees in the I and z-bands down to I ≃ z ≃ 24 with the CFH12K camera. Considering two different cluster ages (100 and 150 Myr), we selected cluster member candidates on the basis of their location in the (I, Iz) CMD relative to the isochrones, and estimated the contamination by foreground late-type field dwarfs using statistical arguments, infrared photometry and lowresolution optical spectroscopy. Results. We find that our survey should contain about 57% of the cluster members in the 0.03 -0.6M ⊙ mass range, including 30-40 brown dwarfs. The candidate's radial distribution presents evidence that mass segregation has already occured in the cluster. We took it into account to estimate the cluster mass function across the stellar/substellar boundary. We find that, between 0.03M ⊙ and 0.6M ⊙ , the cluster mass distribution does not depend much on its exact age, and is well represented by a single power-law, with an index α = 0.69 ± 0.15. Over the whole mass domain, from 0.03M ⊙ to 3M ⊙ , the mass function is better fitted by a log-normal function with m 0 = 0.36 ± 0.07M ⊙ and σ = 0.58 ± 0.06. Conclusions. Comparison between the Blanco 1 mass function, other young open clusters' MF, and the galactic disc MF suggests that the IMF, from the substellar domain to the higher mass part, does not depend much on initial conditions. We discuss the implications of this result on theories developed to date to explain the origin of the mass distribution.

Introduction

Since the first brown dwarf discovery and confirmation in 1995 (Nakajima et al. 1995;Rebolo et al. 1995), new perspectives have opened regarding the formation of condensed objects in molecular clouds. Today more than a thousand brown dwarfs (BDs) are known but their mode of formation is still controversial and the theoretical framework describing the stellar and Send offprint requests to: Estelle.Moraux@obs.ujf-grenoble.fr ⋆ Based on observations obtained at the Canada-France-Hawaii Telescope (CFHT) which is operated by the National Research Council of Canada, the Institut National des Sciences de l'Univers of the Centre National de la Recherche Scientifique of France, and the University of Hawaii. Based on observations collected at the European Southern Observatory, Chile (ESO programmes 68.C-0233 and 71.C-0446). Some of the data presented herein were obtained at the W.M. Keck Observatory, which is operated as a scientific partnership among the California Institute of Technology, the University of California and the National Aeronautics and Space Administration. substellar formation process(es) is far from being satisfactory. How do brown dwarfs form ? Is there a lower mass limit for an object to be formed ? A way to tackle these questions is to determine the mass spectrum resulting from the stellar formation process, i.e. the initial mass function (IMF), down to the substellar regime.

Young nearby open clusters are ideal environments for such a purpose. Their members constitute a uniform population in terms of distance and age, the extinction is usually low and uniform, and their youth ensures that brown dwarfs are still bright enough to be easily detected and followed-up. Moreover, the well-known stellar population of the nearest open clusters complements the discoveries of brown dwarf members to yield a complete mass function from substellar objects up to massive stars.

In this contribution we present a study of Blanco 1, sometimes referred to by the name of its brightest member ζ Sculptoris. It is a young southern open cluster (α J2000 = 0 h 04 m 24 s , δ J2000 = -29 • 56.4 ′ ; Lyngå1987) located at an Hipparcos distance of 260 +50 -40 pc (Robichon et al. 1999) at a high galactic latitude (some 240 pc below the galactic plane, b = -79.3 • ), and has a low extinction (E(B -V) = 0.010). It has a low stellar density with ∼ 200 known stellar members spread over a 1.5 • diameter surface (Hawkins & Favata 1998) and its age is not very well constrained, around 100-150 Myr. Panagi & O'Dell (1997) found it is similar to or slighty older than that of the well-studied Pleiades cluster of 120 Myr, while other authors used an age of 100 Myr (e.g. Pillitteri et al. 2003).

The cluster's youth combined with its large distance from the galactic disc point to an unusual formation history. Ford, Jeffries & Smalley (2005) found also that the abundance pattern in Blanco 1 is quite exceptional with an average [Fe/H] close to solar (= +0.04 ± 0.04) but a subsolar [Ni/Fe], [Si/Fe], [Mg/Fe] and [Ca/Fe]. This lead the authors to suggest that the material from which the cluster formed was different from the local galactic disc ISM. It may have travelled some distance without having the chance to homogenize with the ISM and may also have been polluted by one or two unusual supernova events. All these peculiarities make Blanco 1 a very interesting target to test the dependence of the IMF on environmental conditions when comparing its present day mass function (MF) to those of similar age open clusters.

Blanco 1 was discovered by Blanco (1949) who noticed a small concentration of A0 stars having the same galactic latitude in this region. Since then several photometric studies have been performed, e.g. by Westerlund (1963), Epstein (1968), Eggen (1972), Perry et al. (1978), Lyngå & Wramdemark (1984), Abraham de Epstein & Epstein (1985) and Westerlund et al. (1988). Abraham de Epstein & Epstein (1985) published a relatively complete list of 260 F, G, and K stellar candidates. Kinematic studies based on radial velocities (Jeffries and James 1999) and proper motions (Pillitteri et al 2003) indicate that the contamination level of this photometric sample is about 35 to 40%. X-ray observations performed by Micela et al. (1999) yielded the detection of a few lower mass sources -M dwarfs -confirmed by proper motion. However, the Blanco 1 population at very low masses and in the substellar domain is still unexplored.

We performed a deep wide-field optical survey and followup observations to look for Blanco 1 brown dwarfs (BD) and very low mass stars (VLM) in order to estimate its mass function (MF) down to the substellar domain. The observations and data reduction are described in Section 2 and our results are given in Section 3. In Section 4 we compare our MF estimate to the Pleiades, other young clusters and the galactic disc MF to investigate the IMF universality. We discuss the implications for recently developed IMF theories.

Observations and data reduction

Optical survey

Blanco 1 has been observed using CFHT's 12k×8k optical CCD mosaic camera (Cuillandre et al. 2001) during two separate runs in September 1999 and December 2000. Over the two runs, a total of 7 non overlapping fields, each of size 28 ′ × 42 ′ on the sky, were obtained in the I and z filters (central wave-lengths ∼ 850 and 950 nm respectively). The covered fields are graphically shown in Fig. 1 and their coordinates are listed in Table 1. For each pointing and each filter, the observing sequence included one 10s and two 600s exposures. A total area of 2.3 sq.deg. was thus surveyed down to a detection limit of I ∼ z ∼ 24 on and around the cluster's nominal center. The seeing measured on images obtained in both filters spanned a range from 0.5 to 0.8 arcsec FWHM. 1). The small cross in the middle represents the nominal cluster center and the dashed circles have a radius of 0.5, 1 and 1.5 degrees respectively. Small dots show the location of low mass star candidates selected from the short exposures. Very low mass star and brown dwarf candidates selected from the long exposures are shown as large dots if confirmed by infrared photometry and/or optical spectroscopy (open circles if ruled out).

These observations have been performed with the same instrument and partly during the same nights as the Pleiades cluster survey described by Moraux et al. (2003). The reduction of the CFH12K mosaic images has been done likewise, chip by chip. Briefly, they were first overscan and bias subtracted and a dark map correction was applied. Supersky flats were computed combining all images obtained in each band during the whole course of the observing run, and they were normalized to a reference CCD in order to obtain a uniform gain over the whole mosaic. The same photometric zero-point can thus be used for all the CCDs. Images were then flat-field corrected and fringes were removed (see Magnier & Cuillandre 2004 for details of the reduction procedure). Then the pre-reduced pairs of 600s exposures of each field were finally combined to yield the final long exposure images in the I and z filters.

The photometric analysis was done using the SExtractor package (Bertin & Arnouts 1996). To detect the maximum number of sources we first added all pairs of I and z reduced images together and we used the automatic object-finding algorithm on these frames. We then computed the photometry of each detected object by PSF-fitting on the individual (i.e. I or z) images. Spurious detections (cosmic rays, bad columns) and extended sources were rejected on the basis of their FWHM. For both runs, 1999 and 2000, the zero point of the instrument in the I-band was derived using Landolt photometric standard fields (SA98 and SA113) observed in the same conditions and reduced in the same way as Blanco 1 fields. Several exposures of SA98 were obtained with an incremental offset of several arcminutes in RA so that common sets of standard stars were observed on every CCD of the mosaic. We thus checked that the photometric zero-point in the I-band was the same for each CCD of the mosaic. For the z filter we used unreddened A0 standard stars covered by the Landolt fields assuming their I and z magnitudes were the same. For I < 20, the resulting photometric rms ranges from 0.03 to 0.05 mag for both I and z filters, which corresponds to a colour accuracy between 0.04 and 0.07 mag in Iz. For fainter objects the accuracy is slighty worse but remains better than 0.07 mag in I and z, or 0.10 mag in Iz, down to I ∼ 22.5.

The astrometric solution of each image was calculated with the Elixir software [START_REF] Magnier | [END_REF]. We first computed it for the short exposures using the USNO2 catalog as reference and we used this solution as a first approximation for the long exposures. Since most of the USNO2 stars were saturated in the 600s images we could not use this catalog and we used the stars detected in both short and long exposures as astrometric references instead. The same technique was used by Moraux et al. (2003) and the resulting astrometric accuracy is about 0.3 arcsec which is good considering the CFH12K plate scale of 0.206 arcsec/pixel.

Infrared follow-up photometry

We observed in the near-infrared 90 of the 124 brown dwarf and very low mass star candidates selected from the long exposures of the optical survey (see section 3.2) to confirm their status.

We obtained K s -band images for 77 of them using the SOFI camera (1024 2 pixels, 0.288 arcsec/pixel) at ESO/NTT in August 2003. The total exposure time per object was in the range from 10 to 20 min, from the brightest (K s ∼ 16) to the faintest source (K s ∼ 20). The seeing was typically around 0.8 arcsec during the entire run and the conditions were most of the time photometric.

Thirteen other candidates have been observed with the near-infrared camera CFHT-IR (1024 2 pixels, 0.211 arcec/pixel) in November 2004. The observing conditions were very good, photometric with a seeing of about 0.6 arcsec.

In both cases, individual infrared images obtained for each object at 7 or 9 dithered positions were first median combined in order to obtain a sky image which was subtracted from the original images. Sky subtracted images were then flat field corrected using dome flats and averaged to yield the final K s image of the object. Aperture photometry was performed on the object using the IRAF/APPHOT package. Several photometric standards from Persson et al. (1998) were observed each night throughout the run and reduced in the same way as the Blanco 1 candidates. They were used to derive the K s -band extinction coefficient as well as the zero-point of the instrument (ZP = 22.46 ± 0.02 mag for SOFI and ZP= 22.73 ± 0.01 for CFHT-IR) and the final photometry is given in the K-band Las Campanas system. The resulting internal photometric precision on the K-band measurements of Blanco 1 candidates is typically 0.02-0.03 mag for K ≤ 18 and 0.05-0.1 mag for fainter objects. While most of Blanco 1 targets were observed during photometric conditions, a few were affected by thin cirrus at the end of the 4th and 5th nights of the SOFI run. In order to calibrate the images taken in such conditions, we re-obtained shorter exposure time images of these fields during the following photometric nights. For those few candidates, the photometric accuracy is not as good, of order of 0.10 mag.

Optical spectroscopy

In addition to the photometric data we obtained low resolution optical spectroscopy for the brightest Blanco 1 brown dwarf candidates. The observations were performed using FORS2 (R ∼ 1500) at the VLT in November 2001 and LRIS (R ∼ 1000) at Keck in December 2001 and October 2002. For the three runs the data reduction was done in the same way using the IRAF package. Briefly, the CCD images were overscan, bias and flatfield corrected, the sky background was subtracted, and the spectra were extracted using the command APALL. The spectra were then wavelength calibrated using HgNeAr spectra and corrected from the instrument response computed by dividing observed standard star spectra by their reference spectra.

Results

In the following we present the results we obtained from the analysis of the short (Sect. 3.1) and long (Sect. 3.2) optical exposure images. The corresponding detection limits are I ∼ 13 -21 for t int = 10s and I ∼ 17.5 -24 for t int = 600s which allow our study to be complete from I = 14 to I = 22.2. For both integration times, we constructed the (I, Iz) colourmagnitude diagram (CMD) from which we selected low mass star and brown dwarf candidates. In order to take into account the uncertainty on the cluster age we made two selections, one for an age of 100 Myr and the other for 150 Myr. For each sample, we estimate the contamination by field dwarfs using either statistical arguments or additional near-infrared photometry and optical spectroscopy when available. We then investigate the cluster member spatial distribution and determine the cluster core radius for different mass ranges (Sect. 3.3). We finally use these results to derive the Blanco1 mass function from 30 Jupiter masses to 3M ⊙ (Sect. 3.4).

Low

-mass candidate selection (0.09M ⊙ ≤ m ≤ 0.6M ⊙ )
The (I, Iz) CMD corresponding to the short exposures of the optical survey is presented Fig. 2. The 100 Myr and 150 Myr isochrones given by the NEXTGEN models from [START_REF] Baraffe | [END_REF] 1 (see Moraux et al. 2003). shifted to the cluster distance ((m -M) o = 7.1 ± 0.4; Robichon et al. 1999) are shown as a dashed and dot-dashed line. In order to avoid incompleteness due to detection bias (either saturation on the bright side or too low signal to noise ratio on the faint objects) we restricted our study to the magnitude range from I = 14 to 18.5.

We then made a rather conservative photometric selection for each age (100 Myr and 150 Myr) by shifting the corresponding NEXTGEN isochrone to the maximum cluster distance given by Hipparcos ((m -M) o = 7.5), and on the blue side by the colour photometric error given by:

σ I-z = 0.014 + 2.10 -8 exp(0.8 × I) (1)
All the objects redder than this line are Blanco 1 candidate members. The 150 Myr isochrone is slightly bluer than the 100 Myr one, which means that the second selection will pick up a few more sources. The objects selected by both the 100 Myr and 150 Myr isochrones (by only the 150 Myr isochrone) are shown as triangles (resp. large dots) on Fig. 2. An electronic list of all the candidates is available from the CDS website.

According to the NEXTGEN model, the magnitude range I = 14 -18.5 corresponds to a mass m between 0.09M ⊙ and 0.6M ⊙ at 100 Myr (0.10 -0.6M ⊙ at 150 Myr) for (m -M) o = 7.1. One of the major shortcomings of this photometric selection is the contamination by sources unrelated to the cluster which may lie in the same region of the CMD. While a fraction of the selected objects must be true cluster members, others are merely field objects (galaxies, red giants or older field stars) on the line of sight. We do not expect to find many galaxies as contaminants because we rejected all the extended objects through the photometric analysis. The number of red giants in our sample is also expected to be small because they would have to be at greater than 6 kpc above the galactic plane to overlap in the CMD with Blanco1 members, and the density of giants at that height is small. Instead most of the contaminants are expected to be foreground field M-dwarfs. 1 The NEXTGEN and DUSTY isochrones have been calculated especially for the CFH12K filter set by I.Baraffe using the transmission curves and the CCD quantum efficiency. Their reliability is ensured by the fact that they work well at reproducing the Pleiades single star sequence for an age of 120 Myr and a distance of 125 pc in agreement with the cluster properties We cross-correlated the whole list of 534 candidates with the 2MASS catalogue for K < ∼ 15.5. About 490 of them have a 2MASS counterpart, the others are too faint. The corresponding (I, I -K) CMD and (J -H, H -K) colour-colour diagram are shown in Fig. 3 and Fig. 4 respectively. Most of the candidates lie on or close to the 100 Myr or 150 Myr NEXTGEN isochrone in both diagrams. This indicates that there is little extinction towards Blanco 1 and that all the candidates have infrared colours compatible with being M dwarfs. This supports our contention that there are few giants or galaxies in our candidate sample. It shows also that infrared colours fail to distinguish between field dwarfs and cluster members in this mass range as expected from the models.

To estimate the number of contaminants in our sample, we use instead statistical arguments. For each (I, Iz) position in the CMD, we first estimate the absolute M I magnitude of a field dwarf having the same Iz colour using the 5 Gyr NEXTGEN isochrone, and we derive its distance from the difference I -M I . An apparent magnitude range I, I + ∆I yields a distance range and we compute the corresponding volume containing the contaminants from the surface of the sky covered by our survey. Using the M-dwarf M I -band luminosity function from Zheng et al. (2004) as well as their scale height2 we calculate the number of expected contaminants per M I bin in the magnitude range I, I + ∆I. Taking into account a gaussian photometric error of width σ I-z given by equation 1, we then derive the number of objects per Iz. Finally, to obtain the total number of field M dwarfs per magnitude bin in our sample, we integrate over Iz, for Iz redder than our selection limit. The results are summarized Table 2. Out of the 479 (534) candidates detected in our survey for an assumed age of 100 Myr (150 Myr), we estimate about 211 (292) contaminants, which corresponds to a contamination level of ∼ 45% (∼ 55%).

Brown dwarf and very low mass star selection

(0.03M ⊙ ≤ m ≤ 0.09M ⊙ )
The optical (I,Iz) CMD corresponding to the long exposures of the survey is shown on the left panel of Figure 5, as well as the 100 Myr and 150 Myr isochrones. These isochrones follow the NEXTGEN models [START_REF] Baraffe | [END_REF]) down to T e f f ≃ 2500 K, i.e. down to I ≃ 20.5 for Blanco 1. Then dust starts to appear in the brown dwarf atmosphere and the NEXTGEN models are not valid any more. When T e f f ≤ 2000 K (I ≃ 21.5), all the dust is formed and the isochrones follow the DUSTY models (Chabrier et al. 2000). In between there is a transition and we simply draw a line from the NEXTGEN to the DUSTY models to ensure the continuity of the isochrones.

To select the candidates we proceed in the same way as for the short exposures. For each age, we shifted the isochrones vertically to the largest cluster distance ((m -M) 0 = 7.5) and horizontally on the blue side by the colour photometric error. All the objects redder than this line are then considered as Blanco 1 VLM and BD candidates. This gives us two samples depending on which age we use. The line corresponding to the bluest selection (based on the 150 Myr isochrone) maximizes the number of candidates and is shown on Figure 5. We identi- 3. Their mass spans a range between ∼ 0.03M ⊙ and 0.1M ⊙ and depends on the age of the cluster. Taking into account the distance of the cluster, the stellar/substellar boundary is around I ≃ 19.1 at 100 Myr and 19.6 at 150 Myr according to the NEXTGEN models, which means that about 75-80 of the selected objects are brown dwarf candidates.

As explained in the previous section, this optical photometric sample is contaminated by late M-type field dwarfs and we therefore conducted follow-up observations of the candidates to confirm their membership. We obtained near-infrared photometry with SOFI at ESO/NTT for 77 of them, and with CFHT-IR for 13 others. The measured K-band photometry is given in Table 3. The candidate location in the (I, I -K) CMD (Fig. 5, right panel) compared to the isochrones shifted to (m -M) 0 = 7.5 allowed us to weed out 18 objects, too blue for being Blanco 1 members. They are indicated as open circles in Fig. 5. The I -K colour is a good diagnostic of membership in the DUSTY regime, for I 21, but not as much for the brightest objects because field dwarfs and cluster members have similar infrared colours for Iz 1.0 (see Fig. 6). A statistical estimate using the Zheng et al. (2004) luminosity function (LF) and calculated in the same manner as for the indicates that there should be 10-15 contaminants between I = 17.7 and 18.5. Assuming that we are about 40% complete in this magnitude range (because of saturation), this yields a contamination level of ∼ 30%. Below I = 18.5, we estimate the survey is complete and we find a contamination level of about 30% for I = 18.5 -20. Between I = 20 and 21, we cannot use the same method because we reach the end of the field M-dwarf LF but we do not expect the contamination to be very high as the cluster sequence is well detached from the bunch of field objects in the (I, Iz) CMD. The fact that we do not identify any contaminant (except maybe CFHT-BL-57) in the (I, I -K) CMD in this magnitude range confirms that the con-tamination level is indeed very low. For I 21 and Iz > 1.0 infrared colour is sufficient to assess membership. The I -K colour of a 5 Gyr field dwarf starts indeed to differ from that of a ∼ 100 Myr cluster member in this domain (see Fig. 6). Even though the models may be too blue in Iz, empirical data from Pleiades studies (Moraux et al. 2003, Casewell et al. 2007) indicate that the I -K colour difference is real. Pleiades members confirmed by proper motion are redder in I -K than older field dwarfs with the same Iz > 1.1. Additional spectroscopic data are thus especially needed for candidates brighter than I ∼ 20. We obtained optical low resolution spectroscopy for 17 of them with FORS2 on VLT and/or LRIS on Keck. The spectra are shown Fig. 7 and8.

From these data, we estimate the spectral type of each object using different methods. We first use the PC3 index defined as the ratio of the flux between 8230Å and 8270Å to the flux between 7540Å and 7580Å in order to derive an estimate of the spectral type via the relationship S pT = -6.685 + 11.715 × (PC3) -2.024 × (PC3) 2 given by Martin et al. (1999). We also compare our spectra to those of spectroscopic standard stars obtained in the same conditions which allows us to estimate the spectral type with a precision of one sub-class. By fitting the data with synthetic spectra from the DUSTY models by Allard et al. (2001), we also estimate the effective temperature of our candidates with an accuracy of 100K corresponding to the grid step. An example of such a fit is given in Fig. 9 for BL-50. All these methods are in agreement with each other and the results are given in Table 4.

Knowing the effective temperature of the candidates, we then compare the gravity sensitive Na I doublet (8183Å, 8195Å) profile to synthetic spectra of similar T e f f given by Allard et al. (2001) for various log g. The best fit allows us to estimate the surface gravity of the objects with an accuracy of 0.5 in log g corresponding to the grid step. Since the expected value of log g for ∼ 100 Myr brown dwarfs (log g ≃ 4.8) is lower than the one for older field dwarfs (log g ≃ 5.3), we can use this estimate as a membership criterion. Figure 10 shows the comparison of the Na I doublet profile of a very low mass candidate (BL-29) and a similar spectral type (M6) standard field dwarf GJ-866 with two synthetic spectra of similar T e f f for log g = 4.5 and 5.0. The older field object is better fitted by the spectrum of larger gravity as expected, while the best Table 3. Very low mass star and brown dwarf candidates selected from the long exposure images on the basis of their location in the (I,Iz) CMD (see Fig. 5). Objects selected by the 150 Myr isochrone but not the 100 Myr isochrone are written in italic. An electronic version of this table can be found on CDS. fit of BL-29 profile is obtained for log g = 4.5. This indicates that this object is likely a cluster member. We did this analysis for all the candidates for which we have a spectrum with a good enough signal to noise ratio to estimate their surface gravity. When possible we also measure the equivalent width of the H α emission line (6563Å), indicative of chromospheric activity. All the results are indicated in Table 4. All the candidates observed spectroscopically are late M dwarfs which confirms that the contamination by galaxies or red giants is indeed very low. CFHT-BL-44 has an earlier PC3 spectral type and a hotter effective temperature than the other candidates of similar I-magnitude. Moreover its surface gravity is too large to be a cluster member. It is most likely a foreground field dwarf. We suspect CFHT-BL-28 is also a field dwarf for the same reasons (slighty earlier PC3 spectral type and hotter temperature) even though the analysis of the Na profile did not provide a definitive result. CFHT-BL-34 and -41 have a later spectral type and a lower temperature than ex-Fig. 7. FORS2 spectra of CFHT-BL VLM candidates whose number is given above each spectrum on the left hand side. pected from their I-magnitude and both Iz and I -K colours redder than the NEXTGEN isochrone. Their surface gravity is low and consistent with cluster membership and we therefore believe these two objects are cluster binaries. Note that three other objects (CFHT-BL-52, -53, -54) are also suspected to be cluster binaries from their location above the isochrone in the infrared CMD but we could not get any spectra as they are too faint.

CFHT I I -z I -K s RA J2000 DEC J2000 CFHT I I -z I -K s RA J2000 DEC J2000 -BL- (h m s) (d m s) -BL- (h m s) (d m
Between I = 18.5 and 20, three objects out of 34 photometric candidates have been identified as non-members, either because of their infrared colour (CFHT-BL-47) or because of their spectroscopic properties (CFHT-BL-28 and 44). Eight objects do not have any follow-up observation, neither K photometry nor spectroscopy, and we cannot conclude on their nature. This yields a contamination level between 10 and 30%, in agreement with the statistical estimate given above. We can refine this analysis if we note that all the objects without any follow-up have about the same Iz ≃ 0.74 and form a small group in the CMD (see Fig. 5), slighty bluer than the candidates that are confirmed spectroscopically. Three other objects that have infrared follow-up lie also in this group. All of them are bluer than the 100 Myr (and sometimes also than the 150 Myr) isochrone in I -K, and one is a non-member (CFHT-BL-28) based on its spectroscopy. Therefore the contamination level in this subgroup is probably high, larger than 50%, which gives a total contamination level around 20-30% for I = 18.5 ≤ I ≤ 20.

We can use the same argument to estimate the number of contaminants in the magnitude range I = 21 -22.2 corresponding to our sample completeness limit. We do not have any K-band photometry for nine objects that are amongst the bluest (Iz ≃ 1.0) candidates in this region of the (I, Iz) CMD. Considering that about half of the objects with Iz = 1.05 -1.10 have been identified as non-members (CFHT-BL-73, -79, -84, -95, -96) and that there is no contaminant amongst the redder objects, we expect that about five or more of these blue candidates are field dwarfs.

As for the short exposures, we did the full analysis for both selections based on the 100 Myr and 150 Myr isochrones respectively. The number of estimated contaminants and of probable cluster members per magnitude bin in the long-exposures images is summarized in Table 5 for each age.

Spatial distribution

Now that we have identified the cluster member candidates with 0.03M ⊙ ≤ m ≤ 0.6M ⊙ in the survey, we need to study their spatial distribution. It is necessary to check whether mass segregation has occured in the cluster in order to take it into account when determining the mass function.

Figure 11 shows the spatial distribution of the low mass star candidates selected in the short exposures when assuming a cluster age of 150 Myr. Using the kernel method described by [START_REF] Silverman | Monographs on Statistics and Applied Probability[END_REF] we estimated the isodensity contours. The peak of the distribution allows us to determine the coordinates of the cluster center: α 0 = 0 h 04 m 1.2 s , δ 0 = -30 • 12 ′ . This is about 0.3 • away from the nominal center given by [START_REF] Lyngå | Catalogue of open cluster data[END_REF] but close to the recent estimate given by Kharchenko et al. (2005) who used the maximum of the projected stellar density as the definition of the cluster center.

We then calculate the number of candidates per square degree located within each annulus R, R + dR centered on the new center (α 0 , δ 0 ). The radial profiles obtained for the low mass star candidates (0.09M ⊙ ≤ m ≤ 0.6M ⊙ , short exposures) Table 4. Measurements made on the optical spectra of Blanco 1 candidates. In the few cases where an object has been observed with both FORS2 and LRIS, the two calculated PC3 values are given in column 2 and the spectral type given in column 3 has been derived using the average value. and for the very low mass star and brown dwarf candidates (0.03M ⊙ ≤ m ≤ 0.09M ⊙ , long exposures) are shown Figure 12. We assume that the distribution of the low mass sample contaminants is uniform and we note their spatial density n cont . The low mass star distribution is then fitted by a King profile (King 1962)

plus n cont n(x) = k 1 √ 1 + x - 1 √ 1 + x t 2 + n cont (3)
where k is a normalisation constant, x = (r/r c ) 2 and x t = (r t /r c ) 2 . The distance from the cluster center (α 0 , δ 0 ) is r, r c the core radius and r t the tidal radius.

According to Pinfield et al. (1998), the tidal radius r t of a cluster close to the Sun in a circular orbit is given by

r t = GM c 2(A -B) 2 1/3 = 1.46 M 1/3 c (4)
where M c is the cluster mass and A and B are the Oort constants. Blanco 1 contains about 180 F, G and K stars (Abraham de Epstein & Epstein 1985, Jeffries & James 1999, Pillitteri et al. 2003), i.e. about half as many as the Pleiades. If we assume that both clusters have the same mass distribution, then we expect Blanco 1 is about half as massive as the Pleiades and r t (Blanco1) = 0.794×r t (Pleiades). The Pleiades tidal radius being 13.1 pc (Pinfield et al. 1998), we derive r t (Blanco1) = 10.4 pc or 2.3 • for a distance of 260 pc.

Using this value, the best fit minimizing χ 2 gives n cont = 112 ± 10, k = 430 ± 40 per square degrees and r c = 0.61 +0.16 -0.11 degrees, i.e 2.8 +0.7 -0.5 pc, for the 100 Myr sample. For 150 Myr, we find n cont = 125 ± 10, k = 479 ± 43 and r c = 0.79 +0.15 -0.16 degrees or 3.6 ± 0.7 pc. Note that the values obtained for n cont are consistent with the number of contaminants given in Table 2 divided by the surveyed area (2.3 sq.deg.). We find n cont = 92 ± 13 sq.deg -1 for the 100 Myr selection and n cont = 127 ± 17 sq.deg -1 for the 150 Myr one. This validates a posteriori our statistical contamination estimate. We then compute the total number N of cluster members in the mass range 0.09-0.6M ⊙ by integrating equation 3 for both ages and by adopting the mean value. We find N = 380 ± 50 stars (N 100Myr = 343 +30 -35 and Fig. 9. Spectrum of BL-50 (solid line) fitted by a synthetic spectrum (dashed line) from the DUSTY models of Allard et al. (2001). The best fit gives T e f f = 2700K, which is in agreement with the M7-M7.5 spectral type we find using the PC3 index or the comparison to spectroscopic standards.

N 150Myr = 416 ± 37). According to table2, about 255 ± 60 of them are present in our survey, i.e. ∼ 67%.

The radial distributions of the 49 (for the 100 Myr selection) and 51 (150 Myr selection) probable VLM & BD members remaining after analysis of the follow-up observations are also fitted by a King profile (see Fig. 12). Note that no constant is added to the King distributions as the samples used to plot the histograms have been corrected from the contamination (the candidates without any NIR photometry or optical spectroscopy are not taken into account here). The fit parameters corresponding to the 100 Myr selection are r c = 1.56 +0.16 -0.15 degrees (7.1 ± 0.7 pc) and k = 160 ± 23 sq.deg. -1 . For the 150 Myr sample, we obtain r c = 1.70 +0.16 -0.14 degrees and k = 190 ± 27 sq.deg. -1 . The integration of both fits yields 104 ± 15 and 112±15 objects respectively, indicating that we identified about 47% of them in the CFH12K images. Among the 17 candidates without any follow-up observations, we estimate that ∼ 9 are probable members (see previous section). Assuming that they follow the same radial distribution, this yields a total number of 125 ± 20 VLM & BDs in the cluster.

About 47% of the Blanco 1 very low mass stars and brown dwarfs have been covered by the survey whereas ∼67% of the low mass star members are present in our study. This difference has to be taken into account when estimating the cluster mass function. It also suggests that mass segregation due to two body interactions has already occured and that Blanco 1 is dynamically relaxed. The relaxation time t r is 5) [START_REF] Spitzer | [END_REF], where N c is the total number of cluster members, r h the half mass radius in pc and mc the average stellar mass. The crossing time t cr is then given by King (1980) by

t r = 0.8 √ N c r 3/2 h √ mc (log N c -0.3) Myr (
t cr = t r 31 ln N c 2 N c . ( 6
)
As N c (Blanco 1) < N c (Pleiades), r h (Blanco 1) 3 pc while r h (Pleiades) = 3.66 pc (Pinfield et al. 1998), and mc (Blanco 1) ∼ mc (Pleiades), we find that these two dynamical times are smaller for Blanco 1 than for the Pleiades. In particular, this yields t r (Blanco 1) ≤ 90 Myr which is smaller than its age around 100-150 Myr. This suggests that Blanco 1 is indeed dynamically relaxed. Moreover, when complete equipartition of energy is reached, the cluster core radius varies with stellar mass as m -0.5 , and we find that this is roughly the case here. For

Mass function

The number of probable Blanco 1 members per I-magnitude bin is given in Tables 2 and5, last column. In order to estimate the cluster system mass function -multiple systems are not resolved at the distance of the cluster -we first convert the magnitude bin to a mass bin using the 100 Myr or 150 Myr massmagnitude relationship from the NEXTGEN models [START_REF] Baraffe | [END_REF] for I ≤ 20.0 and DUSTY models (Chabrier et al. 2000) for I ≥ 21.0. The obtained mass ranges are indicated in Table 2 and5, second column. Then we divide the number of probable members by a correction factor corresponding to the percentage of low mass stars and brown dwarfs covered by our survey (i.e. respectively 0.65 and 0.48). The total number of cluster systems per unit mass is thus derived over a continuous mass range from 0.03 to 0.60M ⊙ and is shown in Fig. 13. In this mass domain, the mass function can be fitted by a single power-law dn/dm ∝ m -α . A linear regression through the data points yields α = 0.67 ± 0.14 for the 100 Myr sample and α = 0.71±0.13 for the 150 Myr sample where the uncertainty is the 1σ fit error. (In this representation, Salpeter's IMF corresponds to α = 2.35.)

At larger masses, we completed our study with data from the literature. Pillitteri et al. (2003) established a list of reliable cluster members based on proper motion measurements. Their sample contains objects with mass ranging from 3M ⊙ down to ∼ 0.2M ⊙ but is not complete below ∼ 0.4M ⊙ . Using the mass Fig. 12. Radial distribution of the low mass star candidates (upper panel) and the VLM & BD probable cluster members (lower panel). The cluster center has been estimated using the isodensity contours (Fig. 11). In both cases, the filled dots and solid histogram represent the 100 Myr sample while the empty dots and dashed line histogram correspond to the 150 Myr sample. The best χ 2 fit of a King profile is given for each histogram. Fig. 13. The Blanco 1 system mass function across the stellar/substellar boundary. Note that all the data points are derived from the same survey, using short exposures for the stellar domain and long exposures for the substellar regime. This provides a consistent determination of the slope of the cluster's mass function in the mass range from 0.030 to 0.6M ⊙ . The data points are fitted by a power law dn/dm ∝ m -α with an index α = 0.67 ± 0.14 for the 100 Myr sample (solid line) and α = 0.71 ± 0.13 for the 150 Myr sample (dashed line). estimate given in their paper (table A1)3 , we compute the number of objects per unit mass and multiply it by a constant, so that the number of stars in the range 0.4-0.6M ⊙ is the same as the one we derive from our survey. On the whole mass range, from about 30 Jupiter masses to 3M ⊙ , the Blanco 1 system mass function corrected for contamination is reasonably well fitted by a Scalo-like log-normal distribution

ξ L (m) = dn d log m ∝ exp - (log m -log m 0 ) 2 2σ 2
where m 0 corresponds to the peak of the distribution and σ represents its width (see Fig. 14). A χ 2 fit gives m 0 = 0.34 ± 0.05M ⊙ , σ = 0.58 ± 0.06 for the 100 Myr sample, and m 0 = 0.38 ± 0.05M ⊙ , σ = 0.58 ± 0.06 for the 150 Myr sample. By integration of the log-normal fit, we find a total number of brown dwarf systems (m = 0.01 -0.072M ⊙ ) in the cluster of about 95 (85), for ∼ 660 (690) stars in the mass range 0.072 -3M ⊙4 , which yields a ratio of about 15% (12%) for the 100 Myr (resp. 150 Myr) sample. The cluster total mass is ∼ 410M ⊙ for an age of 100 Myr (∼ 450M ⊙ for 150 Myr) with a substellar contribution corresponding to only ∼ 1%. The stellar mass of the Pleiades is 735M ⊙ (Pinfield et al. 1998), which confirms a posteriori our initial assumption of Blanco 1 being about half as massive as the Pleiades. From equations 5 and 6, we find t r 60 Myr and t cr 15 Myr, which indicates that the cluster is dynamically relaxed indeed.

Discussion

Blanco 1 has an age similar to, or a little older (Panagi et al. 1997) than, the Pleiades but it is less dense (about 30 stars/pc 2 against ∼ 65 stars/pc 2 ), about half as massive and has a different abundance pattern (Ford, Jeffries & Smalley 2005). Comparing the mass functions of these two clusters thus allows us to test the dependence of the IMF on environmental conditions, especially as they have been obtained from similar data sets and analysed in the same way. Both mass functions are plotted on Figure 14. The Pleiades MF is from Moraux et al. (2003), their data points are between 0.03M ⊙ and 0.5M ⊙ and the log-normal fit goes from 0.03M ⊙ to 3M ⊙ .

As seen from both Figures 13 and 14, the Blanco 1 mass function does not depend much on its exact age. The results are similar within the uncertainties. In the following we will do the comparison with the Pleiades MF for the 100 Myr Blanco 1 MF only but the conclusions remain valid for 150 Myr.

A quick look at Figure 14 is sufficient to note that the overall shapes of the Pleiades and Blanco 1 MF are very similar. Both mass functions can be fitted by a log-normal distribution with a characteristic system mass around m 0 ≃ 0.3M ⊙ and σ ≃ 0.5. It is interesting to add that this is also the case for a large number of star forming regions and other open clusters, as well as for the system galactic disk mass function down to 0.1M ⊙ at least (see Chabrier 2003). This suggests either that the IMF does not depend much on environmental conditions or that star formation occurs only for a narrow range of local parameters as discussed in section 4.2.

The M7/8 gap

By looking at the Blanco 1 and Pleiades data points, we also notice that for both clusters the point around 0.035M ⊙ is much above the log-normal fit. This cannot only be due to an underestimate of the contamination in this mass range as it is taken into account in the error bar. It is as if the IMF rises again below 0.04M ⊙ . Similarly, Muench et al. (2002) found a significant secondary peak close to the deuterium limit in the Trapezium IMF. However, the fact that this feature does not occur around the same mass but the same spectral type (later than M8) brought Dobbie et al. (2002) to suggest that this is not real. They argue instead that this structure reflects a sharp local drop in the luminosity-mass (L-M) relationship due to the onset of dust formation in the atmosphere around T e f f = 2700 K. A change in the slope of the L-M relationship implies that the mass of objects with spectral type later than M7-M8 may be significantly underestimated by the current NextGen and Dusty models. Thus, the number of objects in the lower mass bin may be overestimated, resulting in an artificial secondary peak in the IMF.

By applying the empirical magnitude-mass relation given by Dobbie et al. (2002;their Fig.3) to the Blanco 1 and Pleiades luminosity function, we find that the lowest mass point goes down for both MFs whereas higher mass points go slightly up (see Fig. 15), following nicely the log-normal fit. Even if this dust effect has still to be investigated, we consider the lognormal distribution is a good approximation to the observed MF, including in the substellar domain down to 30 Jupiter masses.

Fig. 15.

The Blanco 1 and Pleiades mass function obtained after applying the empirical magnitude-mass relation given by Dobbie et al. (2002). The open squares are the data points which have been modified while the filled circles are the data points which did not change. The log-normal fit for each cluster is shown as a solid line.

Dependency of the IMF on environmental conditions

With an age of about 100 Myr, the Blanco 1 and Pleiades clusters are already dynamically evolved in the sense that they are well relaxed and that mass segregation has occured. However they are still young enough as not to have lost many objects. Adams et al. (2002) estimated that only ∼ 10% of the primordial low mass star population has been lost in the Pleiades, which has been confirmed by other studies combining numerical simulations with observations (e.g. Moraux, Kroupa & Bouvier, 2004). Therefore we can assume that the present day mass functions discussed here represent well the IMF for both clusters.

We have already mentioned that, despite the different local parameters (stellar density, abundance pattern), these 2 IMFs are similar within the uncertainties over 2 decades in mass. This seems also to be the case for other young open clusters such as Alpha-Per (∼ 80 Myr, Barrado et al. 2002), IC2391 (∼ 50 Myr, Barrado et al. 2004a), NGC2547 (∼ 30 Myr, Jeffries et al. 2004) or IC4665 (∼ 30 Myr, De Wit et al. 2006). One could argue that only one "type" of clusters, rich and compact enough, can survive several tens of Myr. Therefore they cannot be very different initially and it is not very surprising to find a similar IMF. However most of the studied star forming regions (λ-Ori, Barrado et al. 2004b;IC348, Luhman et al. 2003a;Trapezium, Muench et al. 2002;σ-Ori, Béjar et al. 2004) also present similar IMF's within the uncertainties. Again, it is true that they are all fairly rich but even in Taurus, a much sparser region, recent results from Guieu et al. (2006) seem to indicate that the IMF is not very different (but see also Luhman et al. 2003b andBriceno et al. 2002 for a differing view). Probably even more remarkable is the similarity with the galactic disk mass function (Chabrier 2003) -although uncertainties remain in the substellar domain -as field stars come from a large variety of environments, from small groups with N star < 10 to rich clusters (e.g. [START_REF] Allen | to appear in "Protostars and Planets V[END_REF]. All these results suggest that the IMF, from the substellar regime to the high mass domain, is similar in all studied regions. The fact that the IMF does not seem to depend much on the environment brings strong constraints on star formation theories. To date, two main models (turbulent fragmentation and accretion/ejection) are proposed to explain the origin of the IMF; they are discussed below.

Turbulent fragmentation

Turbulent fragmentation can generate a distribution of density structures in molecular clouds due to supersonic shocks that compress the gas. Multiple compressions result in the formation of sheets and then filaments, whose density and width are due to the MHD shock conditions. High velocity shocks produce high density but thin filaments resulting in low-mass clumps, whereas low velocity shocks produce low-density but large clumps which account for high mass objects. Using the power spectrum of velocities from numerical simulations of turbulence, Padoan & Nordlund (2002) derive a clump mass distribution resulting in a IMF slope which closely matches Salpeter's. At lower masses however, only a fraction of these clumps are sufficiently dense to be unstable which results in a flattening and a log-normal shape in the substellar regime.

A major success of this model is that it is able to reproduce directly the observed prestellar core distribution (Motte & André 2001) which resembles effectively the stellar IMF at least down to 0.1M ⊙ . However, how the clump spectrum translates into the stellar mass distribution is not fully understood. In particular, it is very likely that the higher mass clumps will subfragment into smaller clumps and form several stars. Moreover, SPH simulations find that most of the small cores are unbound and do not collapse to form stars (Klessen et al. 2005).

An important result of the turbulent fragmentation model is that the shape of the IMF depends on the turbulence via the Mach number or the power spectrum [START_REF] Ballesteros-Paredes | Protostars and Planets V[END_REF], Goodwin et al. 2006) but also on the molecular cloud density (see Padoan & Nordlund 2002, their fig. 2). This is easy to understand using simple considerations. Typically, a denser cloud forms denser but smaller clumps yielding a larger num-ber of low mass cores. Conversely, at lower density, the mass function peaks at higher mass. This result concerns the core mass function and it is not completely clear yet how this will affect the stellar IMF. However it is reasonable to think that the trend will be the same, which is in contradiction with the observed invariance of the IMF.

Accretion/ejection paradigm

Another model which has been developed to explain the origin of the IMF origin of the IMF gives less importance to turbulence and more to gravity. According to Bate & Bonnell (2005), the role for turbulence is to generate structure in molecular clouds which provides the seeds for gravitational fragmentation to occur down to the opacity limit. Stellar masses are then set by a combination of accretion and dynamical ejections which terminate the accretion. In their model based on the results of hydrodynamical SPH calculations, all the objects begin with a mass corresponding to the opacity limit (∼ 0.003M ⊙ ) and grow by accretion until they are ejected from the dense core by dynamical interactions in unstable multiple systems. The objects that end up as brown dwarfs stop accreting before they reach stellar masses because their are ejected soon after their formation. In contrast the higher mass stars are the objects that remain in the dense gas and accrete for longer. The accretion rates of individual objects are drawn from a lognormal distribution and the dynamical ejection of protostars from a multiple system is described using a characteristic time-scale. Taking values for these parameters from the hydrodynamical simulations, the model is able to reproduce a reasonable IMF in the tellar domain but produces more brown dwarfs than are actually observed.

In the accretion/ejection model the characteristic (median) mass of the IMF is approximatively given by the product of the typical accretion rate and the typical time-scale for dynamical ejection. Bate & Bonnell (2005) and Bate (2005) found that that the characteristic mass varies linearly with the mean thermal Jeans mass M Jeans and depends therefore on both density ρ and temperature T . As M Jeans is smaller at higher density, the mass function is expected to peak at smaller mass. Indeed in a denser cloud dynamical ejections are more likely to occur, yielding the formation of a larger number of low mass stars. Variations of the IMF in different environments should then occur in the location of the peak and in the substellar regime, which is not supported by observations.

Thermal physics

The two main models developed to explain the origin of the IMF predict a dependence of the IMF on the environment, in particular on the density and temperature of the molecular cloud via the thermal Jeans mass and/or turbulence.

Observational results however indicate that the IMF is similar in many regions from the substellar regime up to the massive star domain. This suggests that either the global properties of all the molecular clouds in which cluster and field stars form are very similar -which seems very unlikely -or that the phys-ical process for fragmentation by itself imposes specific local conditions. [START_REF] Larson | [END_REF] has suggested that thermal physics in molecular clouds may play this role as it can yield a Jeans mass of order ∼ 1M ⊙ at the point of fragmentation independently of the exact initial conditions. This results from how the cooling rate changes with density. At lower densities the gas cooling is dominated by atomic and molecular line emission whereas at higher densities the gas is coupled to the dust and dust cooling dominates. This results in a barotropic equation of state involving a cooling term at low densities followed by a gently heating term once dust cooling prevails. When used in numerical calculations, this equation of state can set the characteristic mass scale for fragmentation and produce a realistic IMF (Jappsen et al. 2005, Bonnell et al. 2006a).

A prediction of thermal physics is that a different metallicity may yield a different IMF as the cooling equation of state is modified. If a molecular cloud is more metal rich, the cooling is more efficient and the temperature is smaller. Therefore M Jeans is smaller and the characteristic mass is shifted towards the lower mass. However, the thermal physics invoked by [START_REF] Larson | [END_REF] links density and temperature. Thus for a higher metallicity, the Jeans mass occurs at lower T but also lower ρ and whether it should be larger or smaller is not clear yet. It is therefore premature to give any trend for the dependency of the IMF on metallicity but it could be a good observational test to measure and compare the mass function of stellar clusters with different metallicity.

4.2.

4. An hybrid model ? Bonnell et al. (2006b) have suggested an "hybrid" model to explain the origin of the mass distribution. Turbulence is necessary to generates the filamentary structure in the molecular clouds which facilitates fragmentation. Thermal physics sets the mean Jeans mass for gravitational fragmentation which corresponds to the characteristic stellar mass of the IMF. The broad peak can then be understood as being due to the dispersion in gas densities and temperature at the point where fragmentation occurs. The higher mass part is due to the continued competitive accretion in dense cores while the lower mass IMF is ascribed to fragmentation and then ejection. Some of the collapsing regions, especially in filaments and circumstellar discs, sub-fragment due to the increase of gas density. The lower mass clumps are ejected soon after they form from their natal environment by dynamical interactions, stop accreting and remain low mass objects (very low mass stars and brown dwarfs).

This model seems to work at reproducing a realistic stellar IMF as the location of the peak and the higher mass part do not vary with initial conditions as shown by Bonnell et al. (2006a). However some issues remain in the substellar domain. In particular, the lower mass part of the IMF is still expected to depend on the initial density -denser clouds enhance dynamical ejections and thus form more brown dwarfs -which is in contradiction with the current observations.

Conclusion

We performed a deep large-scale photometric survey of the young open cluster Blanco 1 to study its low mass population, from about 30 Jupiter masses to 0.6M ⊙ . We selected cluster candidates on the basis of their location in CMDs compared to theoretical isochrones and we estimated the contamination using statistical arguments, infrared photometry and optical lowresolution spectroscopy when possible. We did this analysis for two cluster ages (100 and 150 Myr) and we find similar results in both cases. We estimate that about 300 cluster members have been covered by our survey, amongst which 30-40 may be brown dwarfs. The study of the candidate radial distribution indicates that this corresponds to ∼ 57% of the cluster low mass population.

It also suggests that mass segregation has already occured in the cluster and we took it into account when estimating the mass function. We find that a single power-law dN/dM ∝ M -α with α = 0.69 ± 0.15 provides a good match to the cluster MF accross the stellar/substellar boundary, in the 0.03 -0.6M ⊙ mass range for both ages. When complementing our survey with literature data from Pillitteri et al. (2003) in the higher mass range, we find that the whole cluster mass function, from 0.03M ⊙ to 3M ⊙ , is well fitted by a log-normal distribution

ξ L (m) = dn d log m ∝ exp - (log m -log m 0 ) 2 2σ 2
with m 0 = 0.36 ± 0.07M ⊙ and σ = 0.58 ± 0.06. This result is very similar to the Pleiades MF given by Moraux et al. (2003), the two clusters having about the same age but different richness. Similar MF shapes are found for other young open clusters, star forming regions and also for the galactic disc population. This suggests that the IMF, from the substellar regime up to the higher mass domain, is fairly insensitive to initial conditions and that there is a characteristic system mass around 0.3M ⊙ for star formation. Theories developed to explain the origin of the mass distribution, such as turbulent fragmentation (Padoan & Nordlund 2002) or the accretion/ejection paradigm (Bate & Bonnell 2005), predict however a dependence of the IMF on density and temperature via the thermal Jeans mass. This seems to be in contradiction with the observational results unless thermal physics as invoked by [START_REF] Larson | [END_REF] sets the characteristic mass for fragmentation independently of the initial conditions.

A clue to this problem may be found in studying the global properties of collapsing clouds. The fact that all known embedded clusters, with N star = 10 to more than 1000, have a stellar surface density which varies by a factor of only a few [START_REF] Allen | to appear in "Protostars and Planets V[END_REF] suggests that the initial conditions may not be very different. Investigation of the mass distribution of very low-mass pre-stellar cores should clarify this issue (e.g. Li al. 2007, Nutter & Ward-Thompson 2007).

Fig. 2 .

 2 Fig. 2. (I, Iz) colour-magnitude diagram for the short exposures of the optical survey. The dashed (dot-dashed) line is the 100 Myr (resp. 150 Myr) isochrone from the NEXTGEN models of Baraffe et al. (1998) shifted to the Blanco 1 distance. The triangles and large dots correspond to our photometric selection of cluster member candidates for both ages (see text). Typical error bars are shown on the right side.

Fig. 3 .

 3 Fig. 3. (I, I -K) colour-magnitude diagram for all the optically selected candidates which have a 2MASS counterpart. The dashed and dot-dashed line are the 100 Myr and 150 Myr isochrones from the NEXTGEN models of Baraffe et al. (1998) shifted to the Blanco 1 distance. The symbols are the same as in Figure 2.

Fig. 4 .

 4 Fig. 4. (J -H, H -K) colour-colour diagram of all the low mass star candidates with a 2MASS counterpart. The different symbols correspond to different flags in the 2MASS catalog, i.e. different photometric accuracy (dots = flag A, open circles = flag B, crosses = flag C). A typical error bar is represented for each symbol. The 100 Myr, 150 Myr and 5 Gyr (corresponding to the field dwarf locus) isochrones are shown as a dashed, dot-dashed and a dotted line respectively and are labelled with mass. The reddening vector is also shown. All candidates are located on or close to the isochrones -within the error barsindicating that there is little extinction and that they all have near-infrared colours consistent with M-dwarf spectral type.

Fig. 6 .

 6 Fig. 6. (Iz, I -K) colour-colour diagram of the very low mass star and brown dwarf candidates. The symbols are the same as in Fig. 5. The 100 Myr, 150 Myr and 5 Gyr isochrones are shown as a dashed line, dot-dashed line and solid line respectively. The asterisks (resp. open squares) represent Pleiades candidates confirmed (resp. rejected) by proper motion(Moraux et al. 2003, Casewell et al. 2007). Only a few error bars are drawn for clarity.

Fig. 5 .

 5 Fig. 5. (I,Iz) and (I,I -K) CMDs for all point-like objects identified in the long exposures of our survey. The 100 Myr and 150 Myr isochrones constructed from the NEXTGEN and DUSTY models (see text) shifted to the distance of Blanco 1 are shown as a dashed line and dot-dashed line respectively. The mass scale (in M ⊙ ) corresponding to a cluster age of 100 Myr is given on the left hand-side of the (I,I -K) CMD. In each CMD, the solid line represents our bluest selection, i.e. for an age of 150 Myr, that maximizes the number of candidates. All the objects redder than this line on the (I, Iz) CMD have been identified as very low mass star and brown dwarf candidates. Sources with I -K bluer than the solid line in the (I,I -K) CMD have then been rejected and are shown as open circles. Furthermore, candidates with optical spectroscopy consistent with cluster membership are shown as filled triangles (open triangles if not consistent).

Fig. 8 .

 8 Fig. 8. LRIS spectra of CFHT-BL VLM candidates.

Fig. 10 .

 10 Fig. 10. Upper panel: Comparison of the Na I doublet profile of BL-29 (solid line) with two synthetic spectra of similar T e f ffor log g = 4.5 (dashed line) and 5.0 (dot-dashed line). Lower panel: Same as before but for GJ-866 which is an older field dwarf having the same spectral type than BL-29 (M6).

Fig. 11 .

 11 Fig. 11. Spatial distribution of all the cluster member candidates detected in the short exposures and selected from the 150 Myr isochrone. Contours of isodensity are drawn. The cross represents the peak of density and is taken as the cluster center whose coordinates are (0.87,-30.2) in degrees.

Fig. 14 .

 14 Fig. 14. The mass function of Blanco 1, from low mass brown dwarfs to the most massive stars, fitted by a log-normal distribution (solid line and dashed line for the 100 Myr and 150 Myr sample respectively). The large filled (100 Myr sample) and open (150 Myr sample) dots are our data points completed in the more massive star domain by Pillitteri et al. (2003) data (open squares). The Pleiades mass function (data points and log-normal fit from Moraux et al. 2003) is shown for comparison (see Section 4).

Table 1 .

 1 J2000 coordinates of the 7 CFH12K fields covered in Blanco 1.

	Field RA(2000)	DEC(2000)
		(h m s)	( • ′ ′′ )
	4 5 6 8 9 10 11	00 : 07 : 27.0 -29 : 52 : 37 00 : 04 : 19.4 -29 : 53 : 21 00 : 01 : 05.0 -29 : 58 : 21 00 : 07 : 37.0 -30 : 18 : 26 00 : 04 : 19.4 -30 : 24 : 38 00 : 01 : 05.0 -30 : 24 : 40 00 : 04 : 19.4 -30 : 50 : 29

Table 2 .

 2 Estimated number of contaminants and probable cluster members for 14 ≤ I ≤ 18.5 for both ages 100 Myr and 150 Myr.

	Assumed I-magnitude mass bin	Nb of selected estimated nb of nb of probable
	age 100 Myr 150 Myr	bin 14.0-15.5 15.5-16.5 16.5-17.5 17.5-18.5 14.0-15.5 15.5-16.5 16.5-17.5 17.5-18.5	(M ⊙ ) 0.60-0.39 0.39-0.23 0.23-0.14 0.14-0.09 0.63-0.41 0.41-0.26 0.26-0.16 0.16-0.10	objects 120 120 138 101 123 133 149 129	contaminants 52.42 ± 8.25 58.22 ± 7.02 58.88 ± 7.17 41.58 ± 6.10 56.62 ± 9.13 74.78 ± 9.43 85.94 ± 10.34 75.19 ± 10.13	members 67.58 ± 13.71 61.78 ± 13.01 79.12 ± 13.76 59.42 ± 11.76 66.38 ± 14.37 58.22 ± 14.90 63.06 ± 16.00 53.81 ± 15.22

  [A colon indicates uncertain detection due to low S/N ratio or blending.]

			SpT from	T e f f from		log g from
	BL-PC3	PC3	reference synthetic	H α	synthetic	Instrument	comments
				spectra	spectra (K) EW(Å) spectra
	16	1.26/1.27 M4.9 M5	2900	-10.0	4.8	FORS2/Keck02
	22	1.33	M5.3 M5	2900		5.0	FORS2
	24	1.48	M6.2 M6	2850			Keck (low s/n)
	25	1.33	M5.3 M5	2850		5.0	FORS2
	28	1.30	M5.1 M5	2950		5.0	FORS2	non-member ?
	29	1.41	M5.8 M6	2850		5.0	FORS2
	34	1.71	M7.4 M7	2700		4.5	FORS2	binary?
	36	1.40/1.40 M5.7 M6	2850	-5.2	5.0	FORS2/Keck
	38	1.49	M6.3 M6	2800	-6.5	4.8	Keck
	41 43	1.77/1.71 M7.6 M8 1.49 M6.3 M6	2650 2800	> -2.0 4.5 4.8	FORS2/Keck FORS2	binary?
	44	1.33/1.39 M5.5 M6	2850	-5.5	5.3	FORS2/Keck	non-member
	45	1.50/1.61 M6.6 M7	2750	-8.0	4.8	FORS2/Keck
	48	1.67	M7.2 M7	2750		5.0	FORS2
	49	1.73/1.62 M7.2 M7	2700	-2:	4.8	FORS2/Keck
	50	1.73	M7.5 M7	2700		4.8	FORS2
	51	1.65	M7.1 M7	2700	-2:	4.5	Keck (low s/n)

Table 5 .

 5 Number of very low mass star and brown dwarf candidates, estimated contaminants and probable cluster members after analysis of the follow-up observations (infrared photometry + optical spectroscopy) for both selections (100 Myr and 150 Myr).

	Assumed magnitude mass bin	number of	nb of	nb of
	age	bin (I) 18.5-19.0	(M ⊙ ) 0.090-0.075	candidates contaminants probable members 14 2±2 12±3
	100 Myr	19.0-20.0	0.075-0.056	15	1±1	14±2
		20.0-21.0	0.056-0.041	12	1±1	11±2
		21.0-22.2	0.041-0.032	29	10±4	19±5
		18.5-19.0	0.104-0.086	15	3±2	12±3
	150 Myr	19.0-20.0	0.086-0.065	19	3±2	16±3
		20.0-21.0	0.065-0.047	12	0±1	12±2
		21.0-22.2	0.047-0.037	30	10±5	20±6

Blanco 1 being at 240 pc below the galactic plane it is indeed necessary to take the scale height into account.

Note that the authors used an age of 100 Myr to estimate masses but using 150 Myr instead would not change the results in this mass range. This is not the case however at lower mass where the 100 Myr and 150 Myr mass-magnitude relationships differ.

Note that the upper mass limit is imposed by the highest mass cluster member with m ∼ 3M ⊙ (M (bol) = -0.74;Westerlund et al. 1988).

Acknowledgements. The authors thank Isabelle Baraffe who calculated the NEXTGEN and DUSTY models for the CFH12K filter set, as well as Sylvain Guieu for computing the isodensity contours.

The authors wish to recognize and acknowledge the very significant cultural role and reverence that the summit of Mauna Kea has always had within the indigenous Hawaiian community. We are most fortunate to have the opportunity to conduct observations from this mountain.


Fabrizio Flacco

Antonio Paolillo

Abderrahmane Kheddar

Residual-based contacts estimation for humanoid robots

The residual method for detecting contacts is a promising approach to allow physical interaction tasks with humanoid robots. Nevertheless, the classical formulation, as developed for fixed-base robots, cannot be directly applied to floating-base systems. This paper presents a novel formulation of the residual based on the floating-base dynamics modeling of humanoids. This leads to the definition of the internal and external residual. The first estimates the joints effort due to the external perturbation acting on the robot. The latter is an estimation of the external forces acting on the floatingbase of the robot. The potential of the method is shown proposing a simple internal residual-based reaction strategy, and a procedure for estimating the contact point that combines both the internal and external residuals.

I. INTRODUCTION

Humanoids must be embedded with capabilities of multicontact interaction with their environment or with the human. Indeed, to operate in confined spaces, it is necessary to relay on multi-contact motion and detect unexpected and undesired contacts on the robot originated from the environment or the colocated workers. Motion planning, with knowledge of the robot dynamics, the environment model, and the contacts, gives the nominal plan to accomplish the desired task, e.g. open doors [START_REF] Banerjee | Human-supervised control of the ATLAS humanoid robot for traversing doors[END_REF], climb stairs or ladders [START_REF] Vaillant | Multi-contact vertical ladder climbing with HRP-2 humanoid[END_REF], move across uneven terrain [START_REF] Fallon | Continuous humanoid locomotion over uneven terrain using stereo fusion[END_REF], etc.

Visual techniques to evaluate the appearance of the environment improved substantially, e.g. positioning and classification of objects [START_REF] Holzbach | Enhancing object recognition for humanoid robots through time-awareness[END_REF], identification of possible contact points [START_REF] Fallon | Continuous humanoid locomotion over uneven terrain using stereo fusion[END_REF], localization of the robot [START_REF] Oriolo | Humanoid odometric localization integrating kinematic, inertial and visual information[END_REF] and mapping of the environment [START_REF] Meilland | On unifying key-frame and voxelbased dense visual SLAM at large scales[END_REF]-thanks to the development of new affordable and effective sensors, such as RGB-D technology.

Meanwhile, the estimation of the interaction forces between the humanoid robot and its surrounding, as well as the identification of how these forces affect the robot internal (joint level) and external (floating-base level) motions, relay on embedded force/torque (F/T) sensors, or on whole-body tactile sensing technologies that are still not yet readily available [START_REF] Dahiya | Robotic tactile sensing-Technologies and system[END_REF] for a whole-body coverage under humanoid constraints. F/T sensors measurements -that are found only at terminal points (feet or wrists), include (interfere with) the forward links load. In such cases, the interaction forces acting on the robot can be indirectly reconstructed, e.g. by looking at the equilibrium status of the robot [START_REF] Bellaccini | Manual guidance of humanoid robots without force sensors: Preliminary experiments with NAO[END_REF].

The authors are with CNRS-University of Montpellier, Laboratoire d'Informatique de Robotique et de Microéectronique de Montpellier (LIRMM), IDH Group, 34095 Montpellier Cedex 5, France. Email: {paolillo,kheddar}@lirmm.fr.

This work is supported by the European Commission, within the H2020 COMANOID project (www.comanoid.eu), and by the ROMEO 2 project, (www.projetromeo.com), bpifrance in the framework of the Structuring Projects of Competitiveness Clusters (PSPC). The measurement of the perturbation on the internal and external state of the robot is of crucial importance for achieving physical interaction tasks. In fact, these information can be used to (i) detect if the contact point is slipping, (ii) detect unforeseen collisions and react to them, (iii) to evaluate the dynamic characteristic of a contact, (iv) etc.

In this work we estimate the joint level effects due to external forces using proprioceptive sensors only: the momentumbased residual [START_REF] De Luca | Collision detection and safe reaction with the DLR-III lightweight robot arm[END_REF], that we name classical residual. This technique, initially developed for fixed-base manipulators, is outlined in II-B. In Sec. II-C, we show that its use for humanoid robots is not proper (see Fig. 1).

Our main contributions are: a novel internal residual that estimates correctly the effect of external forces at the joint level and its computation algorithms, which outperform that of the classical residual (Sec. III); a new external residual that estimates the effect of the sum of all external forces on the motion of the floating-base (Sec. V); the techniques for merging the internal and external residual to estimate the external contact forces (Sec. V) and even the location when there is a single contact point (Sec. VI); an application using the internal residual as feedback in a joint admittance controller to comply with an external perturbation (Sec. IV).

Throughout the manuscript, we present a series of Matlab simulation results using the dynamic model of a 34 degree of freedom robot HRP-4. We use the Spatial Vector and Rigid-Body Dynamics library 1 . The simulated low-level torque control is composed by gravity compensation plus a PD control of the joint positioning error. Contacts with the ground are simulated as a simple spring-damper system.

II. BACKGROUND A. Spatial algebra

We briefly recall spatial vector algebra and notation. The reader may refer to [START_REF] Featherstone | A Beginner's Guide to 6-D Vectors (Part 1)[END_REF], [START_REF]A Beginner's Guide to 6-D Vectors (Part 2)[END_REF] and [START_REF]Rigid Body Dynamics Algorithms[END_REF] for more details.

Consider the i-th (rigid) link of the robot kinematic chain: We can express its motion and the force that acts on it with the so-called spatial motion and force vectors, that are elements of two dual vector spaces, M 6 and F 6 . With reference to Fig. 2, v ∈ M 6 and f ∈ F 6 represent the spatial velocity of the body and the spatial force acting on it. It is convenient to express the coordinate vector of a spatial vector. For example, i v ∈ R 6 and i f ∈ R 6 are the coordinate vectors representing v and f in F i , the i-th link reference frame attached to the i-th joint. The vector i v collects the angular and linear velocity around and along the axes of the frame F i . Similarly, i f is filled with the moment and the force acting on the body, and expressed with respect to F i .

The coordinate vectors can be easily transformed from a frame F A to a frame F B :

B v = B X A A v, and B f = B X * A A f . (1)
By virtue of the duality existing between the motion and force vector spaces, it is B X * A = B X -T A . The transform matrix can be expressed as follows:

B X A = B E A 0 -t AB × B E A B E A (2)
where B E A is the rotation matrix from F A to F B , and t AB × is the skew symmetric matrix associated to the position vector of the origin of F B in F A ; 0 is a 3 × 3 zero matrix 2 . Finally, the spatial inertia I of the link is a 6×6 matrix that includes information about the mass of the body, the position of the center of mass expressed in F i and the rotational inertia around the center of mass. The spatial inertia can be easily transformed from frame F A to F B :

B I = B X * A A I A X B . (3)
The spatial inertia of two bodies i and j can be merged in a composite rigid body, simply by summing their spatial inertia expressed in a common frame A I c = A I i + A I j .

B. Residual estimator

The dynamic model of a rigid fixed-base robot is

H(q)q + c(q, q) + g(q) = τ -τ µ (q, q) + τ c , (4)
where H(q) is the robot inertia matrix, c(q, q) the centrifugal and Coriolis term, and g(q) accounts for the gravity; q, q and q are the vector of the n joint generalized coordinates and its first and second time derivative, respectively 3 . τ are 1 Online: royfeatherstone.org/spatial/v2/index.html 2 We denote respectively with 1 and 0 the identity matrix and the zero matrix -dimension has to be inferred from the context. 3 In the rest of the paper we get rid of the dependence to the generalized coordinates (e.g. H(q) → H) as much as possible. the joint torques applied by the motors, τ µ (q, q) is torques due to friction and τ c represents those due to external contact forces. Considering k external forces f c , with c = 1 . . . k, their action on the joints is given by

τ c = k c=1 J T c f c . (5)
The n×m Jacobian J c can be associated either to the contact point (m = 3), at which the three-dimensional (3-D) contact force f c acts, or to the body where the six-dimensional (6-D) spatial force f c is applied (m = 6). In both cases, the Jacobian matrix takes into account the kinematic chain from the contact location to the robot base frame. Thus, the external forces f c are expressed in the base frame.

If the robot has F/T sensors to measure some of the external forces acting on the robot f c,F/T , it is possible to divide τ c in a measured part τ c and an unknown part τ c

τ c = τ c + τ c = k c=1 J T c f c,F/T + k c=1 J T c f c , (6)
with k = k + k.

The knowledge of the robot dynamic model, and the measurement of q, q and τ m = ττ µ + τ c enables the computation of the so-called classical residual vector [START_REF] De Luca | Collision detection and safe reaction with the DLR-III lightweight robot arm[END_REF]

r = K i H q - t 0 (τ m + C T q -g + r) ds (7)
where K i > 0 is the residual gain, and H q is the generalized angular momentum of the robot. The centrifugal and Coriolis term has been factorized with the form c(q, q) = C(q, q) q, that is built from Cristoffel symbols of second kind. By evaluating the time derivative of r, and exploiting the well known relation

Ḣ = C + C T , it is straightforward to show that the residual r is a filtered version of τ c (ideally, K i → ∞ ⇒ r = τ c
) and can be used for estimating the perturbation of the robot joints due to external forces. Using the gain K i it is possible to control the bandwidth of the residual. With large K i the signal will be more reactive, but also more sensitive to noise. By lowering K i , it is possible to filter the noise, but the response will be slower. A discrete time implementation of the residual is presented in [START_REF] Flacco | Online estimation of variable stiffness in flexible robot joints[END_REF]. The factorized centrifugal and Coriolis term C can be computed with a modified version of the Newton-Euler method [START_REF] De Luca | Exploiting robot redundancy in collision detection and reaction[END_REF]. Note however that the total computational complexity for obtaining C is O(n 2 ), since n instances of the Newton-Euler routine needs to be called. Therefore, this procedure could be time consuming for robots with a large number of joints.

The residual approach proved to be efficient with fixedbase robots. Initially, it was proposed for detecting and isolating actuators faults [START_REF] De Luca | Actuator failure detection and isolation using generalized momenta[END_REF], then it was used for collision detection and reaction [START_REF] De Luca | Collision detection and safe reaction with the DLR-III lightweight robot arm[END_REF], estimating the stiffness of variable stiffness actuators [START_REF] Flacco | Online estimation of variable stiffness in flexible robot joints[END_REF] and the external forces [START_REF] Magrini | Estimation of contact forces using a virtual force sensor[END_REF].

C. The use of the classical residual with floating-base robots

When the robot does not have a fixed base, the dynamic model ( 4) is still commonly used by considering six virtual joints (three prismatic and three revolute) connected to a point of the robot, called floating-base (FB). These virtual joints allow taking into account the attitude of the robot with respect to a virtual base (world) frame.

The virtual dynamic model is then described by

H v qv + C v qv + g(q v ) = τ v -τ v,µ + τ c , (8)
where subscript v indicates the virtual base version of the matrices and vectors characterizing the robot dynamics. In this model, the n + 6 generalized coordinate vector q v is composed by the 6 virtual joint coordinates q f b , and the n robot joint coordinates q, i.e. q v = q T f b q T T . The under-actuation of the floating-base is represented by having zero motor torques associated to the virtual joints τ v = 0 T τ T T . The joint torque due to external contacts τ c is split into measured τ c and unknown τ c torques as in [START_REF] Meilland | On unifying key-frame and voxelbased dense visual SLAM at large scales[END_REF].

The classical residual-based estimator of τ c (7) can be applied with the model ( 8) by adapting the matrices used to the ones of the virtual floating-base model, e.g. H v qv in place of H q. This estimated torque is similar to the fixedbase case [START_REF] Meilland | On unifying key-frame and voxelbased dense visual SLAM at large scales[END_REF]. In this case, the external force acting on the contact point follows the kinematic chain to the FB; hence, it is grounded to the virtual base frame.

Although the dynamic model ( 8) is consistent from an external point-of-view, that is, it correctly models how the robot behaves under the action of external forces, it is not correct from an internal point-of-view. In fact, the perturbations on the joints torque due to the external forces are not considered in a proper way. For instance, assume an external force acting on the right shoulder. According to this model, the torque due to the external force (5) will have nonzero values for those joints of the kinematic chain between the shoulder and the FB, and zero values for the joints of the legs. Clearly, this is not a realistic situation (see Fig. 1a). In fact, the joint torque perturbations depend on the position of the FB point, meaning that if it was positioned right on the contact point, the force would be completely grounded to the virtual joints, and no torques would be accounted for the robot joints. Thus, the estimation computed with the classical residual [START_REF] Dahiya | Robotic tactile sensing-Technologies and system[END_REF] is not able to perceive the effects of the external forces at the joint level, and cannot be used in the same way it is used for fixed-base manipulators.

To show this effect, we simulated a force applied on the right shoulder of the HRP-4. The robot task is to maintain the initial configuration. The external forces due to the ground reactions, applied to the foot are measured using F/T sensors, therefore assumed as known. Figure 3 shows how the classical residual [START_REF] Dahiya | Robotic tactile sensing-Technologies and system[END_REF], applied on the model [START_REF] Bellaccini | Manual guidance of humanoid robots without force sensors: Preliminary experiments with NAO[END_REF], is not able to estimate the real joint torque due to the external contact, that we obtained from the simulation (τ c,real ). In fact, the classical residual estimates: torques/forces on the virtual joints, which are senseless; no torques on the legs, e.g at the knees pitch joints; wrong values for the joints in the kinematic chain between the contact point and the floating-base point, e.g. at the chest yaw joint.

III. A NOVEL RESIDUAL FOR HUMANOID ROBOTS

In the previous Section, we showed that the classical residual estimator, based on the humanoid dynamic model with virtual joints for the floating-base, fails to estimate the internal joint torques due to external contacts. For overcoming this problem we have to consider the under-actuated nature of humanoid robots. To this end, we start from the floating-base dynamic model described in [START_REF]Rigid Body Dynamics Algorithms[END_REF] (p. 181), which decouples the motion of the FB with the dynamic of the robot joints

I c 0 F 0 H fb a 0 q + pc 0 C fb = 0 τ -τ µ + k c=1 1 J fb,T c 0 X * ref f c (9)
where we highlighted the contribution of the k external forces f c , expressed in the world reference frame. Let F 0 be the frame attached to the FB, and a 0 the acceleration of the FB.

H fb = H -F T (I c 0 ) -1 F
a vp 0 = -0 a g for i = 1 to n do [X J , S i ] = jcalc(jtype(i), q i ) v J = S i qi i X λ(i) = X J X T (i) v i = i X λ(i) v λ(i) + v J a vp i = i X λ(i) a vp λ(i) + v i × v J f i = I i a vp i + v i × * I i v i end for f 0 = I 0 a vp 0 + v 0 × * I 0 v i for i = n to 1 do Ĉi = S T i f i f λ(i) = f λ(i) + i X T λ(i) f i end for pc 0 = f 0
F 0 as base, F is the 6 × n matrix whose columns are the spatial forces required at the floating-base to support unit acceleration about each joint variable and I c 0 is the composite inertia of the whole robot. The 6-D vector pc 0 is the spatial bias force for the composite rigid body comprising the whole floating-base system, and C fb = Ĉ -F T (I c 0 ) -1 pc 0 is the floating generalized bias force, with Ĉ joint level bias force 4 . All spatial quantities are expressed in F 0 . For more details on how this model is derived, please refer to [START_REF]Rigid Body Dynamics Algorithms[END_REF].

The last part of the model ( 9) represents the action of external forces that perturbs both the joints state and the motion of the FB. The coordinate transformation 0 X * ref is used to transform external forces from the inertial reference frame (world) to F 0 . The action of an external forces on the joint is guided by the floating Jacobian transposed J fb,T c = J T 0,c -F T (I c 0 ) -1 , where J 0,c is the n × m Jacobian of the contact point (body) with respect to the F 0 .

As with the classical modeling, we can consider the joint torque contribution of known external forces as follows

τ c = k c=1 J fb,T c 0 X * ref f c,F/T . (10)
To estimate the joint torques due to unknown external contacts, we build a new residual using the second line in (9)

r int = K i,int H fb q - t 0 (τ m + Ḣfb q -C fb + r int ) ds (11
) with K i,int > 0. This residual is also a first order filter of the joint torques due to external forces

r int ≈ τ c = τ c,real = k c=1 J fb,T c 0 X * ref f c . (12)
Since this residual estimates the joint level effects of an external force, we call it internal (floating-base) residual.

Algorithm 2 Computing H, F , I c 0 , Ḣ, Ḟ and İc

0 H = 0 , Ḣ = 0 for i = 1 to n do [X J , S i ] = jcalc(jtype(i), q i ) v J = S i qi i X λ(i) = X J X T (i) i Ẋλ(i) = -v J × i X λ(i) I c i = I i , İc i = 0 end for I c 0 = I 0 , İc 0 = 0 for i = n to 1 do I c λ(i) = I c λ(i) + i X T λ(i) I c i i X λ(i) İc λ(i) = İc λ(i) + i X T λ(i) İc i i X λ(i) + i ẊT λ(i) I c i i X λ(i) + i X T λ(i) I c i i
Ẋλ(i)

F i = I c i S i , Ḟ i = İc i S i H ii = S T i F i , Ḣii = S T i Ḟ i j = i while λ(j) = 0 do Ḟ i = j X T λ(j) Ḟ i + j ẊT λ(j) F i , F i = j X T λ(j) F i j = λ(j) H ij = S T i F i , Ḣij = S T i Ḟ i H ji = H ij , Ḣji = Ḣij end while Ḟ i = j X T 0 Ḟ i + j ẊT 0 F i , F i = j X T 0 F i end for
The time derivative of the floating-base inertia matrix Ḣfb can be obtained analytically as

Ḣfb = Ḣ - Ḟ T (I c 0 ) -1 F -F T (I c 0 ) -1 Ḟ - F T -(I c 0 ) -1 İc 0 (I c 0 ) -1 F . ( 13
)
Moreover, Ĉ and pc 0 can be computed using the Algorithm 1, and H, F , I c 0 and their time derivative using Algorithm 2, which is an extension of the method proposed in [START_REF]Rigid Body Dynamics Algorithms[END_REF] (p. 183). In these algorithms, 0 a g is the gravity acceleration expressed in frame 0, λ(i) is the body index of the parent of body i (0 ≤ λ(i) < i). jcalc(•) is a function that according to the type of joint jtype(•) (e.g. revolute or prismatic) returns the joint coordinate transform X J (for revolute joint it is a pure rotation of q i around the z-axis), and the joint motion subspace vector S i that describe the motion axis of the joint (for revolute joint S i = [ 0 0 1 0 0 0 ] T ).

Summing up, the information needed to compute the internal residual [START_REF]A Beginner's Guide to 6-D Vectors (Part 2)[END_REF] are: joint generalized position q and velocity q; FB rotation w.r.t the virtual base frame (for computing 0 a g ) and FB velocity v 0 ; joint motor torque τ , joint friction torque τ µ and joint torques due to known external forces τ c , for computing τ m .

Thanks to Algorithms 1 and 2, the internal residual (11) has a lower computational load w.r.t the classical residual [START_REF] Dahiya | Robotic tactile sensing-Technologies and system[END_REF]. In fact, the algorithm complexity is O(n). For instance, the mean computational time to obtain the classical residual for HRP-4 in the previous simulative example is 0.0795 s, while it is 0.005 s for the internal residual, almost 16 times faster.

We simulated the same sinusoidal external force applied on the right shoulder of the HRP-4 presented in the previous section. Fig. 4 shows the estimated joint torques, using the presented internal residual. The estimated torque at the knees pitch and at the chest yaw joints are now correct, as well as all the other joints (that we omit for lack of space). Note that, compared to Fig. 3, the virtual joints are not presented as their values are identically zero.

IV. AN APPLICATION OF THE INTERNAL RESIDUAL

The estimation of the joint torques due to external perturbation, as computed by the internal residual ( 11) can be used in many ways, for example in the framework of physical human-robot interaction. In order to show the potential and the effectiveness of the proposed method, we propose a possible application where HRP-4 is requested to move on the direction of the external forces without knowing neither their number, the contact points nor the force vector. For this purpose, we implemented a simple joint admittance control that generates the desired joint configuration as

qd = K a τ c,real , q d = q + qd T (14)
where K a is a diagonal matrix of gains and T is the simulation sampling time.

In this example, the robot stands at its position for 15 s; at t = 5 s, we applied two downwards forces on the robot shoulders, with a polynomial profile increasing till t = 8.5 s and then decreasing till t = 12 s. The maximum value of the forces was set to 5 N. As if a human push down the robot expecting that the robot bends for a seating configuration. We compared the results obtained by estimating τ c,real with r, as computed with [START_REF] Dahiya | Robotic tactile sensing-Technologies and system[END_REF], and r int , given by [START_REF]A Beginner's Guide to 6-D Vectors (Part 2)[END_REF].

Fig. 5 shows snapshots of the simulation performing the control with the classic residual, while Fig. 6 shows snapshots of the simulation obtained using the internal residual. The entire motion can be seen in the accompanying video. Reaction strategy to external perturbations using the internal residual: simulation results. The joints effort due to the external forces (red lines) is felt at all the involved joints of the kinematic chain (gradually highlighted in red in proportion to the effort). The admittance control implements a whole-body motion, ensuring an effective reaction strategy.

As explained in Sec. II-B, the classical residual does not allow a whole-body reaction strategy, preventing a proper robot behavior when trying to comply with external forces. On the other side, the simulation that makes use of the internal residual confirms the expectation of Sec. III. It allows to perceive the effort due to the external force at any interested joints of the robot. As a result, the robot can comply with the perturbation with more degrees of freedom and ensure a more effective reaction strategy.

V. EXTERNAL FORCE ESTIMATION

Now, we consider one of the most important use of the residual: the estimation of the force applied on a contact point. In [START_REF] De Luca | Integrated control for pHRI: Collision avoidance, detection, reaction and collaboration[END_REF] it is described how the residual r provides an estimation of the k unmeasured external forces f c from the knowledge of the contact points location:

   f 1 . . . f k    = J T 1 . . . J T k # r . (15)
Successively, how this estimation can be exploited in the control is presented in [START_REF] Magrini | Control of generalized contact motion and force in physical human-robot interaction[END_REF]. The dynamic model of a floating- base robot with virtual joints [START_REF] Bellaccini | Manual guidance of humanoid robots without force sensors: Preliminary experiments with NAO[END_REF] is not proper from an internal viewpoint (Sec. II-C), but the external behavior of the robot is correctly modeled. This is similarly observed in the estimation of the external forces. In fact, while the residual estimator [START_REF] Dahiya | Robotic tactile sensing-Technologies and system[END_REF] is not able to represent properly the internal torques, the estimated external forces using (15) are correct. An example of the use of the classical residual for estimating external forces on a humanoid robot is presented in [START_REF] Manuelli | Localizing external contact using proprioceptive sensors: The contact particle filter[END_REF], where a particle filter takes advantage of the assumption of the knowledge of the friction cone of the contact to estimate the contact point and the applied force. Yet, we show that the use of the floating-base dynamic model ( 9) is more flexible because it allows to decouple the internal and external robot behavior due to external forces. Thus, these information can be used separately or combined for estimating the external forces applied to the robot.

Considering the upper row of the decoupled dynamic model ( 9), we have a relation between the robot motion and the total external forces applied to the robot. Taking into account this 6-D information, we can build a novel residual that we call external (floating-base) residual:

r ext = K i,ext I c 0 v 0 + F q - t 0 ( İc 0 v 0 + Ḟ q -pc 0 + f sum + r ext ) ds (16)
where K i,a > 0 is a residual gain matrix, and

f sum = k c=1 0 X * ref f c,F/T
is the sum of all measured forces. Note that all matrices to compute the external residual ( 16) are obtained with Algorithms 1 and 2. Moreover, this residual is much more reliable, since it does not need the measure of the motors torque and the estimation of joints friction behavior, which are usually more difficult to obtain. As previously, it is simple to show that r ext is a first order filter of the sum of all unknown external forces applied on the robot (expressed

r ext ≈ k c=1 0 X * ref f c . (17)
Using the external residual on the example initially proposed in Sec. II-C, where a single sinusoidal force is applied to the right shoulder, we obtained the estimation of the force (expressed in the floating-base frame) as shown in Fig. 7.

With this new tool, if one need to know just the resultant of all external forces applied to the robot, e.g. for push recovering, you can use it directly. If we assume to have a single unmeasured 3-D force f 1 applied to one contact point, whose position p 1 is known/detected, its estimation is given directly by the last three components of r ext opportunely rotated. In fact, it is

r ext = r ext,m r ext,f ≈ 0 X * ref p 1 × f 1 f 1 (18)
thus,

f 1 ≈ ref E 0 r ext,f . (19)
Figure 8 shows the reconstructed external force applied to the robot shoulders in the previous example. When two or more unmeasured 3-D forces are applied to known robot points, the external residual has to be combined with the internal residual for estimating the forces. The most straightforward method is to use pseudo-inversion

   f 1 . . . f k    = 1 . . . 1 J fb,T 1 . . . J fb,T k # ref X * 0 r ext r int .
(20) When the two residual are combined, we should be aware that they use different sensor information, affected by different source of noise. The quality of the internal residual is mainly disturbed by the noise on the measure of the motor torque and by the estimated values of the joint friction, while the external sensor depends on the estimation of the orientation and the spatial velocity of the floating-base. Thus, we suggest to use different values for K i,int and K i,ext . The intersections of l 1 with the robot links localize the possible contact points (a). These points are used to compute the Jacobian and, along with f 1 , the hypothetic joint effort due to the external perturbation as it was acting on those points (b), (c). The real contact point corresponds to the one producing the joint efforts that matches with the internal residual.

VI. CONTACT POINTS DETECTION AND FUTURE WORKS

One of the more challenging issue is finding the contact position where the external forces are applied, knowing just the number of contact points, or even estimating it. State-ofthe-art methods are based on using external sensors, such as tactile skin [START_REF] Cirillo | A conformable force/tactile skin for physical human-robot interaction[END_REF] or depth camera [START_REF] Magrini | Estimation of contact forces using a virtual force sensor[END_REF].

We explore the possibility to combine the internal and external residual to estimate the position of unforeseen contacts, without using any extra sensors. Meanwhile, we present here how the new two residual can be combined for estimating both the contact position and the applied force, under the assumption of one single unforeseen contact. This situation is the same proposed in the previous Section, with the external residual estimates [START_REF] Magrini | Control of generalized contact motion and force in physical human-robot interaction[END_REF], but in this case the contact position p 1 is considered unknown. Note that eq. [START_REF] Manuelli | Localizing external contact using proprioceptive sensors: The contact particle filter[END_REF] gives an estimation of the applied force, even if the contact position is unknown. At this point, we consider

ref r ext = ref r ext,m ref r ext,f = 0 X T ref r ext , (21)
where, trivially ref r ext,f = ref E 0 r ext,f ≈ f 1 . From the upper row of (18), we know that the contact point lies on the line

l 1 = ref r ext,f × ref r ext,m ref r ext,f 2 + l ref r ext,f , l ∈ R. ( 22
)
This line is oriented as f 1 and passes through p 1 .

Using a CAD representation of the robot, or a simpler representation with primitive shapes, we can individuate all possible contact points on the robot surface that intersect line l 1 . Thus, for all of them we compute the joint torque perturbation we would have observed if the force ref r ext,f had applied on that point. By comparing those joint torques with the internal residual r int , we can evaluate which among the possible contact points is compatible with the estimated internal behavior, thus we can discriminate the correct contact position. This procedure is illustrated in Fig. 9.

For the experimental validation of the presented method on the real HRP-4, we plan to use the two residuals to refine the parameters of the robot dynamic model. In fact, when no unknown contact are present, the residuals give an estimation of the discrepancy between the robot dynamics and the used model. From another side, we plan to improve the estimation of the position and velocity of the floating-base by fusing IMU and visual information.

Fig. 1 .

 1 Fig. 1. Residual computation for humanoid robots. (a) Using a classical dynamic model, the effort due to an external force is unload to the ground through the 6-D virtual joints composing the floating-base. (b) The use of a floating-base dynamic model, more appropriate to the nature of humanoids, allows to correctly compute the joint effort due to the external perturbation.

Fig. 2 .

 2 Fig.2. Schematic representation of the i-th robot link and variables v (the spatial velocity), f (spatial force), and I (spatial inertia matrix), expressed w.r.t. the reference frame F i , placed at the i-th joint.

Fig. 3 .

 3 Fig.3. Estimated joint force/torque due to the external contact (solid) and real values (dashed). The six virtual joints (top), the right and left knees pitch (middle up), the chest yaw (middle down) and norm of the error in estimating the real toque due to the unknown external force (bottom).

Fig. 4 .

 4 Fig. 4. Estimated joint torque due to the external force (solid) and real values (dashed). The right and left knees pitch (top) and the chest yaw (middle) and norm of the error in estimating the real toque (bottom).

t = 4 s t = 6 s t = 11 sFig. 5 .

 115 Fig.5. Reaction strategy to external perturbations using the classic residual: simulation results. The joints effort due to the external forces (red lines) is felt only at the pitch joint of the chest (gradually highlighted in red in proportion to the effort). As a result, the admittance control acts only on this joints, preventing an effective reaction strategy.

Fig. 7 .

 7 Fig. 7. Estimated total torque (top) and force (middle) due to the external contact (solid) and real values (dashed). The estimation error norm (bottom)

1 Fig. 8 .

 18 Fig. 8. External force due to a single external contact (top): estimated (solid) and real values (dashed). Norm of the estimation error (bottom)

Fig. 9 .

 9 Fig.9. Estimation of the external force and the contact point. The external residual estimates the contact force f 1 and the line l 1 where it lies on. The intersections of l 1 with the robot links localize the possible contact points (a). These points are used to compute the Jacobian and, along with f 1 , the hypothetic joint effort due to the external perturbation as it was acting on those points (b), (c). The real contact point corresponds to the one producing the joint efforts that matches with the internal residual.

The hat superscription indicates that the bias action of external forces is not included, since we have expressed it separately.

ACKNOWLEDGEMENT

The authors would like to thank Prof. Roy Featherstone for his advice on the realization of the Algorithm 2.