A&A 486, 283-291 (2008)
DOI: 10.1051/0004-6361:200809672
N. R. Deacon^{1} - G. Nelemans^{1} - N. C. Hambly^{2,}^{}
1 - Department of Astrophysics, Faculty of Science, Radboud University Nijmegen, PO Box 9010, 6500 GL Nijmegen, The Netherlands
2 -
SUPA, Institute for Astronomy, University of
Edinburgh, Royal Observatory, Blackford Hill, Edinburgh EH9 3HJ, UK
Received 28 February 2008 / Accepted 14 April 2008
Abstract
The stellar mass function is one of the fundamental distributions of stellar astrophysics. Its form at masses similar to the Sun was found by Salpeter (1955, ApJ, 121, 161) to be a power-law
with a slope of
.
Since then the mass function in the field, in stellar clusters and in other galaxies has been studied to identify variation due to environment and mass range. Here we use results from previous papers in the SIPS series to constrain the mass function of low mass stars (
). We use simulations of the local low mass stellar population to model the results of the SIPS-II survey. We then vary the input parameters of these simulations (the exponent of the mass function
and a stellar birthrate parameter )
and compare the simulated survey results with those from the actual survey. After a correction for binarity and taking into account potential errors in our model we find that
for the quoted mass range.
Key words: stars: low mass, brown dwarfs - stars: luminosity function, mass function - astrometry
Inferring the form of the mass function from the luminosity function is relatively easy in single age populations such as open clusters. This is simply a conversion using a mass-luminosity relation. However when dealing with the field population - where the ages of stars are spread over the whole history of the Galaxy - this becomes more difficult. Hence we must also consider the birthrate . Schmidt (1959) predicted a declining birthrate based on the assumption that it is related to the density of interstellar gas. However, as noted by Miller & Scalo (1979), because his work came before the widespread acceptance of a hot Big Bang (and hence primordial nucleosynthesis) he used an initial helium abundance of zero - clearly this is not correct. Miller & Scalo (1979) went on to use a continuity constraint, asserting that the mass function should be smooth to derive the birthrate. They rejected Schmidt's declining form instead preferring a roughly constant birthrate. Some studies claim the birthrate is not smoothly varying, for example Rocha-Pinto et al. (2000), who used the chromospheric activity of stars to derive their ages. From this age distribution they derived a birthrate that shows a series of three or four bursts of star formation in the last 10 Gyr. They note some of the bursts may be linked to close encounters with the Large Magellenic Cloud.
Finally any study using only surveys must include a correction for unresolved binarity which will both alter the colours of objects and hide cool, dim companions. The binary fraction of higher mass stars is fairly high. G dwarfs for example have a binary frequency of % (Duquennoy & Mayor 1991). For low mass stars the fraction may be lower. Fisher & Marcy (1992) found % for M 0-M 4 dwarfs and Maxted & Jeffries (2005) found a binary fraction of 32-45% for stars below about 0.15 .
In this paper we model the SIPS-II survey (Deacon & Hambly 2007). This is a proper motion survey combining SuperCOSMOS scans of United Kingdom Schmidt Telescope (UKST) I plates (Hambly et al. 2001) and data from the 2MASS survey (Skrutskie et al. 2006). This produced a sample of roughly 7000 low mass stars with proper motions between a tenth and half an arcsecond per year. In this paper, we outline how we modelled this survey and how these models have been used to constrain the underlying mass function and birthrate.
The method used here is based on the simulations used for Deacon & Hambly (2006). These are similar in principle to those of Burgasser (2004). A short overview of the techniques used in Deacon & Hambly (2006) follows. Individual objects are assigned masses and ages drawn from a mass function (with masses between 0.5 and 0.03 solar masses) and a birthrate (with ages up to 10 Gyr). Models are then used to derive the absolute magnitudes, space positions and space velocities of each object. These can be used to calculate observable parameters such as proper motion, sky position and apparent magnitude. These parameters are then passed through a survey selection mechanism to yield the results of the simulated survey. By varying the input parameters of the simulated survey the results will also vary. These can then be compared with the actual results to constrain those input parameters. Below we outline the calibrations to produce the photometric and astrometric parameters of each simulated object and the selection mechanism imposed on them. Additionally how the input parameters are constrained is described.
Figure 1: The band bolometric correction vs. effective temperature relation. The dots represent data points from Golimowski et al. (2004) and the stars symbols plot data from Berriman & Reid (1987). Halo objects and objects in binaries have been excluded from the sample. The line is a polynomial fit to the data. Note the scatter on these data points, this is probably due to measurement errors both in the determination of apparent magnitudes but more importantly to the errors in the trigonometric parallaxes of the objects. In the Baraffe et al. (1998) models a 0.5 solar mass star would have a temperature of 3650 K and a 0.1 solar mass star a temperature of around 2900 K. | |
Open with DEXTER |
Figure 2: The relation between absolute magnitude and the absolute magnitudes in other passbands. The dots are taken from Reid et al. (2002)'s 8 pc sample, the squares are from Dahn et al. (2002) and the stars are a subset of Cruz et al. (2007)'s sample. In each case there is a polynomial fit to the data plotted. Note the larger scatter (much of it intrinsic) on both the R and the I data. | |
Open with DEXTER |
Figure 3: A plot showing the points of binary data from Delfosse et al. (2000). Along with these is our model for a 5 Gyr old population (solid line), a fit to the Delfosse data points (dashed line) and the raw Baraffe (1998) and Chabrier (2002) models (dot-dash line) also for a 5 Gyr old population. There is no major discrepancy. Note the fit to the Delfosse points is not extrapolated beyond 0.1 as there is no data below this point. | |
Open with DEXTER |
Once the masses and ages for a set of objects have been assigned their photometric characteristics can be calculated. Effective temperatures and bolometric magnitudes for these objects are found from evolutionary models (from Baraffe et al. 2003, for objects below ; and from Baraffe et al. 1998, for objects in the range ). In order to ground our simulations in observational data as well as theoretical models an empirical photometric model was used. An effective temperature vs. bolometric correction relation was derived from the data contained in Golimowski et al. (2004) and Berriman & Reid (1987) with the photometry converted to the 2MASS system using the conversions in Carpenter (2001)^{}. The polynomial fit for the effective temperature - bolometric correction relation is shown in Fig. 1. The absolute magnitudes could then be calculated from the objects' bolometric magnitudes. To test that this method of obtaining was not flawed we plotted M_{K} vs. mass for our model at 5 Gyr and the models of Baraffe (1998) and Chabrier (2003) against the empirically measured masses and M_{K}s from Delfosse et al. (2000). Clearly there is no significant offset. We had decided against using the Delfosse data in our models as it would not allow us to take into account the luminosity evolution of objects. These band absolute magnitudes are then converted to R, I, J and H magnitudes using a relation between these magnitudes and the magnitude. This was found by taking data from Reid et al. (2002)'s eight parsec sample, Cruz et al. (2007) and Dahn et al. (2002) and fitting polynomials to these data. The to other passband fits are shown in Fig. 2 and their coefficients are given in Table 1. Finally each object was given a small Gaussian scatter in R (0.25 mag) and I (0.12 mag) magnitude to simulate the scatter in these passbands on the HR diagram. These values were arrived at by allowing the values of to be free parameters and selecting those which produced the best model-data fit. The offset values are also comparable with the scatter around our R and I to relations once measurement errors in the sample used to produce the photometric model are taken into account. By doing this we can take into account the effect the intrinsic scatter has on our magnitude limited sample. However we do not model a scatter that is metallicity dependent and hence we may miss some age/metallicity dependence.
Table 1: The polynomial fit used in the model.
There are two distinct parts to the astrometric simulation, the velocity simulation and the space positions simulations. The velocity simulations use a simple thin disk model (we see no prominent seperate thick disk population in our sample). These had velocity dispersions (Seabroke & Gilmore 2007) which are typical values for a thin disk model. These are assigned randomly and are not related to space positions. A solar reflex velocity was also added with values (U,V,W)=(10.5,18.5,7.3) (Makarov & Murphy 2007). The positons are also based on a simple thin disk model with a scale height of 300 pc. In the small region of the Galaxy our survey covers (d<300 pc) we neglect the scale length of the disk which is typically 3500 pc (de Vaucouleurs & Pence 1978). Hence we assign random x and y positions in the Galactic plane drawn from flat distributions and a z position out of the plane drawn from a declining exponential with a scale height of 300 pc^{}.
Figure 4: The completeness estimates for the five passbands in the survey. The histogram represents the number of objects in the actual survey and the straight line is a fit to the trend in the area where the survey is complete. The drop-off from this line is the incompleteness estimate used in the survey. The two dashed lines represent the I and J band limits of the SIPS survey. | |
Open with DEXTER |
The astrometric selection fell into two categories: selection by proper motion, and selection based on sky position. The proper motion selection seems simple enough, we just selected on the proper motion of the object after a positional offset had been converted into a proper motion offset using the epoch difference for the particular sky area the object fell into. For selection based on sky position we included the survey area (roughly 20 000 sq. deg. based on the area covered by the UKST plates and a Galactic latitude cut) but excluded crowded areas. The crowding estimate (including the area obscured by bright stars) used the calculations from Deacon et al. (2005). This calculates a probability an object in a particular area will be obscured by a bright star or by crowding. This probability is then used to include or exclude simulated stars from the selected sample.
While simulating the results of surveys is an interesting task for predicting the potential results the real goal is to establish what surveys actually tell us. As stated earlier, the birthrate and the mass function will affect the results of a survey. By varying these parameters and then comparing them to data we can constrain these underlying distributions.
Suppose we take a mass function with a particular value of
and we define a birthrate parameter
such that
(1) |
(2) |
Figure 5: A simulation (dashed line) using the best fit parameters is shown along with the actual observed data (solid line). | |
Open with DEXTER |
Figure 6: The colours of stars in the sample ( left hand panels) and of simulated objects in the best fit model ( right hand panels). The comparison appears good with the exception of a slightly higher scatter in H-K colours in the models. | |
Open with DEXTER |
Figure 7: Histograms for the data ( left hand side) and simulation results ( right hand side). Note the slight difference in the brightness distributions and the higher spread in colours in the simulations. | |
Open with DEXTER |
The grid of simulations was produced as detailed above with the survey results being represented by an I-J histogram. This was originally chosen as I-J colour increases into the T dwarf regime unlike other IR colours. The grid of histograms was then compared to the observed histogram to produce a datacube of values and hence probabilities. This allowed us to produce constraints on the parameters , and . The derived value for is and the value for the space density of stars in the range 0.09-0.1 () is found to be pc^{-3}. The calculated accuracy on the value of is smaller than our grid's resolution (step size). All we can state is that the one limit for lies somewhere between our value and the step size (0.08) away from it. Hence we conservatively choose the value of the error in to be 0.08. The probability distribution for is too noisy to produce any sensible constraint. This is most likely due to the low number of L dwarfs in the SIPS-II sample (fourteen in total) as the characteristics of the mainly hydrogen burning M dwarfs which dominate this sample do not change rapidly over time, unlike and L dwarfs below the substellar limit. The results for and assume no errors in the model beyond simple counting errors and do not include binarity. Hence they should be assumed to apply to the system mass function (a mass function calculated from a luminosity function with no binary correction applies, see Chabrier 2002). The best fit has a per degree of freedom of 1.64 and is shown in Fig. 5. Additionally various colours of the individual objects in the SIPS-II sample and in the best fit simulation are plotted in Fig. 6 while histograms for various colours, apparent magnitudes and the proper motion are shown in Fig. 7.
Clearly we cannot assume that our model is perfect. The most likely source of potential errors is the photometric model. To attempt to quantify this we examined the effects of a simple offset in the I-J colour. We did this for two reasons, firstly I-J is the colour used in the histogram for comparing results to surveys and secondly the I magnitude fit has much more scatter than the J or H fits (see Fig. 2). Note we do not examine the R band fit as, while it is as noisy as the I band fit, it is not used as a cut for all stars (as some will not have a paired R magnitude due to photometric incompleteness and stars moving beyond the SuperCOSMOS pairing radius of 6 arcsec).
A series of simulations with the best fit value of and a range of I-J offset between -0.1 and 0.1 mag were produced. These were then compared with the simulated grid and the probability distributions in produced by this process were added. The standard deviation of this combined distribution was then measured to determine the error due to the offsets. This was found to be 0.21. Hence we find that our value for ( ) for the system mass function to be . Additionally our density normalisation parameter has an error of 0.0008 from potential photometric offsets.
Figure 8: The effects on the mass function exponent and the normalisation factor from the inclusion of different binary fractions ( ) in the simulations. Note that values of are due to noise in the simulations. | |
Open with DEXTER |
Figure 9: A plot showing the number of objects of different masses passing through the survey selection mechanism in the best fit simulations. The dashed line represents the best fit model after binary correction and the dotted line the best fit model with no binary correction. The solid line shows the best fit mass spectrum (a power law with an index ). The change in shape between the mass spectrum and the range of detected masses indicates the incompleteness. | |
Open with DEXTER |
As stated above the mass function is strongly affected by the inclusion of unresolved binary systems. In order to remove this effect and to find the form of the individual object mass function we undertook a series of simulations. In each a proportion of the sample (the binary fraction
)
was assigned an unresolved binary companion. These had their masses drawn from the same mass function as the primary objects. In each passband the luminosities of the objects were added and then converted to magnitudes. This allows us to include both the brightening due to the luminosity of the unresolved companion and the effect on the colour of the unresolved object. In each of these simulations both primary and secondary objects counted towards our density normalisation factor .
After initial test runs a value of
was selected to give a result after the effects of binarity to the measured value of
for the system mass function.
Figure 10: A plot showing the different values of found by this and other studies. The vertical error bars represent errors in while the horizontal error bars represent the range of masses over which the value of is valid. | |
Open with DEXTER |
Once these simulations were carried out they were compared to the grid of simulations with zero binarity. The value of alpha measured then had the input value of subtracted to yield the change in which we shall call . The correction factor for the normalisation factor was also calculated for each different value of . Figure 8 shows the effect on both parameters for a range of values of . Note that we assume a value of that does not change with distance. Clearly this is not correct as nearby objects will be more likely to be resolved. As the lower mass objects in our sample tend to be nearer by this may introduce a bias into the measurement.
In order to make the correction to yield the individual object mass function we take the binarity estimate of Maxted & Jeffries (2005). As our survey uses 2MASS data and photographic plates which have fairly low resolution we shall assume that all binaries are unresolved. The Maxted & Jeffries (2005) estimate of 32-45% varies due to different underlying distributions such as the distribution of separations. For our correction we shall take the middle of this range to be our binary estimate and use the extremities of it as the one error. Hence we use . Interpolating between points we get a value for of . The uncertainty in the value of adds an additional error of 0.09 to the determination of . Hence . Applying this correction to the system value of ( ) we find the value of the exponent for the individual object mass function ( ) to be .
Using the assumed value of above we can also find the density of individual objects with masses between 0.1 and 0.09 solar masses ( ). We find that . This gives pc^{-3}.
Now we have an estimate for the value of we must estimate over what range of masses it is valid. As we have very few L dwarfs we cannot claim we are probing the substellar regime. Hence a lower limit of seems sensible. As for an upper limit, excluding photometric errors the I-J histogram cuts out at around 1.5. This equates to . Using the mass-luminosity relations from Delfosse et al. (2000) this gives us an upper limit of . Additionally we can use our simulations to estimate the region of mass function applicability. Examining the range of masses passing through the simulation (see Fig. 9) we find a sharp drop below the hydrogen burning limit, a peak at around 0.1 solar masses followed by a steady decline towards higher masses.
In order to compare our calculated value of with those of other studies an plot was produced. This is shown in Fig. 10. This value of differs just outside the quoted errors from those of Kroupa (2001) at the high end and middle of our mass range. However it agrees with Allen et al. (2005) at the low end and is in good agreement with Zheng et al. (2001). In the middle of the mass range the gradient of the Chabrier (2005) lognormal form differs from ours just outside the error bounds. However assuming that the errors on his parameters are of similar magnitude to those in Chabrier (2001) the two determinations agree within one sigma. We differ significantly from the studies of Reid et al. (2002) and Martini & Osmer, these studies however cover a much larger (and mostly higher) mass range than our sample.
Burgasser (2004) uses the studies of Reid et al. (2002) and Chabrier (2001) to derive a value for the number density in the region 0.1 to 0.09 solar masses. He estimates this to be pc^{-3}. Our estimate of pc^{-3} agrees with this within one sigma. The mass function calculated by Chabrier (2005) gives a value for the number density in this region of 0.0036. While no errors are quoted on the parameters of this mass function it is clear that this value agrees well with that calculated in this work.
We have used simulations of the low mass star population to attempt to constrain the birthrate and mass function. Unfortunately no clear constraint could be set on the birthrate. However after a correction for binarity and taking into account the potential errors in our model we found a value of for the exponent of the mass function power law (). Additionally we find a constraint on the number density of stars with masses in the region 0.1-0.09 solar masses of . Both these results are consistent with some studies in the field.
The obvious next step for such work is to extend it to cover other surveys for low mass stars and brown dwarfs such as the UKIDSS Large Area Survey (Lawrence et al. 2007). That survey has substantially more accurate photometry compared to the UKST I data used in this study. Additionally the scope of the simulations could be extended to include different stellar populations such as low metallicity halo objects, allowing a more accurate model of the local low mass population.
Acknowledgements
N.R.D. is supported by NWO-VIDI grant 639.042.201 to P. J. Groot, G.N. is supported by NWO-VENI grant 639.041.405; The authors would like to thank Sue Tritton and Mike Read for their help in selecting plates, Harvey MacGillivray and Eve Thomson for their scanning of the material on SuperCOSMOS and David Bacon, Andy Taylor, Thomas Kitching, Hugh Jones and John Cooke for their helpful discussions. This publication makes use of data products from the Two Micron All Sky Survey, which is a joint project of the University of Massachusetts and the Infrared Processing and Analysis Center/California Institute of Technology, funded by the National Aeronautics and Space Administration and the National Science Foundation. SuperCOSMOS was funded by a grant from the UK Particle Physics and Astronomy Research Council. This publication makes use of the SLALIB positional astronomy library (Wallace 1998).