© ESO, 2016

1. Introduction

The disagreement between theoretical and observed parameters of young magnetically active and of fully convective or almost fully convective low-mass stars remains one of the main long standing problems in stellar physics. Current investigations focus on the inhibition of the convective transport by interior dynamo-generated magnetic fields and/or by the blocking of flux at the surface by cool magnetic starspots (e.g. Mullan & MacDonald 2001; Chabrier et al. 2007; Feiden & Chaboyer 2013, 2014; Jackson & Jeffries 2014), which produce an increase in stellar radius and a decrease in T_eff. The same effect is also thought to be linked to the observed correlation between Li abundance and rotation (e.g. Somers & Pinsonneault 2014, 2015b,a; Jackson & Jeffries 2014). The consequences of these discrepancies are manifold. These include, for example, the determination of the mass and the radius of exoplanets, whose accuracy depends on that of the hosting star (e.g. Henry 2004; Mann et al. 2015), the age estimate of young open clusters (e.g. Soderblom et al. 2014; Somers & Pinsonneault 2015a), and the mass-luminosity relationship for magnetically active low-mass stars.

Fundamental determinations of stellar masses and radii with a 3 percent accuracy or better are provided by the light-curve analysis of detached eclipsing binaries (e.g. Torres et al. 2010; Feiden & Chaboyer 2012). Interferometric angular diameter measurements of single stars are available today for tens of stars (e.g. Boyajian et al. 2012) with diameters measured to better than 5 percent.

Statistical methods based on the product of P and vsini, which produces the projected radius Rsini (Sect. 3.1), and the assumption of random orientation of the spin axis (e.g. Jackson et al. 2009) have the advantage of providing mean radii estimates for a large number of (coeval) single stars independently of distance, interstellar absorption, and models. No evidence of preferred orientation of the spin axis in open clusters has been found so far (e.g. Jackson & Jeffries 2010, and references therein), and therefore the method seems to be sound in this respect. The main difficulty is that the data sample is always truncated at a combination of sufficiently low inclination angle i and low equatorial velocities v_eq. In these cases, depending also on the spectral resolution, vsini cannot be derived and only an upper limit can be given. A low i may also cause difficulties in measuring P and therefore there may be cases in which either one or both vsini and P cannot be measured. At the other extreme, ultra-fast rotator spectra can be so smeared by the rotational broadening that in some cases only a lower vsini limit can be given.

To take a low Rsini truncation into account, Jackson et al. (2009) considered a cut-off inclination such that stars with lower inclination yield no Rsini, and they corrected the average sini accordingly. Mean radii are then derived by taking the average of the ratio Rsini/⟨sini⟩ in suitable magnitude bins.

Here we present a new method, based on the survival analysis concept (Klein & Moeschberger 2003), that makes use of the whole information content of the dataset by also considering upper and lower limits and data truncation. Data may also come from inhomogeneous estimates, like those in which vsini upper and lower limits are obtained from different analyses and instrumentation, as long as they are not affected by significant biases. Uncertainties due to surface differential rotation (SDR) are also estimated, with the most likely values derived from the recent work of Distefano et al. (2016). The method is applied for the first time to the rich dataset available for the Pleiades.

In Sect. 2 we present the data used in this work. The method is described in Sect. 3. The results obtained for the Pleiades dataset are discussed in Sect. 4. We draw our conclusions in Sect. 5.

2. Data

For this work rotational periods and memberships from Hartman et al. (2010) are used. Measurements of vsini are taken from Stauffer & Hartmann (1987), Soderblom et al. (1993), Queloz et al. (1998) and Terndrup et al. (2000). Magnitudes are adopted from Stauffer et al. (2007).

Fig. 1 Hartman et al. (2010)P dataset of the Pleiades (blue filled circles). Upper panel: open red squares outline targets for which vsini is available and green stars those considered in this work. Lower panel: open red triangles indicate vsini from Queloz et al. (1998), green squares the values adopted from Terndrup et al. (2000), plum pentagons those from Stauffer & Hartmann (1987), and olive hexagons those from Soderblom et al. (1993). See text for details.

The P dataset of Hartman et al. (2010) comprises 383 stars and is 93 percent complete in the mass range 1.0 ≳ M/M_⊙ ≳ 0.7. Of these, 227 have measured vsini. Stars flagged as binaries in Hartman et al. (2010) were excluded from the sample. A total of 217 stars constitute the final working sample (Fig. 1).

Theoretical mass-radius-magnitude relationships were taken from Baraffe et al. (2015). Together with the adopted distance to the Pleiades of 136.2 pc (Melis et al. 2014), an age of 120 Myr (Stauffer et al. 1998), and an extinction of A_K = 0.01 mag (Stauffer et al. 2007), the models of Baraffe et al. (2015) were used to build bins in M_K magnitude corresponding to approximately regular intervals of mass, and then in the comparison of our results with the theoretical R-M_K relationship. The calculations reported in Sect. 4 were repeated, also assuming a distance to the Pleiades of 120.2 pc (van Leeuwen 2009). Our results, however, are more consistent with a distance of ≈136.2 pc, and therefore we report only the results obtained assuming this value (see also Soderblom et al. 2005).

We note that the fraction of stars with both P and vsini considered in the analysis with respect to the whole P dataset of Hartman et al. (2010) is very close to one down to M_K ≈ 4.3, corresponding to M ≈ 0.8 M_⊙. For fainter magnitudes this fraction decreases progressively to M_K ≈ 6 (M ≈ 0.5 M_⊙), below which there are no vsini measurements and very few and sparse P measurements (Fig. 1, upper panel). Furthermore, vsini values for fast rotators were mostly adopted from Soderblom et al. (1993), for slow-rotators with M_K ≲ 4.3 mostly from Queloz et al. (1998), while at fainter magnitudes the measurements are mostly those reported by Terndrup et al. (2000). Possible consequences of this inhomogeneity for our results are discussed in Sect. 4.

3. Method

3.1. Projected radius

When we assume spherical symmetry, the relationship between stellar radius R, stellar equatorial rotational period P_eq, and equatorial velocity v_eq is $\mathit{R}\mathrm{=}\mathrm{𝒞}{\mathit{P}}_{\mathrm{eq}}{\mathit{v}}_{\mathrm{eq}}\mathit{,}$(1)where when P_eq is in days, R in solar units, and v_eq in km s^-1. The analysis of photometric time-series of stars showing rotational modulation produces the measured period P (see Sect. 3.2), while the estimated rotational broadening from spectroscopic analysis provides projected rotational velocities vsini. For most of our sample deviations from spherical symmetry are expected to be negligible (e.g. Collins & Truax 1995). Furthermore, from our SDR estimate (Sect. 3.2) and the work of Reiners (2003), Reiners & Schmitt (2003), and von Eiff & Reiners (2010), who also took limb darkening into account, we estimate that the SDR effects on vsini are not larger than ≈1 percent and are therefore significantly smaller than our conservative estimate of the period uncertainties associated with SDR (see Sect. 3.2). For our purposes we therefore neglected the SDR effects on vsini and adopted the relationship $\mathrm{\left(}\mathit{v}\sin \mathit{i}\mathrm{\right)}\mathrm{=}{\mathit{v}}_{\mathrm{eq}}\sin \mathit{i,}$(2)where i is the inclination of the spin axis from the line of sight. Combining the measured vsini and P_eq, we obtain the projected stellar radius $\mathrm{\left(}\mathit{R}\sin \mathit{i}\mathrm{\right)}\mathrm{=}\mathrm{𝒞}{\mathit{P}}_{\mathrm{eq}}{\mathit{v}}_{\mathrm{eq}}\sin \mathit{i}\mathrm{=}\mathrm{𝒞}{\mathit{P}}_{\mathrm{eq}}\mathrm{\left(}\mathit{v}\sin \mathit{i}\mathrm{\right)}\mathit{.}$(3)The inclination angle i is unknown and therefore it is not possible to derive either v_eq from vsini or R from Rsini for each individual star. However, when the underlying probability density functions of R and i are known, it is possible to estimate the expected value of R, ℛ, from an ensemble of Rsini measurements (see Sect. 3.3).

3.2. Surface differential rotation

Surface differential rotation implies that magnetically active regions, associated with dark spots and bright faculae, rotate with different frequencies, depending on their latitudes. Multiple active regions at different latitudes broaden the periodogram peak and contribute to the uncertainty in the period, as discussed by Hartman et al. (2010), for instance. On the other hand, SDR also leads to a systematic error in determining the rotation period of the star since the degree of differential rotation and the latitude of the dominant active region are not known, which prevents us from relating the measured period to the equatorial period P_eq. For young rapidly rotating stars like the Pleiades, Hartman et al. (2010) assumed that the dominant active region groups can be assumed to be isotropically distributed, from which they estimated that the mean rotational period is ⟨P⟩ = 1.03P_eq. Following Kitchatinov (2005), for a solar-like SDR¹ in which spots are confined to latitudes |β| < 30°, Hartman et al. (2010) estimated ⟨P⟩ = 1.07P_eq.

Fig. 2 Fit to the mean period vs. equatorial period for the AB Dor young loose association.

Distefano et al. (2016) made use of long-term photometric monitoring to estimate SDR lower limits in members of young loose association from the P modulation itself. Considering the AB Doradus young loose association, whose age is very similar if not identical to the Pleiades, assuming a solar-like SDR and that the minimum period corresponds to P_eq, we fitted the observed ⟨P⟩ vs. P_eq with a linear relationship obtaining ⟨P⟩ = 1.04P_eq (Fig. 2), which is very similar to that estimated by Hartman et al. (2010).

We therefore assumed that the Hartman et al. (2010)P are representative of ⟨P⟩ and derived the equatorial period from the ⟨P⟩ = 1.04P_eq relationship. In this way, we obtained our best estimate of P_eq. Furthermore, to evaluate the uncertainties associated with the SDR, we considered at one extreme solid-body rotation, ⟨P⟩ = P_eq, and at the other extreme the ⟨P⟩ = 1.07P_eq relationship. Given the information derived from the AB Dor mean vs. equatorial period fit, such uncertainties likely overestimate the true uncertainties due to the SDR. The corresponding Rsini uncertainties obtained using Eq. (3)are then also overestimated and larger, overall, than expected from the propagation of P_eq and vsini uncertainties due to SDR (see Sect. 3.1). These can therefore be considered to include both P_eq and vsini uncertainties due to the SDR.

3.3. Distribution functions

The statistical analysis presented here is similar to that used extensively in the past to derive the expected distribution of v_eq from vsini measurements (e.g. Gaige 1993; Queloz et al. 1998; Spada et al. 2011). The application of these concepts to the analysis of the Rsini distribution in open clusters has some advantages with respect to the analysis of vsini. Stars of similar mass in a stellar cluster, which therefore have approximately the same age, are expected to also have a similar radius, and this can justify the assumption that radii are distributed normally around the expected value, with a standard deviation that includes the intrinsic spread and measurement uncertainties. In contrast, there are no such constraints on v_eq, which makes reconstructing its distribution conceptually more difficult. Following Chandrasekhar & Münch (1950), the y ≡ Rsini probability density function (p.d.f.) for a random orientation of spin axis can be therefore written as $\mathit{\phi }\mathrm{\left(}\mathit{y}\mathrm{\right|}\mathrm{ℛ}\mathit{,\sigma }\mathrm{)}\mathrm{=}\mathit{y}{\int }_{\mathit{y}}^{\mathrm{\infty }}\frac{\mathrm{𝒩}\mathrm{\left(}\mathit{x}\mathrm{\right|}\mathrm{ℛ}\mathit{,\sigma }\mathrm{)}}{\mathit{x}\mathrm{(}{\mathit{x}}^{\mathrm{2}}\mathrm{-}{\mathit{y}}^{\mathrm{2}}{\mathrm{)}}^{\mathrm{1}\mathit{/}\mathrm{2}}}\mathrm{d}\mathit{x,}$(4)where x ≡ R, and is the normal distribution with expected value ℛ and standard deviation σ; the first two moments of the distributions are related by

$\begin{array}{ccc}& & \mathit{x̅}\mathrm{=}\frac{\mathrm{4}}{\mathit{\pi }}\mathit{y̅}\\ & & \overline{{\mathit{x}}^{\mathrm{2}}}\mathrm{=}\frac{\mathrm{3}}{\mathrm{2}}\overline{{\mathit{y}}^{\mathrm{2}}}\mathit{.}\\ & & \end{array}$(5)As there is no known analytical solution to the integral in Eq. (4), this is evaluated numerically. From Eq. (4)the cumulative distribution function (c.d.f.) $\mathrm{\Phi }\mathrm{\left(}\mathit{y}\mathrm{\right|}\mathrm{ℛ}\mathit{,\sigma }\mathrm{)}\mathrm{=}{\int }_{\mathrm{0}}^{\mathit{y}}\mathit{\phi }\mathrm{(}{\mathit{y}}^{\mathrm{\prime }}\mathrm{\left|}\mathrm{ℛ}\mathit{,\sigma }\mathrm{\right)}\mathrm{d}\mathit{y}\prime$(6)is also estimated numerically. Hereafter the dependence on the parameters ℛ and σ is considered implicitly, that is, φ(y) = φ(y | ℛ,σ) and Φ(y) = Φ(y | ℛ,σ).

3.4. Censoring and truncation

The finite wavelength resolution of spectrographs and the intrinsic broadening of a non-rotating star set limits on the capability of measuring vsini below a certain threshold. This can be due to low v_eq, low i, or both. When i is sufficiently low, we can also expect that P cannot be measured because the star is seen almost pole-on and the rotational modulation induced by surface inhomogeneities is therefore undetectable. In practice, P is still measurable also when only a vsini upper limit can be estimated, and therefore the limits on vsini in general dominate those on P. Other cases when P is expected not to be measurable include uniformly distributed surface inhomogeneities and unfavourable photometric sampling.

A survival analysis (Klein & Moeschberger 2003; Feigelson & Jogesh Babu 2012) can be applied to recover ℛ and σ from a { Rsini } set in presence of censoring and truncation, in which cases Eq. (5)become invalid. When only a vsini upper limit, (vsini)_lim, is available for the object, this is translated into an Rsini upper limit using the relation $\mathrm{(}\mathit{R}\sin \mathit{i}{\mathrm{)}}_{\lim }\mathrm{=}\mathrm{𝒞}\mathit{P}\mathrm{(}\mathit{v}\sin \mathit{i}{\mathrm{)}}_{\lim }\mathit{,}$(7)which is then considered as a left-censored data point.

At the other extreme, ultra-fast rotators may have such a high vsini value that its measurement is very uncertain or impossible, and therefore only a lower limit is available. These measurements can be treated as right-censored data points, for which a (Rsini) lower limit can be defined in analogy with Eq. (7).

Left- and right-censored data points are taken into account through a survival function (Klein & Moeschberger 2003), which gives the probability that an object has a value above some specified level. For the case at hand, the survival function is $\mathit{S}\mathrm{\left(}\mathit{y}\mathrm{\right)}\mathrm{=}\mathit{P}\mathrm{(}\mathit{R}\sin \mathit{i}\mathit{>}\mathit{y}\mathrm{)}\mathrm{=}\mathrm{1}\mathrm{-}\mathrm{\Phi }\mathrm{\left(}\mathit{y}\mathrm{\right)}\mathit{,}$(8)where Φ(y) is given by Eq. (6). Using Eq. (8), we assign a probability S(y) if y is a right-censored (lower limit) data point and (1−S(y)) if y is a left-censored (upper limit) data point.

The incompleteness of the sample does not represent a limitation for the analysis as long as selection effects do not depend on the variates of interest. We expect, however, that when i is sufficiently low, neither vsini nor P can be measured and that the data are left-truncated in most cases. Data points like this are not present in the { Rsini } dataset, but we do not know the exact truncation value (Rsini)_trunc in each bin a priori. However, for a random orientation of the spin axis the probability of occurrence of the inclination i between i and i + di is sini di, and therefore the distribution favours high values of i. The results are therefore quite insensitive to the accuracy in determining the left-truncation level, and we adopt the practical approach of considering the lowest Rsini value as an approximation of (Rsini)_trunc in each bin. The way in which truncation is taken into account in the analysis is described in Sect. 3.5.

Right-truncation is ignored because we assumed that for ultra-fast rotator for which the vsini measurement is uncertain or impossible a lower limit is reported instead of omitting the measurement.

3.5. Estimating the mean radius from the projected radius distribution

After the p.d.f. and the corresponding c.d.f. and survival function (Eqs. (4), (6), and (8)) were defined, we defined a likelihood function as $\mathrm{ℒ}\mathrm{=}\underset{\det }{\mathrm{􏽙}}\mathit{\phi }\mathrm{\left(}{\mathit{y}}_{\mathit{j}}\mathrm{\right)}\underset{\mathrm{lcens}}{\mathrm{􏽙}}\mathrm{[}\mathrm{1}\mathrm{-}\mathit{S}\mathrm{(}{\mathit{y}}_{\mathit{j}}{\mathrm{\left)}}^{\mathrm{\right]}}\underset{\mathrm{rcens}}{\mathrm{􏽙}}\mathit{S}\mathrm{\left(}{\mathit{y}}_{\mathit{j}}\mathrm{\right)}\mathit{,}$(9)where the products are over detected, left-censored, and right-censored data points, respectively. Since our dataset is always left-truncated, in Eq. (9)we replace φ(y_j) by φ(y_j) /S(Y_L) and S(y_j) by S(y_j) /S(Y_L) (Klein & Moeschberger 2003), with Y_L = (Rsini)_trunc the left-truncation y value in the dataset. The negative log-likelihood function is therefore $\begin{array}{ccc}\hspace{0.17em}\begin{array}{c}\begin{array}{c}\end{array}\mathrm{-}\ln \mathrm{ℒ}\mathrm{=}\mathrm{-}\left(\sum _{\det }\ln \mathrm{\left(}\mathit{\phi }\mathrm{\right(}{\mathit{y}}_{\mathit{j}}\mathrm{\left)}\mathrm{\right)}\mathrm{+}\sum _{\mathrm{lcens}}\ln \mathrm{(}\mathrm{1}\mathrm{-}\mathit{S}\mathrm{(}{\mathit{y}}_{\mathit{j}}\mathrm{\left)}\mathrm{\right)}\mathrm{+}\sum _{\mathrm{rcens}}\ln \mathrm{\left(}\mathit{S}\mathrm{\right(}{\mathit{y}}_{\mathit{j}}\mathrm{\left)}\mathrm{\right)}\right)\mathit{.}\end{array}& & \end{array}$(10)The parameters ℛ and σ are finally evaluated by minimising the negative log-likelihood function Eq. (10) using the L-BFGS-B optimisation method of Byrd et al. (1995) as implemented in R.

3.6. Simulations

The method was extensively tested using numerical simulations. For brevity, only a summary of these tests is reported here. More details are reported in Appendix A.

For a sufficiently low level of censoring and truncation, that is, (Rsini)_trunc ≲ 0.1ℛ and (Rsini)_cens ≲ 0.3ℛ, and σ/ ℛ ≈ 0.1 the method is capable of reproducing ℛ within ≈2 percent and σ within ≈2 percent with n ≳ 10 data points (median over 100 realisation – see Appendix A).

Maintaining a low level of censoring and truncation, the accuracy degrades for extremely high or low values of σ/ ℛ. For σ → 0, the numerator in Eq. (4)tends to a Dirac delta function and φ develops a singularity at y = ℛ (see Chandrasekhar & Münch 1950). As a consequence, for σ/ ℛ ≲ 0.01 the numerator of the integrand in Eq. (4)is a very steep function of x, which causes numerical instabilities in the evaluation of the integral. Simulations for σ/ ℛ → 0.01 show that in this limit ℛ can still be recovered within ≈5 percent, while σ can be overestimated by a factor of several because of the smoothing implied by the integral in Eq. (4). At the other extreme, for σ/ ℛ ≳ 0.3 the theoretical distribution would imply a non-negligible probability of having negative (unphysical) values and the numerical procedure fails. Simulations for σ/ ℛ → 0.3 show that in such extreme cases ℛ can be recovered within ≈4 percent and σ within ≈50 percent (median over 100 realisations).

Fixing σ/ ℛ ≈ 0.1, an increase in the level of censoring and truncation does not significantly decrease the accuracy in reproducing ℛ and σ as long as the core of the distribution remains unaffected. In practice, both ℛ and σ are recovered within a few percent up to (Rsini)_trunc and (Rsini)_cens ≈ 0.6ℛ. Above this level, censoring and truncation affect the core of the distribution, making it difficult to recover its original shape.

We note that these simulations outline limitations that are intrinsic to the problem at hand and cannot be overcome using a simplified approach, and they set the boundaries for a meaningful estimate of the mean radius from an { Rsini } set.

4. Mean stellar radii in the Pleiades

Fig. 3 Comparison of the expected radius ℛ with the models of Baraffe et al. (2015; dotted line) and Spada et al. (2013; dashed line). ℛ is reported as a function of the average MK in each bin. The three values of ℛ correspond to our best estimate of the SDR (blue diamond) and the estimated upper and lower SDR limits (black and red diamonds). The horizontal bars encompass the MK range in each bin. The vertical bars represent the standard deviation σ for each bin centred on our best estimate of the SDR. Right panels show the ratio with the Baraffe et al. (2015) model. Panels a) and a′) are obtained using all data, panels b) and b′) for the slow-rotator sequence, panels c) and c′) for stars with periods between fast and slow rotators, panels d) and d′) for fast rotators only (P< 2 d). Bins are ΔM ≈ 0.1 M⊙ wide except for the fast rotators set, for which ΔM ≈ 0.2 M⊙ to ensure that at least ≈10 stars fall in each bin.

The method described in Sect. 3 was applied for the first time to the dataset available for the Pleiades (Sect. 2). The data were organised in M_K bins corresponding to fixed ΔM according to the theoretical relationship reported by Baraffe et al. (2015) for M_K vs. M at the Pleiades age. For each bin we compute ℛ and σ from the Rsini distribution. In Fig. 3 we compare ℛ vs. ⟨M_K⟩ with the models of Baraffe et al. (2015) and Spada et al. (2013). The comparison is carried out for the whole P range, for the slow-rotator sequence (as identified in Lanzafame & Spada 2015), for the fast rotators (P< 2 d), and for the stars with periods in between the fast rotators and the slow-rotator sequence. According to the definition of Barnes (2003), these correspond to the I sequence, the C sequence and stars in the gap². Bins are ΔM ≈ 0.1 M_⊙ wide, except for fast rotators, for which ΔM ≈ 0.2 M_⊙ to have at least ten stars in each bin, and which are allowed to overlap to provide ℛ estimates at steps of ≈0.05 M_⊙. For 1.2 ≳ M/M_⊙ ≳ 0.85, using the whole P range, ℛ agrees with the theoretical models within ≈10 percent (≈5 percent for 1.2 ≳ M/M_⊙ ≳ 1.00) with no significant bias. The noise in the data, however, increases with decreasing mass, resulting in increasingly larger σ. For 0.85 ≳ M/M_⊙ ≳ 0.65, ℛ becomes systematically larger than theory, although still within ≈1σ. Restricting the calculations to stars belonging to the slow-rotator sequence, ℛ still agrees with the theoretical models within ≈10 percent with no significant bias. The results for the fast rotators have larger σ and ℛ scatter, up to ≈15 percent, mostly because larger bins are required to have at least ten data points in each bin, but no systematic deviations from the models are found. The discrepancy with the theoretical model is confined to stars with mass just below the low-mass end of the slow-rotator sequence and P> 2 d (panels (c) and (c′) in Fig. 3), for which ℛ is inflated at 2σ level. According to the scenario described by Barnes (2003), these are stars that converge on the slow-rotator sequence (gap).

Possible alternative explanations for this behaviour are

• 1.

  systematic deviations due to some outliers;

• 2.

  P-dependent observational biases not taken into account by our procedure;

• 3.

  biases in the vsini datasets.

To consider the first alternative explanation, we repeated the calculation using different binning and obtained essentially the same results. We also repeated the calculations excluding m = 3 stars in each bin and using all n−m combinations, with n the number of stars in each bin, with no significant change with respect to the original results. In this way we verified that there is indeed a group of measurements that produces the observed deviation and not just some isolated outliers, to which our procedure is rather insensitive in any case.

The second alternative explanation is deemed rather unlikely since the decrease in fraction of stars with both P and vsini measurements is rather uniform in P down to at least M ≈ 0.6 M_⊙ (see Fig. 1) and the P dataset of Hartman et al. (2010) is 93 percent complete in the mass range 1.0 ≳ M/M_⊙ ≳ 0.7.

The third alternative explanation is therefore the only one that could be of some concern. Hartman et al. (2010) estimated that for M ≲ 0.85 M_⊙ the vsini dataset could be affected by a bias Δvsini ≈ −1.5 km s^-1, but this was still insufficient to explain the sini distribution, so that the authors invoked other factors like the radius inflation and a rather high SDR. By constraining the SDR effects as described in Sect. 3.2, we estimate that to explain the ℛ vs. ⟨M_K⟩ discrepancy at M ≈ 0.7 M_⊙, a bias of at least Δvsini ≈ −4.0 km s^-1 limited to a rather restricted M_K range would be necessary. Considering the expected uncertainties in the intrinsic width of non-rotating stars (e.g. Queloz et al. 1998), which are expected to be the main cause of such systematic deviations, this bias seems too high and its dependence on spectral type much steeper than deemed plausible.

As discussed in Sect. 2, vsini of fast rotators are mostly adopted from Soderblom et al. (1993), while for the other stars they are mostly taken from Queloz et al. (1998) for M ≳ 0.85 M_⊙ and from Terndrup et al. (2000) for M< 0.85 M_⊙. Excluding the fast-rotator measurements of Soderblom et al. (1993), the comparison for the five stars in common between the two remaining datasets does not point to any significant bias. To investigate the possible vsini bias in more detail, we repeated the calculations on a ΔM ≈ 0.2 M_⊙ bin centred on M ≈ 0.7 M_⊙ using all available vsini, only the Q98, and only the T00 vsini. The results, reported in Table 2, show that in all three cases the expected mean radius is larger than the theoretical mean radius by 1 or 2σ. We note that the discrepancies are smaller than those shown in Fig. 3 (panels c) and c′)) because to have at least ten stars for each set, the bin is larger and the mean mass in two sets (Q98 and T00) differs by 0.06 M_⊙. In conclusion, we have no evidence of an observational vsini bias at the level required to explain the radius discrepancy found for stars that are converging on the slow-rotator sequence.

From a different perspective, it can be argued that the results obtained for the fast-rotator sequence may be affected by the possible omission of high vsini values that have not been reported as lower limits as our method requires. This aspect can be of concern as more than one hundred stars of the Hartman et al. (2010) periods sample have no vsini measurement and a significant fraction of them have short periods. In the Pleiades case, however, we can reasonably assume that the lacking vsini data do not depend on its value as high values of vsini and lower limit are reported. The treatment of datasets in which high vsini values are lacking would require sufficient information to allow considering right-truncation in Eq. (9).

5. Conclusions

We have set up a new method for deriving mean stellar radii from rotational periods, P, and projected rotational velocities, vsini, based on the survival analysis concept (Klein & Moeschberger 2003). This method exploits the whole information content of the dataset with an appropriate statistical treatment of censored and truncated data. Provided censoring and truncation do not significantly affect the peak of the Rsini distribution and that there is no significant bias in the data, the method can recover the mean stellar radius with an accuracy of a few percent with as few as n ≈ 10 measurements. The total standard deviation, σ, which cumulatively takes the data noise and the intrinsic R standard deviation into account, can also be estimated with an accuracy of a few percent except in extreme cases where the distribution is too broad (σ/ ℛ ~ 0.3) or too narrow (σ/ ℛ ~ 0.01).

The method has been applied for the first time to the dataset available for the Pleiades. We found that deviations of the empirical ℛ vs. M_K relationship from standard models (e.g. Spada et al. 2013; Baraffe et al. 2015) do not exceed 5 percent for 1.2 ≥ M/M_⊙ ≥ 1.0 and 10 percent for 1.0 > M/M_⊙ ≥ 0.85, with no significant bias. Evidence of a systematic deviation at 1−2σ level of the empirical ℛ vs. M_K relationship from standard models is found only for stars with M ≈ 0.7 ± 0.1 M_⊙ that are converging on the slow-rotator sequence. Deviations of the ℛ vs. M_K relationship for fast rotators (P < 2 d) do not exceed ≈15 percent in the whole mass range with no evidence of a systematic deviation from standard models. No evidence of a radius inflation of fast rotators in the Pleiades is therefore found.

Acknowledgments

F.S. acknowledges support from the Leibniz Institute for Astrophysics Potsdam (AIP) through the Karl Schwarzschild Postdoctoral Fellowship. Research at the Università di Catania and at INAF – Osservatorio Astrofisico di Catania is funded by MIUR (Italian Ministry of University and Research). The authors warmly thank Sydney Barnes (Leibniz Institute for Astrophysics Potsdam) for valuable discussions and an anonymous referee for useful comments. This research made use of the NASA Astrophysics Data System Bibliographic Services (ADS) and the R software environment and packages (https://www.r-project.org).

References

Appendix A: Accuracy of the method

In this appendix we present some examples and tests carried out to estimate the method accuracy.

Fig. A.1 Illustrative tests in recovering ℛ from a synthetic {Rsini} set with censoring. The simulations reproduce n = 30 slow-rotator sequence stars in the Pleiades with M ~ 1 M⊙ and R ~ 1 R⊙ at varying censoring levels: (vsini)lim = 0 (panel a)); 6 (panel b)); 9 (panel c)); 10 (panel d)) km s-1. Blue and red histograms are used for censored and not-censored data. The kernel density estimate is shown as a thin blue line for not-censored data and as a red thin line for censored data, the function φ corresponding to the given synthetic distribution parameters is plotted as a thick black line and that corresponding to the reconstructed parameter as a thick red line. Histograms are representations of frequency.

We considered as an illustrative example a group of stars with M ~ 1 M_⊙ in the Pleiades slow-rotator sequence (Lanzafame & Spada 2015). This group of stars has a normal distribution of periods with d and σ_P ≃ 0.35 (see Table 2 in Lanzafame & Spada 2015). Assuming R ~ 1 R_⊙, it follows that the equatorial velocities have a normal distribution with $\overline{{\mathit{v}}_{\mathrm{eq}}}\mathrm{\simeq }\mathrm{12.7}$ and σ_v ≃ 1.3 km s^-1. Assuming that this group is composed of n = 30 stars, we generated synthetic { Rsini } datasets with these parameters by applying different levels of censoring and truncation. Figure A.1 reports some examples in the censoring no-truncation cases. These tests show that for the Pleiades dataset considered here, where Rsini truncation does not exceed 0.4, the maximum censoring is around (vsini)_lim ~ 7 km s^-1, and σ/ ℛ ≈ 0.1, we expect that ℛ can be recovered with a precision better than 2% and σ with a precision better than 1%. We note that censoring at (vsini)_lim ≳ 9 km s^-1 affects the core of the distribution, making it difficult to recover its original shape.

To evaluate how the accuracy depends on the number of measurements and on the σ/ ℛ ratio, we applied the method to groups of 100 Rsini synthetic distribution realisations, each one with a different number of stars in each bin and different σ/ ℛ values. For brevity, here we compare the results obtained with no-censoring and no-truncation with the worse levels of censoring and truncation in the Pleiades dataset. Furthermore, we compare the results obtained for σ/ ℛ = 0.1, which is representative of non-extreme values of the radius dispersion, with those obtained with σ/ ℛ = 0.03 and 0.3, which are representative of extremely low and extremely high values of the radius dispersion. We recall that σ takes both the intrinsic radius dispersion and the observational uncertainties into account.

Fig. A.2 Box-and-whisker plots for the reconstructed ℛ and σ from 100 random realisations of the Rsini distribution with no censoring and no truncation at different σ/ ℛ values and different number of measurements per bin.

Figure A.2 shows the results for no censoring and no truncation. For σ/ ℛ = 0.1 the median of the ratio of the expected value to the true value of ℛ ranges from 0.976 for n = 10 to 0.993 for n = 100. For narrow distributions (e.g. σ/ ℛ = 0.03) this is essentially unity with n = 15 or more while for extremely broad distributions (e.g. σ/ ℛ = 0.3) it amounts to 0.963 in the worst case. The scatter in the expected-to-true ℛ ratio decreases with decreasing width of the distribution and with increasing number of observations. The median of the reconstructed σ behaves in a similar way, although its relative accuracy decreases more significantly with increasing width of the distribution. The intrinsic skewness of the Rsini distribution leads to a general tendency of underestimating the true ℛ and σ for small n and large σ/ ℛ.

Fig. A.3 Box-and-whisker plots for the reconstructed ℛ and σ from 100 random realisations of the Rsini distribution with censoring and truncation at different σ/ ℛ values and different number of measurements per bin.

The results of the simulation with the worse level of censoring and truncation in the Pleiades dataset are summarised in Fig. A.3. For sufficiently narrow distributions the accuracy does not degrade significantly with respect to the no-censoring and no-truncation cases. Only for large σ/ ℛ, that is, when censoring and truncation affect the peak of the Rsini distribution, the median of the ratio of the expected value to true value of ℛ is significantly below unity. For the σ/ ℛ = 0.3 case shown in Fig. A.3, the radius is underestimated by ≈10 percent (median) even when the number of measurements is increased to n = 100. We note, however, that this latter condition is never met in the Pleiades dataset we analysed in this paper and it is presented here to outline a condition in which it is not possible to recover the average radius reliably.

In summary, the golden rule for evaluating the mean radius is that censoring and truncation must not affect the core of the Rsini distribution. We argue that this is a general requirement, which is not due to a limitation of this particular method, but to the lack of sufficient information when the data cannot define the core of the Rsini distribution with sufficient detail.

All Tables

All Figures

Fig. 1 Hartman et al. (2010)P dataset of the Pleiades (blue filled circles). Upper panel: open red squares outline targets for which vsini is available and green stars those considered in this work. Lower panel: open red triangles indicate vsini from Queloz et al. (1998), green squares the values adopted from Terndrup et al. (2000), plum pentagons those from Stauffer & Hartmann (1987), and olive hexagons those from Soderblom et al. (1993). See text for details.

In the text

	Fig. 2 Fit to the mean period vs. equatorial period for the AB Dor young loose association.
In the text

Fig. 3 Comparison of the expected radius ℛ with the models of Baraffe et al. (2015; dotted line) and Spada et al. (2013; dashed line). ℛ is reported as a function of the average MK in each bin. The three values of ℛ correspond to our best estimate of the SDR (blue diamond) and the estimated upper and lower SDR limits (black and red diamonds). The horizontal bars encompass the MK range in each bin. The vertical bars represent the standard deviation σ for each bin centred on our best estimate of the SDR. Right panels show the ratio with the Baraffe et al. (2015) model. Panels a) and a′) are obtained using all data, panels b) and b′) for the slow-rotator sequence, panels c) and c′) for stars with periods between fast and slow rotators, panels d) and d′) for fast rotators only (P< 2 d). Bins are ΔM ≈ 0.1 M⊙ wide except for the fast rotators set, for which ΔM ≈ 0.2 M⊙ to ensure that at least ≈10 stars fall in each bin.

In the text

Fig. A.1 Illustrative tests in recovering ℛ from a synthetic {Rsini} set with censoring. The simulations reproduce n = 30 slow-rotator sequence stars in the Pleiades with M ~ 1 M⊙ and R ~ 1 R⊙ at varying censoring levels: (vsini)lim = 0 (panel a)); 6 (panel b)); 9 (panel c)); 10 (panel d)) km s-1. Blue and red histograms are used for censored and not-censored data. The kernel density estimate is shown as a thin blue line for not-censored data and as a red thin line for censored data, the function φ corresponding to the given synthetic distribution parameters is plotted as a thick black line and that corresponding to the reconstructed parameter as a thick red line. Histograms are representations of frequency.

In the text

	Fig. A.2 Box-and-whisker plots for the reconstructed ℛ and σ from 100 random realisations of the Rsini distribution with no censoring and no truncation at different σ/ ℛ values and different number of measurements per bin.
In the text

	Fig. A.3 Box-and-whisker plots for the reconstructed ℛ and σ from 100 random realisations of the Rsini distribution with censoring and truncation at different σ/ ℛ values and different number of measurements per bin.
In the text