© ESO, 2017

1. Introduction

Young stars are subjected to several phenomena that affect their overall evolution. Accretion, ejection, and angular momentum evolution are connected in a complex way (see e.g. Bouvier et al. 2014). For solar-type stars, the disk seems to play an important role in preventing the star from spinning up to its break-up speed, but it has a short duration. Even so, the disk influence on the rotational properties of young stellar clusters is remarkable. Rebull et al. (2004) and Herbst & Mundt (2005) investigated with observational data and simple models the influence of accretion disks on the evolution of rotational rates in young clusters. In Vasconcelos & Bouvier (2015, hereafter Paper I), Monte Carlo simulations were used to show that indeed the rotational properties of clusters younger than 20 Myr can be explained by the co-existence of stars with disks prevented from spinning up and diskless stars that are being rotationally accelerated while conserving angular momentum. The double-peaked period distributions seen in the Orion nebula cluster (ONC), in NGC 2264 (Cieza & Baliber 2007), and in h Per (Moraux et al. 2013) show both populations present in the peaks, but there is a predominance of stars with disks in the long-period peak. This was seen not only in the period distributions, but also when we examine disk fractions and mass accretion rates as a function of the period. In Paper I we also explain the existence of slowly rotating diskless stars as stars that have long-lasting disks and that had no time to spin up. Another finding of Paper I is that the specific angular momentum evolution is not very different for disk and diskless stars in a young cluster. The reason for this is that the disk lifetime is not the same for all stars and thus there is a sequential release from the disk locking condition for a number of stars with a disk. Individual diskless stars will start to spin up and conserving angular momentum, but the average specific angular momentum of the diskless sample as a whole will decrease.

However, is this rotational behaviour the same for very low-mass (VLM) stars or brown dwarfs (BDs)? Several works show that these low-mass objects rotate faster than their solar mass counterparts (Scholz et al. 2015). Also, there is a debate in the literature about the disk fraction at this mass range, with some studies pointing to a higher disk fraction (Luhman et al. 2008; Luhman & Mamajek 2012 and Downes et al. 2015 – see also Table 2). There is also evidence that disk locking is not as efficient (Lamm et al. 2005).

In this work, we investigate the main variables that can influence the spin rate evolution of a cluster of very low-mass stars and brown dwarfs using Monte Carlo simulations, and compare the results to observations available in the literature. We consider the evolution of disk and diskless stars and brown dwarfs from 1 Myr to 12.1 Myr. In Sect. 2, we explain the main assumptions of the simulations. In Sect. 3, we present and discuss the different models considered in this work. In Sect. 4, we draw our conclusions.

2. Monte Carlo simulations

We follow the same set-up as described in detail in Paper I. However, we consider another mass range, from 0.01 M_⊙ to 0.4 M_⊙, encompassing the BD (0.01 M_⊙ ≤ M_∗ ≤ 0.07 M_⊙) and the VLM (0.07 M_⊙<M_∗ ≤ 0.4 M_⊙) regimes, separated into bins of 0.01 M_⊙ from 0.01 M_⊙ to 0.1 M_⊙, and into bins of 0.1 M_⊙ for higher masses. We simulate more than 200 000 stars and for each of them we assign a mass M_∗, a mass accretion rate Ṁ_acc, and a rotational period P and we evolve the system from 1 Myr to 12.1 Myr. The number of stars per mass bin is calculated from the canonical initial mass function (IMF) by Kroupa et al. (2013). The mass accretion rate values are randomly chosen from 13 different log-normal distributions, one per each mass bin, all with σ = 0.8 dex and with mean values given by (1)where t₀ = 1.0 Myr is the initial age of the simulations. The dispersion of the mass accretion rate was chosen in order to reproduce the observed spread of more than 2 orders of magnitude in Ṁ_acc seen, for example, in NGC 2264 by Venuti et al. (2014) and in the ONC by Manara et al. (2012). The mass dependency of the mean mass accretion rate follows from recent results obtained for Lupus and ρ – Ophiucus (Alcalá et al. 2014; Manara et al. 2015) for stars less massive than 0.5 M_⊙. In Paper I we have used from the results of Venuti et al. (2014) for NGC 2264. However, as stated there, our results – namely the rotational distributions, the disk fraction as a function of time, and even the behaviour of the mass accretion rate as a function of the rotational period – do not depend on the choice of this exponent. In Table 1 we show the different values of ⟨ Ṁ_acc(t₀,M_∗) ⟩ and also the maximum mass accretion rate value, which is between 4σ–5σ, at each mass bin.

Since we want to disentangle the contributions of disk and diskless stars on the rotational evolution of young VLM stars and BDs, a key ingredient of our work is the disk lifetime of individual stars. This parameter is related to the stellar mass accretion rate value, a quantity that changes with time. As in Paper I, we set a mass accretion rate threshold Ṁ_acc,th given by (2)where t_th is the threshold time scale. The time dependency of the mass accretion rate comes from the self-similar accretion theory by Hartmann et al. (1998). If a star has Ṁ_acc ≥ Ṁ_acc,th it will have a disk. Otherwise, it will be diskless. In Table 1 we show the threshold values obtained assuming t_th = 2.2 Myr.

The evolution equations for the stars assumed to be disk locked are as follows, for t<t_disk, where t_disk is the disk lifetime, I(t) = M_∗k(t)R(t)² is the moment of inertia, k is the gyration radius, R is the radius, and ω = 2π/P is the angular velocity of the star. As in Paper I, the stellar parameters are obtained from the stellar evolutionary models of Baraffe et al. (1998). Equation (3) expresses the constancy of the rotational period for a star with a disk, a behaviour known as the “disk locking hypothesis” proposed initially in the context of T Tauri stars by Koenigl (1991) and required by stellar angular momentum evolution models (e.g. Gallet & Bouvier 2013, 2015). According to this hypothesis, even though the star is accreting it does not spin up owing to its magnetic interaction with its surrounding accretion disk; however, the exact mechanism behind this interaction is not well known (see e.g. Bouvier et al. 2014). Equation (4) gives the temporal evolution of the mass accretion rate derived from fits of different young clusters (see e.g. Hartmann et al. 1998; Venuti et al. 2014; and Antoniucci et al. 2014). With this time-dependency, we were able to reproduce well the observed disk fraction as a function of time (see Fig. 4 of Paper I). However, there are controversies about the right value of the time exponent (e.g. Manara et al. 2012; and Rigliaco et al. 2011). Equation (5) gives the specific angular momentum for a star with a disk. Since the stellar angular velocity is constant, j decreases with time as a result of the contraction of the stellar radius. We are also assuming solid body rotation up to 12.1 Myr.

For diskless stars (with Ṁ_acc<Ṁ_acc,th), the equations are valid for t ≥ t_disk. We impose conservation of angular momentum (Eq. (8)) since there is no mechanism able to remove it until the star reaches the zero age main sequence (Scholz et al. 2009) at a time beyond the maximum time considered in our simulations (12.1 Myr). The star is no longer locked to its disk and it is free to spin up (Eq. (6)). The stellar rotational period will decrease because of the radius contraction. The mass accretion rate is set equal to the mass accretion rate threshold for the corresponding mass bin (Eq. (7)). This is an arbitrary choice, but it mimics the lower mass accretion rate limit of the observations. In the initial mass accretion rate distributions there are stars that have Ṁ_acc<Ṁ_acc,th, and for these stars we also impose Eq. (7). This causes a bump in the initial mass accretion rate distributions (see Fig. 3 in Paper I).

In the next section, we analyse the disk fraction as a function of time obtained with the set-up presented above.

3. Results

3.1. Disk fraction

The disk fraction as a function of time is an important diagnostic for disk dispersal mechanisms, disk lifetimes, and hence planetary formation. Disk fractions have been measured in different young clusters. The disk detection is based on accretion indicators, such as Hα emission, and/or dust indicators, such as near- and mid- IR excess, with a recent predominance of mid-IR IRAC colours. According to Mamajek (2009), it can be described by an exponential law with a characteristic time of 2.5 Myr. Pfalzner et al. (2014) raise several concerns¹ about the relevance and the reliability of the data used to build the disk fraction and consider that disk lifetimes can be much longer than estimated from the observational samples. Also, the disk fraction as a function of time depends on the correct determination of the age of the star or of the cluster to which it belongs, which can be quite uncertain (Bell et al. 2013; Kraus et al. 2015; Herczeg & Hillenbrand 2015).

With these uncertainties in mind, we compare disk fractions obtained from our simulations with available observational values shown in Table 2. We have assembled disk fractions for nine young stellar clusters, associations, and groups with ages in the range from 1 Myr to approximately 10 Myr. We focus on M stars mostly later than M3 and we show disk fractions for BDs, VLM stars, and all the stars (BD + VLM) together. We note that the total disk fraction for the ONC obtained by Rodríguez-Ledesma et al. (2010) is much smaller than the values for the ONC itself (found by Davies et al. 2014), Taurus, and Cha I, all regions with ages around 2 Myr. This is probably due to the use of near-IR excesses as a disk indicator which could lead to a low level of disk detections for VLM stars and BDs (see e.g. the discussion in Cody & Hillenbrand 2010). The values for 25 Ori and TW Hya, on the other hand, can suffer from low number statistics, and for Upper Scorpius the problem is the great uncertainty in its age.

The relevant parameters to constrain the disk fraction at a given age in our simulations are the time exponent of the mass accretion rate (Eq. (4)), the value of the mass accretion rate threshold (Eq. (2) and Table 1), and the dispersion of the mass accretion rate distribution (0.8 dex in our simulations). In Fig. 1 we plot two numerical disk fractions as a function of time superimposed on the total disk fractions from the regions listed in Table 2 and on two exponential laws with short (2.5 Myr) and long (5.5 Myr) characteristic times. We note that they lie between the two theoretical curves although they cannot be fitted by a simple exponential law. One of them was obtained taking t_th = 2.2 Myr. Its initial disk fraction is 74% and a disk fraction of 50% is reached at 2.2 Myr. The curve is very similar to that obtained in Paper I since the relevant parameters are the same (the time exponent of the mass accretion rate and its dispersion) or are scaled in the same manner (mass accretion rate values and mass accretion rate threshold).

Fig. 1 Disk fraction as a function of time obtained from our simulations taking tth = 2.2 Myr (solid black line) for VLM stars and BDs, and assuming different mass accretion rate threshold time scales for VLM stars tth,VLM = 2.2 Myr and for BDs tth,BD = 3.2 Myr (black dash-dotted line). The coloured dashed lines are exponential decay laws expected from disk e-folding times 2.5 Myr (red line) and 5.5 Myr (blue line). The coloured dots are data for nine young nearby clusters, associations, and groups taken from Table 2.

Fig. 2 Disk fraction as a function of time for VLM stars (top panel) and BDs (bottom panel). Symbols are for disk fraction of BDs (circles) and VLM stars (triangles) in nine young clusters.

When we observe the disk fractions for BDs and VLM stars shown in Table 2 we note that for the ONC, σ Ori, 25 Ori, and TW Hya the values differ significantly. In the literature there are some claims that low-mass stars actually have lower disk lifetimes when compared to BDs (e.g. Luhman et al. 2008a; Luhman & Mamajek 2012; and Downes et al. 2015) and some claims against this idea (e.g. Luhman et al. 2005; Damjanov et al. 2007; and Scholz & Jayawardhana 2008).

We take this into account in one of our models (model M3, see next sections), setting t_th,VLM = 2.2 Myr for VLM stars (the same value considered above) and t_th,BD = 3.2 Myr for BDs. The resultant disk fraction is the other black curve shown in Fig. 1. The initial disk fraction in this case is slightly higher, 77%, and a disk fraction of 50% is reached at ~2.4 Myr. In Fig. 2 we plot separately the disk fractions of VLM stars and BDs. As expected, the BD disk fraction curve falls off more slowly than the VLM disk fraction curve. In the analysis of model M3, we verify whether this difference in disk lifetimes is able to better describe the rotational periods observed for the lowest mass population of young clusters.

3.2. Models

We analyse four models (Table 3). These models differ by their period distributions, the disk lifetimes of VLM stars and BDs, and the disk locking hypothesis. Our final goal is to reproduce the rotational properties of the least massive population of young clusters which rotate faster in general than their solar-type counterparts. We want to investigate how the initial period conditions and the presence of a locking disk influence these properties. We compare our results with the period distributions and median period values of the ONC (Rodríguez-Ledesma et al. 2010; Davies et al. 2014), the σ Ori cluster (Scholz & Eislöffel 2004; Cody & Hillenbrand 2010), NGC 2362 (Irwin et al. 2008), and the Upper Sco association (Scholz et al. 2015). We also compare the disk fraction as a function of the period to verify whether VLM stars and BDs with disks rotate more slowly than diskless objects, as is the case for solar-type stars. We also investigate whether the period-mass relation holds and how the angular momentum evolves with time at this mass range. In general, the rotational properties of very low-mass objects are different from those of their more massive counterparts. They rotate faster than solar-type stars and their period distributions are not bimodal (e.g. Scholz et al. 2015 and Downes et al. 2014), although Rodríguez-Ledesma et al. (2010) find two distinct period peaks for the VLM stars of the ONC. Therefore, it is not clear whether the disk locking is as efficient at these mass ranges as it is for more massive stars (Lamm et al. 2005).

3.2.1. Model M1: reference model

We start with the simplest possible model, M1, which has analogous properties to those of the M1 model from Paper I (cf. Table 3). Both disk and diskless stars share the same initial period distribution, the disk lifetime is the same for VLM stars and BDs, and all the disk stars are locked, meaning that their periods are constant. However, we took different period distributions for VLM stars and BDs since it seems that they rotate faster (Bouvier et al. 2014).

In Fig. 3 we show the period distributions for BDs and VLM stars at 1.0 Myr, 2.1 Myr, 3.1 Myr and at 5.1 Myr and 10.1 Myr. We superimpose observational period distributions obtained from several works in the literature. We use data from Rodríguez-Ledesma et al. (2010 – RL10) and Davies et al. (2014 – D14) for the ONC; for the σ Ori cluster, we take data from the combined sample of Scholz & Eislöffel (2004 – S04) and Cody & Hillenbrand (2010 – C10); for NGC 2362 we use the sample from Irwin et al. (2008 – I08); and for the Upper Sco association we use data from Scholz et al. (2015 – S15). The data are summarized in Table 4.

Fig. 3 Period distributions obtained from model M1 for disk (light grey), diskless (grey), and disk and disless stars (bluish grey) for BDs (left panels) and VLM stars (right panels). The distributions are shown at 1.0 Myr (1st row), 2.1 Myr (2nd row), 3.1 Myr (3rd row), and at 5.1 Myr and 10.1 Myr (4th row). The binsize is 1.0 day. Superimposed on the histograms are observational data – squares for disk objects, crosses for diskless objects, and filled circles for objects not classified as disk or diskless – at 2.1 Myr from RL10 (red symbols) and D14 (blue symbols) for the ONC; at 3.1 Myr using the combined sample of S04 and C10 (blue symbols) for σ Ori; at 5.1 Myr from I08 (blue circles) for NGC 2362; and at 10.1 Myr from S15 (blue symbols) for Upper Sco. Error bars for all observational data were calculated following Burgasser et al. (2003).

At 1.0 Myr, the fraction of disk objects relative to diskless objects is the same for both mass regimes, ~ 74%, as established in the model (see also Fig. 1). The overall properties of the rotational period distributions (see Table 3) are the same for the two classes of objects (disk and diskless), but are different for BDs and VLM stars; the bulk of BDs rotate faster than VLM stars. The period medians for disk and diskless objects have the same value, 2.0 days for BDs and 3.7 days for VLM stars. As the system evolves stars with disks maintain their rotational rate, but diskless stars are expected to spin up. Then the diskless distribution should move to the left toward low periods, but the median of the period for disk objects should not change. The number of diskless stars is supposed to increase and the rate of disk to diskless stars should decrease.

Fig. 4 Disk fraction as a function of period (black circles) for BDs (left panels) and VLM stars (right panels) at 2.1 Myr (top panels) and at 3.1 Myr (bottom panels). The period bins are P ≤ 2.5 d, 2.5 d < P ≤ 5.0 d, and P > 5.0 d for BDs and P ≤ 2.5 d, 2.5 d < P ≤ 5.0 d, 5.0 d < P ≤ 7.5 d, and P > 7.5 d for VLM stars. Our simulations at 2.1 Myr are compared to observations of the ONC by RL10 (red squares) and D14 (blue triangles) and, at 3.1 Myr, they are compared to the combined σ Ori sample of S04 and C10 (red squares) The model’s error bars are equal to 1σ above and below the fraction value and they were calculated through resampling. The RL10 error bars were obtained from Rodríguez-Ledesma et al. (2010). For the observational sample at 3.1 Myr, the error bars were calculated following the Bayesian approach of Burgasser et al. (2003).

Our results (Fig. 3) confirm all these features for both mass regimes. The period distribution of BDs at 2.1 Myr is compared to ONC data from RL10. Their median values are 1.6 days for disk objects and 1.7 days for diskless objects against 2.0 days and 1.8 days in our simulations (cf. Table 5). We ran a χ² test between our data and the RL10 data, and the results show that the histogram values for stars with disks match, but they do not agree for diskless stars. At the VLM mass regime, they obtain very different values for the period medians of disk (4.2 days) and diskless stars (2.4 days), which disagrees with our results and with the data from D14. Our simulations give period medians of 3.7 days for disk stars and 2.7 days for diskless stars, while D14 obtain 3.7 days and 2.6 days for disk and diskless stars, respectively. The χ² test between the model and the RL10 VLM data show that the probability of the points matching is smaller than a 0.01 confidence level. We ran K-S tests that give probabilities of 59% and 18% that our disk and diskless samples and those from D14 are derived from the same population. It should be noted that although the period distributions are not bimodal and the peak values for disk and diskless objects are the same, the medians present different values but not as separated as the results from RL10 indicate.

At 3.1 Myr, the medians for diskless stars have moved to even smaller values (1.6 days for BDs and 2.2 days for VLM stars). When we compare our results with data for the combined catalogue of the σ Ori cluster we see a poorer agreement. The K-S tests give probabilities of 58% for disk and 10% for diskless BD objects and 2% for disk and 8% for diskless VLM stars that the samples share the same population. They also obtain lower medians than we do (Table 5).

At 5.1 Myr, data from I08 for VLM stars are not classified in terms of disk and diskless stars. The median value of their sample is 1.9 days, the same value we obtain for the median of the combined population of disk and diskless VLM stars. However, the probability that our data and I08 for NGC 2362 share the same population is only 16%.

Finally, at 10.1 Myr, the median values obtained from our results, 2.0 days for disk and 0.8 day for diskless BDs, are respectively around 20% smaller and 13% greater than the median values obtained by S15 for the Upper Sco association of 1.8 days for stars with disks and 0.96 day for diskless. The K-S tests indicate probabilities of 64% (for stars with disks) and 91% (for diskless stars) that the two samples have their origin in the same population.

In Fig. 4 we analyse the dependency of disk fractions on the rotational period. We show our results at 2.1 Myr and at 3.1 Myr and compare them to the observational data of RL10 and D14 for the ONC at 2.1 Myr and to the combined sample of S04 and C10 for σ Ori at 3.1 Myr. The period bins are the same as those proposed by RL10. The disk fractions of model M1 increases smoothly towards longer periods implying that stars with disks indeed rotate more slowly than diskless ones on average.

From χ² tests, we find that at 2.1 Myr our results agree with the ONC sample from RL10 for BDs and from D14 for VLM stars, but not with the results of RL10 at this mass range, within a 0.01 confidence level most probably owing to their much lower disk fraction values at shorter periods. This can be due to their use of near-IR excess as an indication of the presence of disks which, as we have mentioned before, misses many of the detections.

When we compare our results at 3.1 Myr with the σ Ori data the agreement is not as good, although χ² tests do not reject the hypothesis that the two data sets match within a 0.01 confidence level. In spite of this, the maximum value of the period observed in σ Ori is 3.1 d, while in our simulations there are stars with longer periods. For VLM stars, we also observe the same trend of increasing disk fractions with period, again except for the 3rd data point which, however, exhibits a very large error bar².

Furthermore, no clear bimodality appears when we analyse the rotational period as a function of the mass accretion rate. In Fig. 5 we plot these quantities at an age of 3.1 Myr to make the comparison with the data of the σ Ori cluster from C10 easier. The mass accretion rate values are normalized to the mass accretion rate threshold, which has a different value depending on the mass of the object (see Eq. (2)). Because of this and because Ṁ_acc(M_∗,t) = Ṁ_acc,th(M_∗) for all diskless stars, these will be clustered at Ṁ_acc/Ṁ_acc,th = 1; however, for stars with disks this ratio will be higher than 1.0. We see no correlation of Ṁ_acc with the period even for stars with disks, the same conclusion reached by C10 based on their observations. We observe no paucity of stars with disks at short periods, as can be seen in Fig. 9 of Paper I where there are no stars with disks below P ~ 1.0 d at the same age.

Fig. 5 Rotational periods versus mass accretion rate normalized to the mass accretion threshold at 3.1 Myr. Circles are for disk objects, while crosses are for diskless objects. Each colour corresponds to a different object mass, shown on the colourbar at the right of the plot.

We also investigate the period-mass relation for model M1. In Fig. 6, we plot mass xlog P for 200 objects randomly chosen from model M1. We also show the 75th percentiles of log P, the best linear fits for these percentiles, and the slopes of the fits from 1 Myr to 100 Myr. These percentiles and the fits were calculated for all objects with 0.1 M_⊙ ≤ M ≤ 0.4 M_⊙. The error bars of the slopes are equal to the standard errors of the estimate of the least-square fits. We see that the slopes are approximately equal to 0 at all ages and we conclude that there is no dependency of the period with mass in model M1. For comparison, we plot the slopes obtained by Henderson & Stassun (2012) for NGC 6530 (1.65 Myr), ONC (2.0 Myr), NGC 2264 (3.0 Myr), NGC 2362 (3.5 Myr), IC 348 (4.5 Myr), NGC 2547 (40 Myr), and NGC 2516 (150 Myr) for the low-mass stars from 0.1 M_⊙ to 0.5 M_⊙, and we see a clear positive correlation that increases with age. Cody & Hillenbrand (2010) also found a positive correlation between period and mass in their sample of σ Ori stars (see also Bouvier et al. 2014). In Paper I we analysed the period-mass relation for stars from 0.3 M_⊙ to 0.5 M_⊙ and we concluded that a correlation only exists if we relax the disk locking hypothesis for the lowest mass stars in those simulations, i.e. the 0.3 M_⊙ stars. In model M4 (Sect. 3.2.4) we relax the disk locking hypothesis for all stars except those of 0.4 M_⊙ and we verify whether we can obtain a correlation similar to that obtained in Paper I.

Interestingly, when we construct a period-mass relation for BDs, the slope we obtain is the opposite observed for stars, as can be seen in Fig. 7. The slope in this case is the fit of the 75th percentiles of the log (P) for M_∗ ≤ 0.07 M_⊙. This is caused by the more rapid decrease in the moment of inertia I of higher mass BDs than of lower mass ones, as illustrated in Fig. 8. On the other hand, this does not happen at the 0.1–0.4 M_⊙ mass interval where I decreases at the same rate for all masses.

Fig. 6 Left panels: mass versus the logarithm of the rotational period of a sample of 200 stars extracted from model M1 (black circles) at four different ages superimposed by the 75th percentiles of the log (P) (cyan diamonds). The best linear fits of the percentiles for 0.1 M⊙ ≤ M∗ ≤ 0.4 M⊙ are shown as red lines. Right panel: slope of the best linear fits plotted as a function of age (red circles) superimposed with the values obtained by Henderson & Stassun (2012) for seven clusters of different ages (black circles).

In Fig. 9 we show the specific angular momentum distributions for BD and VLM stars at 2.1 Myr. We also show data from D14 for stars in the ONC with spectral type later than M2. These were calculated using period values from D14, but gyration and stellar radii were obtained from Baraffe et al. (1998) models interpolated for the ages and masses of the D14 ONC sample.

The observational distributions show peaks at higher log j values when compared to the simulations. The medians of log j are also displaced by ~ 0.1. Since the agreement between the model’s and D14’s period medians is very good (as seen in Fig. 3) the discrepancies stem from the broader range of masses and ages found in the sample compared to those present in the model. In the numerical distributions, which for BDs is single-peaked with a maximum at log j ≃ 15.9–16.1. For VLM stars, the peaks of the distributions of disk and diskless objects are different: log j = 16.1–16.3 for stars with disks and log j = 16.3−16.5 for diskless stars. This is related to the different internal structure of BD and VLM objects since the period distributions for both are single-peaked (cf. Fig. 3). Indeed, Fig. 8 shows that low-mass BDs do not contract much over the first few Myr so that their moment of inertia remains constant, while this is not the case for VLM stars which steadily contract over time and therefore start to spin up earlier.

From K-S tests, we obtain that the specific angular momentum distributions for the stars with disks from the ONC and from this work have a 6% probability of sharing the same origin, while for the diskless distributions the probability is 20%.

According to the conditions imposed in our simulations (Eqs. (3)–(6)), the specific angular momentum is constant for diskless objects, but it decreases for disk ones owing to the stellar radius contraction at the pre-MS phase. This can be seen in Fig. 10 where we plot j as a function of age for 100 objects extracted from model M1. We overplot the specific angular momenta of D14 ONC sample taking into account the different ages of the stars. Although initially the specific angular momentum has the same value for disk and diskless stars, which is expected since they share the same initial rotational conditions, in the course of the system’s evolution the median values become lower for disk objects, which decrease with time as ⟨ j ⟩ ∝ t^γ with γ = −0.71 ± 0.02, while for diskless ones we obtain γ = −0.20 ± 0.01, a much weaker time dependency. This is not the same result we obtained in Paper I for solar-type stars where γ_d = −0.65 and γ_dl = −0.53 for disk and diskless stars, respectively. Also, the time exponent for stars with disks is smaller than is expected from polytropic models (j ∝ t^-0.66). For diskless stars it is the opposite; the value shows a much weaker time dependency, closer to the constancy expected individually. As can be seen in Fig. 8, for M_∗ ≥ 0.2 M_⊙ the moment of inertia falls off steadily, while for less massive stars and BDs it remains constant for a time scale that increases with decreasing stellar mass. Then, a VLM star with a disk will experience a stronger decrease in its specific angular momentum than a BD disk object. As said before, this is due to the different internal structure evolution of BDs and VLM stars. We also examine separately the slopes for BDs and stars with M_∗ ≥ 0.2 M_⊙. We obtain a slope for stars with disks equal to − 0.690 ± 0.005, not very far from the value expected from the polytropic case, but this exponent is − 0.56 ± 0.05 for BDs. We thus observe a flatter j decrease for accreting BDs since their momentum of inertia hardly changes up to 3 Myr.

Fig. 7 Slopes of the best linear fits to M∗ × log (P)75th plots as a function of age only for BDs from model M1.

Fig. 8 Moments of inertia obtained from Baraffe et al. (1998) normalized to the moment of inertia at 1.0 Myr for BDs (top panel) and stars from 0.1 to 0.4 M⊙ (bottom panel). Different colours show different mass values.

Fig. 9 Distribution of the logarithm of the specific angular momentum of BD (left panel) and VLM (right panel) disk (light grey bars) and diskless (grey bars) objects at 2.1 Myr obtained from model M1. We superimpose data from the D14 ONC sample for disk (squares) and diskless (crosses) stars of spectral type later than M2.

Fig. 10 Left: specific angular momentum evolution for 100 objects randomly chosen from model M1. Simulated data are shown as light grey symbols, squares for disk objects and crosses for diskless ones. Two stars are highlighted as magenta filled circles (when in the disk phase) or crosses (when in the diskless phase) to illustrate the evolution of a single star. The black symbols are for the D14 ONC sample of stars of spectral type later than M2 whose j values were recalculated using gyration and stellar radii from the stellar evolutionary models of Baraffe et al. (1998). Squares represent Class II stars, crosses Class III stars. Right: logarithm of age versus the logarithm of median j for disk objects (black squares) and diskless objects (black crosses). Error bars (black for disk objects, red for diskless ones) show the rms of log ⟨ j ⟩ values around the median. The γ values express the time dependency ⟨ j ⟩ ∝ tγ.

In summary, model M1 is able to reproduce (i) the period distributions and medians of the ONC provided by D14, of NGC 2362, and of the Upper Sco association; (ii) the period × disk fraction for the ONC obtained by RL10 for BDs and by D14 for VLM stars; and (iii) the absence of correlation between the mass accretion rate and period for stars with disks in agreement with C10. However, it fails to reproduce (i) the period distributions and medians of the ONC provided by RL10 and of the combined sample of σ Ori; (ii) the period × disk fraction for the same sample; and (iii) the period-mass relation. In order to try to improve the agreement with these three points, we will investigate other models that address these discrepancies.

3.2.2. Model M2: initial bimodal period distributions

In order to try to reproduce the bimodal period distribution of VLM stars in the ONC seen by Rodríguez-Ledesma et al. (2010) we ran model M2, which differs from model M1 by the initial period distributions (see Table 3). We set a bimodal period distribution at 1.0 Myr, but only for VLM stars. For BDs the initial period distribution is the same used in model M1. For VLM stars the mean period of the diskless distribution is 2.0 days and the minimum period is 0.2 day. For stars with disks, the mean period is 3.0 days and the minimum period is 0.5 day. In Fig. 11 we show the distributions of rotational periods for VLM stars at 1.0 Myr, 2.1 Myr, 3.1 Myr, and 5.1 Myr. The bimodality imposed at the beginning of the simulations is clearly visible at later times. In terms of medians (see Table 5), there is an improvement between the model and the observational data from RL10, but this is not the case at 3.1 Myr and for the σ Ori data. At 5.1 Myr there is no change in what was seen with model M1.

The χ² test between model M2 and the RL10 data shows that the fit between the diskless VLM histogram data is now acceptable, but this is not true for the disk VLM histogram data which again do not match the RL10 observations. There was no change for the BD distribution. The K-S test for D14 data now show probabilities of 40% that the disk distribution and 0.4% that the diskless distribution and the numerical distributions are drawn from the same population against 59% and 18%, respectively, obtained from model M1 (see discussion at Sect. 3.2.1).

At 3.1 Myr, the K-S probability degraded for the VLM disk distribution which now is only 0.3%, while with model M1 it is 2%. For the diskless distribution, to the contrary, it is greatly improved, reaching 49% agreement with the observational σ Ori data (against 8% using model M1). At 5.1 Myr, the probability that the numerical VLM period distribution and the sample from I08 are the same is 6%, smaller than the 16% value obtained with model M1. Since we have not changed the distributions for the BD objects, the probabilities at 10.1 Myr have not changed.

Fig. 11 Period distributions for the VLM range for model M2 at 1.0 Myr, 2.1 Myr (two upper panels), 3.1 Myr, and 5.1 Myr (two lower panels). Bars, symbols, and error bars have the same colour and form coding of Fig. 3.

When we analyse disk fractions as a function of the period for VLM stars we note that – compared to model M1 – the agreement is better with the RL10 sample and worse with the D14 one, although the RL10 data still do not fit the numerical data. The match with the D14 data is still better than it is with the RL10 data that justified model M2.

At 3.1 Myr, the χ² value slightly decreases and the hypothesis that the two values match cannot be discarded.

Thus, model M2 reproduces the bimodality at the VLM range at 2.1 Myr, but the agreement with the RL10 data is not as good as is desirable. The bimodal period distribution of VLM stars is observationally controversial since neither D14 nor C10 found it for their samples. Therefore, we consider that model M2 with its initial bimodal distribution is not adequate. Instead, in models M3 and M4 we keep the same initial period distributions for disk and diskless stars. In these models we investigate different disk lifetimes for BD and VLM stars and address the disk locking hypothesis.

3.2.3. Model M3: different disk lifetimes for BD and VLM stars

Several observations point to a longer lifetime for disks around BDs (e.g. Luhman et al. 2008a; Luhman & Mamajek 2012 and Downes et al. 2015; see Sect. 3.1 and Table 2). To examine its impact on our simulations, we ran model M3 (see Table 3). The disk fraction of all objects as a function of time is shown in Fig. 1 and the disk lifetimes for BD and VLM stars are depicted in Fig. 2. Except for this property, the set-up of model M3 is the same as model M1.

In comparison to model M1 in Fig. 3, the period distributions are practically unaltered except that there are more disk than diskless BDs. When we analyse the disk fractions as a function of the period, there is no change for VLM stars compared to model M1 (cf. Fig. 4), but the disk fractions change at the BD regime (Fig. 12). Now they are higher than in model M1. At 2.1 Myr, they reach 70% which is well above the RL10 fractions, which are around 40%. At 3.1 Myr, the numerical fractions show a smooth rising towards longer periods and present values between 45% and 65%, far from the σ Ori fractions. This is a direct consequence of the longer disk lifetime at the BD mass regime. Again, we do not obtain any correlation between the mass accretion rate and the period.

Fig. 12 Disk fractions as a function of the period for BDs from model M3.

We conclude that the significant increase in the number of BDs with disks with short periods predicted by model M3 contradicts the observations. Model M3 does not adequately represent the rotational evolution of BDs.

3.2.4. Model M4: no disk locking for stars less massive than 0.3 M_⊙

Model M4 imposes no disk locking for all objects except the 0.4 M_⊙ stars. As shown in Table 3, this model has the same initial period distributions as model M1 and also the same mass accretion rates and disk lifetimes.

Owing to the results of models M1 to M3, for model M4 we should expect the total disk fraction to follow the curve shown in Fig. 1, which is indeed the case. The period distributions, however, are different from those of previous models (Fig. 13). The free spin-up of all BDs erases the differences between the medians of disk and diskless BD objects and the two curves show the same values at all ages (see Table 5). As can be seen, in comparison with the observations, the agreement is improved at 2.1 Myr with the RL10 BD medians, but at 3.1 Myr and at 10.1 Myr it becomes worse. At 3.1 Myr, the diskless BD median is greater than the observational value, while at 10.1 Myr the medians have moved to shorter periods and thus farther away from the observational ones.

Concerning the period distributions, the agreement with RL10 data seen through the χ² tests is below the confidence level. The K-S tests show a probability of only 3% of compatibility with the D14 period distribution for stars with disks and of 11% with the same distribution but for diskless stars.

At 3.1 Myr the probabilities of sharing the same population are 90% and 16% among observational and numerical BD disk and diskless period distributions, and 47% and 30% for disk and diskless VLM stars, respectively.

At 5.1 Myr the disagreement between the distributions is great: the K-S test provide a probability of 0.02% that the observational and theoretical samples are the same.

Finally, at 10.1 Myr we obtain a probability of 0.2% that the disk BD period distributions agree and of 46% that the diskless BD distributions are the same.

So, except at 3.1 Myr, model M4 is not as good as model M1 in reproducing the observational period distributions.

Fig. 13 Same as Fig. 3 but for model M4.

In Fig. 14 we show plots of disk fraction versus period. The results are different from what was seen in model M1 and in model M3. The disk fractions for BDs remains approximately constant around 50% at 2.1 Myr in agreement with the observations and around 40% at 3.1 Myr. For VLM stars, we do not observe the high disk fractions at longer periods that are seen in model M1, but this is not very far from what is observed. The χ² tests show that model M4 fits the BD and VLM disk fractions from the σ Ori sample at 3.1 Myr. At 2.1 Myr, the agreement with RL10 BD data is inside the confidence level and the same occurs with the fit between our data and the D14 sample. However, this is not the case with VLM data from RL10.

The specific angular momentum evolution on the other hand is drastically different. The time dependency exponents are now very low, and for the BD sample they are practically zero in accordance with the simulation’s assumption of constant specific angular momentum for diskless stars. In contrast with model M1 for which the sequential release of a group of stars from the disk causes a decrease in the mean specific angular momentum, in model M4 – since all stars except those of 0.4 M_⊙ are already free to spin up – the group as a whole evolves at almost constant angular momentum. In Table 6 we summarize the γ exponents obtained so far with our models. The probability that the numerical and D14 specific angular momentum distributions are similar rises to 31% for disk VLM stars and 41% for diskless VLM stars, against 6% and 20%, respectively, obtained with model M1.

Fig. 14 Same as Fig. 4 but for model M4.

We analysed the period-mass relation in model M4, with no disk locking for stars less massive than 0.4 M_⊙. In Fig. 15 we show the slopes of plots M × log P₇₅ (cf. Fig. 6). In comparison to what was obtained with model M1, we observe a correlation between the rotational period and the mass for ages older than 3.0 Myr in agreement with the Henderson & Stassun (2012) trend with, however, a shallower slope. We thus confirm the results of Paper I that suggest that VLM stars and BDs are not as efficiently disk locked as solar-type stars when accreting from their disks.

In summary, model M4 is worse in reproducing the period distributions, but is equivalent in accounting for the disk fractions as a function of the period in comparison to model M1. Model M4 predicts practically no angular moment evolution, which is in disagreement with the observations of D14. However, the stronger correlation between the mass and the period favours this model.

4. Conclusions

We have analysed four Monte Carlo models in order to investigate the rotational behaviour of VLM stars and BDs with masses below 0.4 M_⊙ during the early pre-main sequence.

Stars in this mass range rotate differently from their solar mass counterparts. In order to fit the period distributions for BDs, we had to consider that they rotate faster than VLM stars, whose period distribution is also centred at a shorter period compared to that of the solar mass stars considered in Paper I. There is no need to consider bimodal distributions for BDs and VLM stars at the beginning of the simulations, contrary to what was obtained for solar mass stars. All this suggests a weaker influence of the disk on the angular momentum evolution of BD and VLM stars. It is possible to speculate that the different magnetic field topology found among the lowest mass stars is the cause of their faster rotation rate.

Fig. 15 Slopes of mass versus the logarithm of the 75th percentile of the rotational period plots at ages ranging from 1.0 Myr to 100 Myr obtained from model M4.

Of all four models outlined here, model M1 (with an initial single-peak period distribution for disk and diskless stars and BDs as faster rotators) best reproduces the observations. Considering an initial bimodal period distribution as was done in model M2 does not improve the agreement with the observational data. Model M1 nevertheless has shortcomings; most notably, it does not reproduce the period-mass relationship observed for VLM stars. Changing the average disk lifetime of BDs (model M3) does not improve the results. Relaxing the disk locking hypothesis for the lowest mass stars, as was required in Paper I, yields a period-mass trend in the same direction as observed, but is poorer in reproducing the evolving rotational distributions. Hence, none of the models developed here can reproduce all available observational constraints simultaneously, though they offer useful guidelines regarding initial period distributions, disk lifetimes, and the disk locking efficiency at very low masses. We also note that the observational results do not always seem consistent, and much remains to be done to properly characterize the rotational properties of young VLM stars and BDs.

Acknowledgments

M.J.V. would like to thank the financial support provided by CAPES (fellowship No. 2565-13-7) under the program “Science without borders” and by the project PROCAD – CNPq/CAPES number 552236/2011-0. J.B. acknowledges the support of ANR grant 2011 Blanc SIMI5-6 020 01 Toupies: Towards understanding the spin evolution of stars (http://ipag.osug.fr/Anr_Toupies/). The authors thank the referee for very useful suggestions that helped to improve the quality of the paper.

References

All Tables

All Figures

Fig. 1 Disk fraction as a function of time obtained from our simulations taking tth = 2.2 Myr (solid black line) for VLM stars and BDs, and assuming different mass accretion rate threshold time scales for VLM stars tth,VLM = 2.2 Myr and for BDs tth,BD = 3.2 Myr (black dash-dotted line). The coloured dashed lines are exponential decay laws expected from disk e-folding times 2.5 Myr (red line) and 5.5 Myr (blue line). The coloured dots are data for nine young nearby clusters, associations, and groups taken from Table 2.

In the text

	Fig. 2 Disk fraction as a function of time for VLM stars (top panel) and BDs (bottom panel). Symbols are for disk fraction of BDs (circles) and VLM stars (triangles) in nine young clusters.
In the text

Fig. 3 Period distributions obtained from model M1 for disk (light grey), diskless (grey), and disk and disless stars (bluish grey) for BDs (left panels) and VLM stars (right panels). The distributions are shown at 1.0 Myr (1st row), 2.1 Myr (2nd row), 3.1 Myr (3rd row), and at 5.1 Myr and 10.1 Myr (4th row). The binsize is 1.0 day. Superimposed on the histograms are observational data – squares for disk objects, crosses for diskless objects, and filled circles for objects not classified as disk or diskless – at 2.1 Myr from RL10 (red symbols) and D14 (blue symbols) for the ONC; at 3.1 Myr using the combined sample of S04 and C10 (blue symbols) for σ Ori; at 5.1 Myr from I08 (blue circles) for NGC 2362; and at 10.1 Myr from S15 (blue symbols) for Upper Sco. Error bars for all observational data were calculated following Burgasser et al. (2003).

In the text

Fig. 4 Disk fraction as a function of period (black circles) for BDs (left panels) and VLM stars (right panels) at 2.1 Myr (top panels) and at 3.1 Myr (bottom panels). The period bins are P ≤ 2.5 d, 2.5 d < P ≤ 5.0 d, and P > 5.0 d for BDs and P ≤ 2.5 d, 2.5 d < P ≤ 5.0 d, 5.0 d < P ≤ 7.5 d, and P > 7.5 d for VLM stars. Our simulations at 2.1 Myr are compared to observations of the ONC by RL10 (red squares) and D14 (blue triangles) and, at 3.1 Myr, they are compared to the combined σ Ori sample of S04 and C10 (red squares) The model’s error bars are equal to 1σ above and below the fraction value and they were calculated through resampling. The RL10 error bars were obtained from Rodríguez-Ledesma et al. (2010). For the observational sample at 3.1 Myr, the error bars were calculated following the Bayesian approach of Burgasser et al. (2003).

In the text

	Fig. 5 Rotational periods versus mass accretion rate normalized to the mass accretion threshold at 3.1 Myr. Circles are for disk objects, while crosses are for diskless objects. Each colour corresponds to a different object mass, shown on the colourbar at the right of the plot.
In the text

Fig. 6 Left panels: mass versus the logarithm of the rotational period of a sample of 200 stars extracted from model M1 (black circles) at four different ages superimposed by the 75th percentiles of the log (P) (cyan diamonds). The best linear fits of the percentiles for 0.1 M⊙ ≤ M∗ ≤ 0.4 M⊙ are shown as red lines. Right panel: slope of the best linear fits plotted as a function of age (red circles) superimposed with the values obtained by Henderson & Stassun (2012) for seven clusters of different ages (black circles).

In the text

	Fig. 7 Slopes of the best linear fits to M∗ × log (P)75th plots as a function of age only for BDs from model M1.
In the text

	Fig. 8 Moments of inertia obtained from Baraffe et al. (1998) normalized to the moment of inertia at 1.0 Myr for BDs (top panel) and stars from 0.1 to 0.4 M⊙ (bottom panel). Different colours show different mass values.
In the text

	Fig. 9 Distribution of the logarithm of the specific angular momentum of BD (left panel) and VLM (right panel) disk (light grey bars) and diskless (grey bars) objects at 2.1 Myr obtained from model M1. We superimpose data from the D14 ONC sample for disk (squares) and diskless (crosses) stars of spectral type later than M2.
In the text

Fig. 10 Left: specific angular momentum evolution for 100 objects randomly chosen from model M1. Simulated data are shown as light grey symbols, squares for disk objects and crosses for diskless ones. Two stars are highlighted as magenta filled circles (when in the disk phase) or crosses (when in the diskless phase) to illustrate the evolution of a single star. The black symbols are for the D14 ONC sample of stars of spectral type later than M2 whose j values were recalculated using gyration and stellar radii from the stellar evolutionary models of Baraffe et al. (1998). Squares represent Class II stars, crosses Class III stars. Right: logarithm of age versus the logarithm of median j for disk objects (black squares) and diskless objects (black crosses). Error bars (black for disk objects, red for diskless ones) show the rms of log ⟨ j ⟩ values around the median. The γ values express the time dependency ⟨ j ⟩ ∝ tγ.

In the text

	Fig. 11 Period distributions for the VLM range for model M2 at 1.0 Myr, 2.1 Myr (two upper panels), 3.1 Myr, and 5.1 Myr (two lower panels). Bars, symbols, and error bars have the same colour and form coding of Fig. 3.
In the text

	Fig. 12 Disk fractions as a function of the period for BDs from model M3.
In the text

	Fig. 13 Same as Fig. 3 but for model M4.
In the text

	Fig. 14 Same as Fig. 4 but for model M4.
In the text

	Fig. 15 Slopes of mass versus the logarithm of the 75th percentile of the rotational period plots at ages ranging from 1.0 Myr to 100 Myr obtained from model M4.
In the text