© ESO 2018

1 Introduction

Gaia all-sky scanning and multi-epoch observations offer unprecedented opportunities to detect and characterise stellar variability. One of the variability phenomena that can be studied using Gaia data is the flux modulation induced by surface inhomogeneities (such as spots and faculae) combined with rotation. An analysis of this modulation allows us to infer the stellar rotation period. In this paper we refer to dwarf, sub-giant, and T Tauri stars showing such flux modulation as BY Dra variables, after the archetype BY Draconis.

As for the Sun, the variability in late-type stars also includes more irregular effects that are due to the emission originating from the chromosphere/corona, and to sporadic outbursts that in turn are due to reconnection of magnetic flux tubes, that is, flares. Numerous classes of stars (single and multiple systems) are affected by this variability, including BY Dra, UV Cet, and RS CVn types. As in classical papers on the subject, this classification is phenomenological: when flares are observed, the star is classified to be of UV Cet type; a BY Dra can be a single star or a star in a binary system; and when the star is in a close binary system, a condition that enhances magnetic activity, the system is classified to be of RS CVn type.

Variability due to spots and faculae can reach amplitudes of the order of some tenths of magnitudes for the most active objects (e.g. T Tauri stars), but can be of the order of milli-magnitudes or lower for late-type dwarfs of solar age. Inferring rotation periods and modulation amplitude in a large and heterogeneous sample of stars, both in clusters and in the field, for which other astrophysical parameters can also be inferred homogeneously, as in the Gaia case, opens new frontiers in our understanding of stellar angular momentum evolution. Together with information from asteroseismology, the Gaia data and analysis provide previously unattainable insights that enable us to effectively study issues like the dependence of surface rotation on age, the existence of different dynamo regimes, and the role of multiplicity, planetary systems, and environment at large in the stellar angular momentum evolution.

As part of the activities carried out by the variability Coordination Unit (CU7) of the Gaia Data Processing and Analysis Consortium (DPAC), we have devised and implemented methods for the identification and analysis of rotational modulation in low-mass stars. In this paper we present these methods and give an overview of the periods and modulation amplitudes obtained on the first 22-month Gaia baseline, which are available in Gaia DR2. In Sect. 2 we present the methods we developed, in Sect. 3 we describe the results, and in Sect. 4 we draw our conclusions.

2 Method

The problem of inferring stellar rotation periods from Gaia photometric time series has been discussed in Distefano et al. (2012). In brief, on the one hand, the Gaia sampling depends on the star’s position on the sky and is irregular, in the sense that the time interval between two consecutive transits is not constant and rather sparse. There are, of course, regularities like those due to the consecutive passages of the star’s image on the two satellite’s FoVs and the precession of the satellite, but the time step is not constant. Gaia DR2 also contains a subset in which regions around the ecliptic poles were observed with a regular sampling. On the other hand, in analogy with the solar case, the distribution of inhomogeneities on the surface of a low-mass star surface changes in time as a result of the evolution of spots, faculae, and active region complexes, each with its own characteristic timescale, making the rotational modulation signal, in general, non-coherent over long timescales. Furthermore, the stochastic variability due to flaring is superimposed on the rotational modulation.

Distefano et al. (2012) adopted the detailed light-curve model of Lanza et al. (2006b) for their simulation of Gaia observations of solar twins. Extrapolating these models to stars of different mass and/or magnetic activity is still very uncertain and not justified for the Gaia case because of the sparseness of its sampling. We therefore do not attempt to fit a detailed model to the data, but rather search for periodicities in the light curves that can be produced by rotational modulation by assuming a sinusoidal model. Since this can be applied to any (quasi-)periodic variability, it is necessary to constrain our classification using other information available in the Gaia data.

To achieve a reliable identification of stellar variability classes like that of BY Dra, we plan to combine the variability parameters with the Gaia inference of astrophysical parameters (Bailer-Jones et al. 2013). However, only preliminary astrophysical parameters based on integrated photometry are available in DR2 (Andrae et al. 2018), and therefore, we here adopted a simplified approach (Sect. 2.1).

Within the Gaia variability pipeline, rotational modulation is processed by two modules: the special variability detection (SVD) solar-like, devised for a first selection of candidates, time-series segmentation, and outlier detection (possible flares), and the specific object study (SOS) rotational modulation, devised for period search, modulation amplitude estimate, and filtering of the final results. The whole procedure is described in the following.

2.1 Input data

We refer to Holl et al. (2018) for a summary of the general variability processing, analysis, and the observation filtering and flagging applicable to the results published in Gaia DR2.

The simulations presented in Distefano et al. (2012) indicate that the analysis of rotational modulation in solar-like stars requires at least a dozen samplings in at least one interval not longer than ~150 d (see also Sect. 2.2). The rotational modulation processing, therefore, follows the geq20 path described in Holl et al. (2018), that is, the path starting from the 826 million sources with ≥20 field-of-view(FoV) transits.

From the initial geq20 catalogue, an initial list of BY Dra candidates was compiled by selecting sources around the observed main sequence in the (M_G, (G_BP − G_RP)) diagram. The approximate absolute magnitude M_G was estimated from parallax and G, ignoring reddening, to build the diagrams shown in Fig. 1. Only sources with a relative parallax uncertainty better than 20% were considered. The selected region (coloured region in Fig. 1) broadly embraces the main sequence. The leftmost boundary, (G_BP − G_RP) = 0.6, limits the selection to stars of spectral type approximately later than F5.

The initial list of BY Dra candidates was then further limited to sources whose time series could be divided into at least two segments with a number of transits N_G ≥ 12. The segmentation algorithm is described in Sect. 2.2.

Fig. 1 Region of the HRD in which rotational modulation was found. The region is highlighted with colours (rainbow). The colour code indicates the relative density of data points (red for lower density, and purple for higher density). The HRD shown is built using the AGIS02.2 solution for the parallax (left panel, see Lindegren et al. 2018, for details) uncorrected for interstellar extinction.

2.2 Time-series segmentation

Lanza et al. (2004) analysed the total solar irradiance (TSI) time series collected by the VIRGO experiment on the SoHo satellite and found that period-search algorithms can reproduce the correct solar rotation period only if the effects due to the evolution of the active region complexes are limited by dividing the long-term time series into sub-series of ~150 d. The active region evolution timescale in young stars, which have a higher magnetic activity than the Sun, can be significantly shorter (~60–90 d, Messina et al. 2003).

These effects are illustrated in Fig. 2, where we show the expected G magnitude variations for a solar-like star with an effective temperature T_eff = 5000 K, a rotation period P = 5 d, and an average T_spot = 4200 K, according to the simulation code described in Distefano et al. (2012). The evolution of the surface magnetic field produces changes in area, temperature, and location of the active region, which gradually distort the signal produced by the rotational modulation.

The segmentation of long-term time series also allows detecting period variations induced by the spot latitude migration combined with the stellar surface differential rotation (SDR). When a star is characterised by SDR and by a solar-like cycle, where the spots gradually migrate from intermediate latitudes to the equator, the detected rotation period changes slightly according to the cycle phase at which the star is observed (see e.g. Distefano et al. 2016).

The advantage of working with a more coherent pattern of flux variations in segments than in the whole time series comes at a price, as the segments contain fewer samplings and the period extracted from them may consequently have a lower statistical significance. Therefore, a compromise has to be found between the length of the segments and the number of points in the segments. Distefano et al. (2016), for instance, segmented the long-term photometric time series of a sample of young solar-like stars into 50 d long sub-series by means of a sliding-window algorithm. Unfortunately, except for the observations of the ecliptic poles, this cannot be applied in general to the Gaia time series as many segments would have only two samplings ≃106 min apart, corresponding to two consecutive transits on the two telescope’s FoVs (see e.g. Eyer & Mignard 2005; Distefano et al. 2012). In the nominal scanning law (NSL) observations,¹ segments 50 d long would comprise ~10–45 samplings only at ecliptic latitudes close to β = ±45°.

We therefore modified the Distefano et al. (2016) procedure into an adaptive segmentation algorithm that takes the peculiar features of the Gaia NSL into account. In this new segmentation algorithm, the length of the sliding window is not fixed, but changes according to the sampling of each time series. The algorithm first searches groups of consecutive samplings separated by less than 6 h. Each group is identified by the set , the observation time of the first point, the observation time of the last point, and the total number of observations per group. Then all possible segments with extremities and with j ≥ i, that is, all segments starting from the first point of a group to the final point of the same or any other subsequent group, are considered, and those satisfying the conditions (1) (2)

are selected for subsequent processing. Nested segments are discarded, in the sense that if one segment satisfying the above conditions is entirely included in a longer segment, only the longer segment is considered. In this way, the window length is set so that each segment covers an interval not exceeding 120 d and includes at least 12 samplings.

The upper limit of the segment length (Eq. (1)) may in some cases hamper a period detection in young stars, but this limit arises from an unavoidable compromise between the typical timescale of active region evolution and the characteristics of the Gaia sampling. At any rate, possible spurious periods that may results from active region evolution with timescales shorter than the segment length are filtered out, as discussed in Sect. 2.6.

Fig. 2 Simulated light curve for a solar-like star with an effective temperature Teff = 5000 K, a rotation period P = 5 d, and an average Tspot = 4200 K.

2.3 Magnitude–colour variation correlation and outliers

Magnitude and colour variations induced by rotation are linearly correlated, insofar as surface inhomogeneities maintain approximately the same location, area, and temperature, and their surface distribution is not uniform. This is confirmed by observations(Messina et al. 2006) and can be adopted as a sufficient, but not necessary, condition for rotational modulation detection in late-type stars. The evolution of active regions may gradually change the characteristics of such a correlation so that when an isolated sampling is superimposed on a group of close samplings at a distant epoch, it may appear as an outlier. More significantly, flares may also appear as outliers in such a correlation, since they can be viewed as sudden changes of the surface inhomogeneity configuration produced by magnetic energy dissipation with consequent transient plasma heating.

Outliers are identified using the Theil–Sen robust linear regression between G and (G_BP − G_RP) (see Gilbert 1987; E-Shaarawi 2001, for details). Briefly, this algorithm computes the slopes and the intercepts of the straight lines passing through each couple of data points (G, (G_BP − G_RP)) measured inthe segment. The median slope and the corresponding intercept are taken as a first estimate of slope andintercept, and outliers are identified as points that are significantly distant from the straight line that best fits the data. Furthermore, since flares can produce significant blueward changes in (G_BP − G_RP) without greatly affecting theG magnitude,² we add transits with (3)

to the list of outliers. The identification of flare events amongst the outliers is discussed in Distefano et al. (in prep.). After the outliers are removed, the final slope and intercept are computed using a simple linear regression.

The Pearson correlation coefficient r between the G and (G_BP − G_RP) variations is also computed. The closer r is to ± 1, the higher the probability that rotational modulation is occurring, but period search (Sect. 2.4) is carried out in all cases.

2.4 Rotation period

Period search is performed using the generalised Lomb-Scargle periodogram method as implemented by Zechmeister & Kürster (2009). This was previously adopted in the test analysis of Gaia DR1 data (Eyer et al. 2017) and is best suited for the Gaia case (Distefano et al. 2012).

The search is carried out on each segment and on the whole 22-month baseline. Calling T the time interval spanned by the segment (or by the whole time series), the frequency range in which the search is performed is (4)

with a frequency step d⁻¹. The upper limit 3.2 d⁻¹ is chosen in order to avoid the aliases corresponding to the 6-h rotation period of the satellite.

The period with the highest power in the periodogram is selected, and the false-alarm probability (FAP), indicating the probability that the detected period is due to just noise, is associated with it. In principle, methods based on Monte Carlo simulations would give the most reliable estimate of the FAP associated with a given period (see e.g. Distefano et al. 2012), but they are computationally expensive and prohibitive given the size of the Gaia sample. Süveges et al. (2015) demonstrated that the Baluev formulation (Baluev 2008) is computational inexpensive and quite effective also in presence of strong aliases, and therefore this is adopted for the case at hand.

The choice of the FAP threshold has to balance the rates of false positives and true negatives. For DR2, valid periods are chosen as those with FAP ≤ 0.05, which corresponds to a false-detection rate of about 5% (Süveges et al. 2015).³

When a significant period P is detected, a sinusoidal model (5)

is fitted to the time-series segment using the Levenberg–Marquardt method (Levenberg 1944; Marquardt 1963).

Periods inferred in different segments of the same source are finally combined using the mode as statistical estimator. For each source, a double loop is performed on the detected periods, counting for each period P_i, the number N_i of periods that differ by less than 20% (which includes both uncertainties and the effects of differential rotation, see e.g. Distefano et al. 2016) from P_i. We then consider the groups of similar periods (i.e. periods similar within the same group and different from one group to another considering the 20% tolerance) with the highest frequency N_max and select the group with the lowest ⟨FAP⟩. The final period is given as the average of the periods in the selected group P = ⟨P_i⟩.

2.5 Magnetic activity

The amplitude of the rotational modulation is related to the non-axisymmetric part of the active region distribution, which is a complex function of the projected area and the contrast in temperature between active regions and the unspotted stellar surface (see e.g. Lanza et al. 2006a; Distefano et al. 2012). The area and temperature of active regions, in turn, depend on the strength and topology of the surface magnetic field, which makes the amplitude of the rotational modulation a useful proxy for studyingthe stellar magnetic activity (see e.g. Rodonò et al. 2000; Ferreira Lopes et al. 2015; Lehtinen et al. 2016).

Given the sparseness of the Gaia sampling and the presence of outliers (flares), we derive two quantities related to the rotational modulation amplitude, to be compared amongst each other for filtering the final results. These are the differences between the 95th (G_95th) and 5th (G_5th) percentiles of the G magnitudes measured in the segment, (6)

and the amplitude derived from the sinusoidal fitting of the light curve (Eq. (5)) (7)

Both quantities give a measurement of the modulation amplitude that takes the presence of outliers into account, but they can be quite different when the light curve is not sampled with sufficient uniformity or when the sinusoidal model is a very poor representation of the data. This is used in the final candidate filtering as described in Sect. 2.6.

In DR2, an estimate of the unspotted magnitude is given as the brightest magnitude G observed in each segment. A more robust estimate, however, can be obtained from the Eq. (5) fit coefficients as (8)

2.6 Final filtering

The algorithm described above found periodicity in some 7 × 10⁵ sources in the Hertzsprung–Russel diagram (HRD) region described in Sect. 2.1. These may include spurious detections derived from the sampling incompleteness as well as other classes of periodic variability; eclipsing binaries and ellipsoidal variables are the only other expected classes. Limiting the analysis to objects with a relative parallax uncertainty better than 20% (Sect. 2.1) makes the expected contamination from pulsating variables negligible, as only outlier in parallax would fall in the selected sample.⁴ To minimisecontamination as much as possible, we impose constraints on phase and amplitude in the period folded light curves without assuming a strictly sinusoidal variation.

The constraints on phase are based on two parameters: the phase coverage (PC), and the maximum phase gap (MPG). Data in a given segment are folded using the period detected in that segment and binned in equally spaced phase intervals.⁵ The number of bins containing at least one data point is then divided by the total number of bins to obtain the phase coverage, which is therefore a real number in the [0, 1] interval. A uniform phase coverage produces PC → 1, while a very concentrated phase coverage gives PC → 0. Accepting only sources with PC ≃ 1 would leave out many sources with reliable period estimates, but extending the accepted sample to PC ≳ 0.5, for instance, would include some extreme cases with unacceptably large phase gaps, as the sampled phases can be concentrated in only half of the period-folded time series, for example. In order to balance these two aspects, the MPG parameter, (9)

where (10)

with ϕ_i and ϕ_j the phases of the i^th and j^th phase in theperiod-folded segment, is also taken into account. In DR2, accepted sources are those satisfying the two requirements (11) (12)

in at least one segment, which ensures a good phase coverage with acceptable gaps.

The constraints on amplitude are based on the (A_fit∕A_per) ratio and on the expected maximum amplitude. In order to filter out light curves that are far from sinusoidal and/or have poor sampling, the (A_fit ∕A_per) distribution is cut at the 5^th and 95^th percentile, which corresponds at selecting sources with (13)

only. Regarding the expected amplitudes, the analysis of ≈30 000 main-sequence stars observed by Kepler (McQuillan et al. 2014) places them in the 0.006–0.025 mag range (white light). Messina et al. (2011) found a ΔV ≈ 0.02−0.5 rotational modulation range in stars belonging to the young loose associations ϵ Cha and η Cha. Strassmeier et al. (1997) reported an amplitude ΔV = 0.65 mag for the T Tauri star V410 Tau and an amplitude ΔV = 0.48 mag for the giantstar XX Tri. Based on previous knowledge, therefore, we conservatively filter out light curves with amplitude larger than 1 mag, considering it unlikely that such large variations can be due to just rotational modulation induced by surface inhomogeneities.

A further external filter is applied a posteriori by cross-matching the final BY Dra sample with all other Gaia variability classes whose analysis can exploit more regular and stable light curves, and therefore their classification is considered more reliable than that of BY Dra. The comparison with the still unpublished eclipsing binaries is particularly interesting because it gives an order of magnitude of the remaining eclipsing binaries contamination after filtering. Such a comparison resulted in ≈ 100 sources in common between the BY Dra and the eclipsing binary classifications.

An illustrative example of a final BY Dra candidate is shown in Fig. 3 (Gaia DR2 source id 894812322514469120). In this case, periods of ≈3.9−4.6 d are found in three segments with A_fit∕A_per ≈ 1 and a quasi-sinusoidal shape in all of them.

For each source that passed the final filtering, the whole solution is contained in Gaia DR2. This includes the source’s time-series, the start and end observation times defining each segment, the solution found in each segment (P, A_per, and the coefficients in Eq. (5)) and the source’s final period and amplitude derived from all segments.

Fig. 3 Period-folded light curves for the Gaia DR2 source id 894812322514469120. Three similar significantperiods were found in two segments and in the whole time series with similar Afit and Aper.

3 Results

3.1 Coverage

The analysis of the first 22 months of Gaia observations with the method and filtering described in Sect. 2 resulted in the identification of 147 535 BY Dra variable candidates whose output parameters are included in Gaia DR2. This sample covers 38% of the whole sky when divided into bins of ≈0.84 square degrees (level 6 HEALPix, see Górski et al. 2005). Figure 4 shows the sky distribution of period measurements. The distribution is restricted to regions that are scanned more frequently in the 22-month baseline. The filtering described in Sect. 2 discards poorly sampled period-folded time series, so that the final sample has a range of 18 to 236 observations, with a median of 43, a mean of 45.65, a 5^th percentile of 28, and a 95^th percentile of72 observations per target.

As anticipated by Distefano et al. (2012), the Gaia frequency sensitivity favours the detection of shorter periods (≲ 10 d), as shown in Fig. 5. As expected, detection of longer periods is concentrated around ecliptic latitudes ± 45°. Figure 5 also shows the presence of some moderate aliases, the more perceivable lie at ≈23.5, 25.3, and 31.6 d. We aim at some deeper understanding of these features in future analyses.

Fig. 4 Sky distribution of Gaia DR2 BY Dra candidates. P is given in day units.

Fig. 5 Rotation period distribution of Gaia DR2 data for BY Dra stars (P is given in day units). A 3 d binning P histogram is shown in the left panel, together with the Gaussian kernel density estimate with a 3 d width and 3-sigma truncation. The P – ecliptic latitude scatter plot (right panel) shows the higher number of period measurements at ± 45° ecliptic latitudes, particularly regarding longer periods. The colour code (rainbow) is an indication of the relative density of data points (red for lower density, and blue for higher density). Both panels give an indication of the moderate degree of aliasing, the most evident being at ≈23.5, 25.3, and 31.6 d.

3.2 Completeness

Estimating the completeness of the BY Dra sample in Gaia DR2 is hampered by the lack of a well-established model for the expected number in the Galaxy. It is expected that all stars with a convective envelope develop surface inhomogeneities, although in some cases these may not produce an observable rotational modulation. Detections of new low-mass dwarfs displaying rotational modulation is continuously increasing with the ever larger span and increasing sensitivity of modern surveys, none of which, however, has the full sky coverage capabilities of Gaia.

At this stage, it is possible to perform some meaningful comparison only with the HATNet (Hartman et al. 2010) and K2 (Rebull et al. 2016) surveys of the Pleiades, with which, nevertheless, the DR2 geq20 sky coverage still overlap only marginally (see Appendix A for details). The Hartman et al. (2010) catalogue, in particular, is useful for a comparison with the Gaia data because of its uniform coverage of an approximately square sky region including the Pleiades cluster and an area close to it. Furthermore, even though the Kepler (K2) pixel size is smaller than that ofHATNet and Kepler can observe fainter stars than HATNet, the Hartman et al. (2010) catalogue contains a larger number of non-members in a smaller area, making our rough estimate easier and more realistic (see Appendix A). Assuming that the Hartman et al. (2010) catalogue lists all the BY Dra in its FoV down to G ≈ 14.5, we estimatethat the completeness of the BY Dra sample is 14% in the overlapping field. Assuming that this value is uniform over the whole sky, we estimate a completeness upper limit of 5%. At the other extreme, we may assumethat all low-mass dwarfs are BY Dra variables. Then, comparing with all stars observed by Gaia in the same sky region and magnitude range, we estimate a lower completeness limit of 0.7%.

3.3 Contamination

An estimate of the contamination can be based on the type and amount of expected contaminants left after our pre-selection, based on stars around the main sequence, and our filtering (Sect. 2), which discards results with FAP > 0.05 and light curves that are far from being sinusoidal. Neglecting outliers in parallax and the expected false-detection rate of about 5%, the only remaining expected main contaminants are close grazing binaries with a period ≲ 10 d, whose light curve is not far from being sinusoidal. When the binary period is long, then the grazing eclipse covers a limited phase interval, and it does not look like rotation modulation.

We assume a upper limit period of 10 d for this contaminating effect. The fraction of stars with such short-period companions is only about 2–5% (Raghavan et al. 2010), out of which only 2–5% are grazing binaries (at most), which means 0.04–0.25% of the stars, while we did not yet take into account that the secondaries are substantially smaller in most binaries. Even though we cannot consider the rotation modulation sample a random sub-set of G and later type stars (which is the assumption in above estimate), itseems unlikely that this type of contamination would exceed 1%. The same conclusion is reached from our internal validation against known and newly identified (but not published) eclipsing binaries, based on which we estimate an upper limit of 0.5% grazing binaries contamination.

3.4 Validation of rotation periods

Comparison with the rotational periods of stars in common with the Hartman et al. (2010) and Rebull et al. (2016) dataset is shown in Fig. 6. There are 36 stars in common with Hartman et al. (2010), 4 of which have a period close to the half or double period, and 25 agree to better than 10%. Of the 40 stars in common with Rebull et al. (2016), 2 are close to the half-period, and 36 agree to better than 10%. We therefore conclude that the robustness and accuracy of the Gaia rotational period recovery is comparable to other data published in the literature.

Fig. 6 Left panel: Period comparison with the Hartman et al. (2010) Pleiades data. Black short dashed lines are the half and double period loci, red long dashed lines show the spacecraft rotation alias loci, blue dot-dashed represent the alias loci associated with Gaia spin axis precession, and green short dashed lines show the alias associated with the time delay between the two FoV (basic angle). Right panel: Period comparison with the K2 observations of the Pleiades (Rebull et al. 2016). P is given in day units.

3.5 Period - colour and modulation amplitude distributions

The lowest values of the modulation amplitude have a dependence on the apparent G magnitude that closely mimics the dependence of the photometric sensitivity on G (Fig. 7). In order to simplify the description as best possible, that is, to avoid complications arising from a G-dependent completeness level, we select a subsample in the range 13.5 ≤ G ≤ 15.0 mag, for which the modulation amplitude sensitivity is approximately constant over G.

The selected subsample is further divided into stellar mass bins using M_G as an approximate proxy, which is firstconverted into M_V (Evans et al. 2018, and references therein) and then into mass using the relationships given in Barnes& Kim (2010). In Fig. 8 we show the distributions of the modulation amplitude⁶ over P in ΔM ≈ 0.10M_⊙ bins from 0.65 to 1.25 M_⊙. The data show a clear bimodality with a boundary around maxA_fit ≈ 0.04−0.05. For P < 2 d, there is also a clear gap between the low-amplitude and the high-amplitude population from ≈0.01 to 0.05 mag. The high-amplitude population is more prominent at lower mass and is gradually depleted as we consider higher mass, until it almost disappears in the 1.15–1.25 M_⊙ bin. Conversely, the low-amplitude population is very poorly populated for M ≲ 0.85M_⊙ in the P ≲ 2 d region, but this region is increasingly populated as we consider higher mass.

The low- and high-amplitude populations identified in the amplitude-period diagram also tend to occupy distinctive areas in the period-colour diagram. The period-colour diagrams of these two populations (Fig. 9) show that the low-amplitude population tends to mostly occupy the region where we expect to find a collection of Barnes’ I-sequences (or slow-rotator sequences, Barnes 2003) at different ages, while the high-amplitude population tends to mostly occupy the region where we expect a collection of Barnes’ C-sequences and the gap between the C- and the I-sequence.

We investigated possible instrumental or analysis artefacts that could spuriously produce these distributions, but we found no evidence of this. The bimodal distributions described above can indeed be an indication of two different modes of operation of the hydromagnetic dynamo in stars with the same mass, effective temperature, and rotation period. The Sun has time intervals of prolonged low activity, called grand minima, lasting for several decades. They are separated by intervals of normal activity characterised by the 11-yr cycles and greater variability amplitudes. Grand minima occur for about 17–20% of the time (e.g. Usoskin et al. 2007; Zechmeister & Kürster 2009). In the case of moreactive stars, the possibility that the dynamo is operating in different regimes has also been suggested, for example, to account for the distributions of the rotation periods of late-type stars in open clusters (e.g. Barnes 2003; Brown 2014). Previous observationsdid not clearly detect the bimodality because of the lower photometric precision of ground-based photometry or a rather low percentage of low-amplitude fast-rotators in Kepler data (McQuillan et al. 2014; Reinhold & Gizon 2015).

Contamination by grazing eclipsing binaries (Sect. 3.3) is expected to have negligible effects on the distributions described above. Given the short periods where the bimodality is mostly concentrated, the contaminating eclipsing binaries are expected to have synchronised components. Eclipsing binaries with deep eclipses, larger than about 0.015 mag, for instance, are rejected by our filtering criteria, but grazing binaries or ellipsoidal variables with small amplitude (<0.01 mag) can escape filtering and might be misclassified as rotational variables, especially because the intrinsic variability of their components due to their magnetic activity effectively increases the level of noise in their light curves, hiding the true shape of the modulation. However, since the expected percentage of such objects is well below 1%, they have no significant effects on the distributions shown in Figs. 8 and 9.

Fig. 7 Modulation amplitude (in mag.) dependence on apparent magnitude G. The period-colour diagram and amplitude distribution reported in this paper are built using only the subsample highlighted in orange.

Fig. 8 Amplitude-period diagrams in bins of approximate mass showing a bimodality in modulation amplitude. The colour coding gives an indication of the relative density of the data points (red for lower-density, and purple for higher density). Amplitude units are mag., and P is given in day units.

Fig. 9 Period-colour diagrams for the low-amplitude population (maxAfit ≤ 0.04, upper panel, and for the high-amplitude population maxAfit > 0.04, lower panel). P is given in day units.

4 Conclusions

We have presented the method we used to identify and analyse rotational modulation variables of the BY Dra class in the second Gaia data release. This contains 147 535 BY Dra candidates, together with their photometric time series and the essential parameters related to their flux modulation induced by surface inhomogeneities and rotation. Although this sample represents only a few percent of all BY Dra that can be studied with Gaia, this is already by far the largest BY Dra catalogue available to date. False positives are estimated to be of the order of 5%. The main contaminants are essentially limited to low-mass short-period grazing binaries or binary ellipsoidal variables whose light curves are not far from being sinusoidal. Their expected fraction, however, is expected to be below 1%, making their effects on the distributions over the parameters of interest negligible.

The richness of the sample, the range of periods, and the photometric precision make it possible for the first time to unveil a clear bimodality in the amplitude-period diagrams binned by mass, with a gap between the two populations at P < 2 d. The high-amplitude population is particularly evident at lower stellar mass (≲ 1M_⊙) and is gradually depleted at increasing mass. Conversely, the low-amplitude population at P < 2 d is particularly rich at higher mass (≳1M_⊙) and is increasingly depleted at lower mass. Furthermore, the low-amplitude population is mainly located in regions of the period-colour diagram resembling a collection of Barnes’ I-sequences at different ages (Barnes 2003), while the high-amplitude population is mainly located in regions resembling C-sequences and the gaps between the two. A deeper analysis would require the inclusion of stellar age and G-dependent completeness estimates, but this is beyond the scope of this paper. However, the very existence of this bimodality is further confirmation of the existence of different modes of operation of the dynamo in stars with the same mass, effective temperature, and rotation period.

The possibility of distinguishing such new details in the sample distribution over the parameters of interest, which independently confirm a scenario that is in line with previous studies, is by itself an indication of the validity and quality of the BY Dra Gaia DR2 data. Comparison of the Gaia rotation periods with literature data, although limited to a few tens of cases, confirms that the accuracy and robustness is at the same level as those of other surveys, including, notably, those conducted with Kepler.

Overall, the preliminary data contained in Gaia DR2 illustrate the vast and unique information that the mission is going to provide on stellar rotation and magnetic activity, opening new and unique opportunities in our understanding of the evolution of stellar angular momentum and dynamo action.

Acknowledgements

This work presents results from the European Space Agency (ESA) space mission Gaia. Gaia data are being processed by the Gaia Data Processing and Analysis Consortium (DPAC). Funding for the DPAC is provided by national institutions, in particular the institutions participating in the Gaia MultiLateral Agreement (MLA). This work was supported by the Italian funding agencies Agenzia Spaziale Italiana (ASI) through grants I/037/08/0, I/058/10/0, 2014-025-R.0, and 2014- 025-R.1.2015 to INAF and contracts I/008/10/0 and 2013/030/I.0 to ALTEC S.p.A and Istituto Nazionale di Astrofisica (INAF) (PI M.G. Lattanzi). For Switzerland this work was supported by the Swiss State Secretariat for Education, Research and Innovation through the ESA PRODEX program, the “Mesures d’accompagnement”, the “Activités Nationales Complémentaires”, the Swiss National Science Foundation, and the Early Postdoc. Mobility fellowship. We acknowledge the use of TOPCAT (http://www.starlink.ac.uk/topcat/, Taylor 2005). The Gaia mission website is https://www.cosmos.esa.int/gaia. The Gaia Archive website is http://gea.esac.esa.int/archive/. Thanks to Sydney Barnes, AIP Potsdam, for inspiring discussions.

Appendix A Estimating completeness

Limiting our considerations to pre- and main-sequence single stars or wide binaries, there are several factors to be considered:

• In the main sequence, the amplitude of the flux modulation decreases while the period increases with stellar age; rotational modulation is therefore increasingly difficult to detect as the star evolves.

• Solar-like activity is characterised by magnetic cycles, so that rotational modulation is more easily detectable at epochs close to a maximum of the cycle and may be completely undetectable at epochs close to a minimum.

• Stars seen pole-on or quasi-pole-on are expected to produce negligible rotational modulation.

Bearing these considerations in mind and in view of what we showed in Sect. 3.1, in order to estimate the completeness of the BY Dra sample in Gaia DR2, we start by defining a detection efficiency (A.1)

as a functionof the ecliptic coordinates⁷ (λ, β), the period P, the apparent magnitude G, and colour (G_BP − G_RP), such that (A.2)

where N = N(λ, β, P, G, (G_BP − G_RP)) is the number density of detected BY Dra and the true number density of BY Dra. Since we do not know , we make use of two estimates that may reasonably bracket the real value: at one extreme, we assume that all low-mass stars selected for the analysis in the HRD are BY Dra (say ), at the other, we assume that a reference survey (RS) has detected all BY Dra in the field, explored down to its completeness magnitude limit (e.g. ).

To a certain extent, the Gaia DR2 BY Dra sample overlaps in a region around the Pleiades cluster with the HATNet (Hartman et al. 2010) and K2 (Rebull et al. 2016) surveys. The Hartman et al. (2010) catalogue is particularly useful for a comparison with the Gaia data because of its uniform coverage of an approximately square sky region including the Pleiades cluster and an area close to it. Furthermore, this catalogue contains very many non-members in a smaller area than K2, which makes our estimate easier and more realistic. We therefore take this survey as reference for our completeness estimate. Moreover, since the coverage in (λ, β) is still sparse, Eq. (A.2) is evaluated in the direction of the Pleiades and the value obtained taken as representative for the whole sky for the time being, deferring the full analysis implied by Eq. (A.2) to a later stage when a more uniform BY Dra sky coverage will become available with a more extended Gaia baseline.

A sky map of the Hartman et al. (2010) data with the DR2 rotational modulation variable sample overlapped is shown in Fig. A.1. In order to estimate the ratio, the whole sphere in equatorial coordinates has been divided into HEALPix pixels with resolution index 6, that is to say, a division of the whole sphere into 12 × 6² = 49152 pixels of size of 0.83929366 square degrees. In Fig. A.1, HEALPix with stars belonging to both samples are outlined with a colour code indicating the integrated (in P, G, and (G_BP − G_RP)) ratio in each HEALPix. In this case, extends from 0.02 to 5, with a median of 0.29. The extreme values correspond to the borders of either the Hartman et al. (2010) field or the Gaia strips, and therefore we take the median value as a representative estimate.

Fig. A.1 Upper left panel: sky map of the Hartman et al. (2010) sample (red filled dots) with the Gaia rotational modulation variable sample (blue open dots) overlapped. The HEALPix in which we have stars from both samples are outlined, and the colour code indicates the starnumber ratio in each HEALPix. Upper right and bottom panels: Comparison of the P, G, and (GBP − GRP) distributions in Gaia and the Hartman et al. (2010) sample in the overlapping HEALPix. P is given in day units.

Then we take the P, G, and (G_BP − G_RP) distributions in the overlapping HEALPix of the two samples into account. The Gaia sample contains fainter sources than the Hartman sample. On the other hand, the Gaia sample has a more restrictive lower cut of 0.3 d in P and of 0.6 at the bluer end in (G_BP − G_RP). By restricting the comparison samples to similar P, G, and (G_BP − G_RP) ranges, we obtain the comparison map shown in Fig. A.2. The count ratio, in the remaining HEALPix with stars in common, ranges from 0.02 to 2, with a median of 0.14, in 36 common level 6 HEALPix. From this we roughly estimate that the Gaia completeness upper limit is 14% in the common HEALPix. Furthermore, if we assume that this completeness percentage is representative for the whole sky, considering that the Gaia sample covers 18593 out of 49152 level-6 HEALPix, we estimate an upper limit of 5% completeness for the whole sky down to G ≈ 16.5.

Fig. A.2 Comparison of the P, G, and BP − RP distributions between the Gaia DR2 BY Dra and the Hartman et al. (2010) samples in the overlapping HEALPix and in common P, G, and BP − RP ranges. P is given in day units.

At the other extreme, for what matters here, we assume that all late-type dwarfs observed by Gaia have surface inhomogeneities and therefore should display rotational modulation, but Gaia does not reveal it. We then compare the Hartman et al. (2010) sample with all late-type dwarfs observed by Gaia in the same FoV (same level 6 HEALPix). In this way, we find a median count ratio of 0.14, from which we estimate a lower limit of 14% for the completeness of the Hartman et al. (2010) sample and therefore a lower limit of 0.7% for the Gaia DR2 BY Dra sample over the whole sky.

The expected main sources of uncertainty in this simplified completeness estimate are those related to the dependence of the detection efficiency (Eq. (A.1)) on the ecliptic coordinate and on P. Figures 4 and 5 give an idea of the problem: the distribution of detected P depends on the ecliptic coordinates, and the preliminary data contained in Gaia DR2 are still insufficient for a meaningful empirical determination of the Q dependencies. Regarding the dependence on P, taking an existing survey as reference would require, in principle, that such a survey has at least a known P detection efficiency. This is not our case, but taking a survey that explored just a sky region containing the Pleiades cluster (≈120 Myr old) is expected to have no great consequences on our overall percentage completeness estimate. Considering both cluster members and non-members, as available in the Hartman et al. (2010) catalogue, provides a sample sufficiently rich to make our rough estimate meaningful even without a detailed knowledge of the detection efficiency of the reference catalogue P. Further work is required in this direction before the final Gaia data release.

Appendix B Catalogue retrieval

Here we present some examples of ADQL queries to be used in the web interface to the Gaia DR2 archive⁸ to retrieve the BY Dra data.

Get the first 1000 stars in the rotational modulation table:

SELECT TOP 1000 ^★ FROM

gaiadr2.vari_rotation_modulation

Get the position, parallax, and link to the light curves for the first 1000 sources in the rotational modulation table:

SELECT TOP 1000 gs.source_id,gs.ra, gs.dec,

gs.parallax,

epoch_photometry_url FROM gaiadr2.gaia_source

as gs INNER JOIN

gaiadr2.vari_rotation_modulation as rotmod ON

gs.source_id=rotmod.source_id

Get the time-series statistics for the first 1000 sources in the rotational modulation table:

SELECT top 1000 tsstat.^★ FROM

gaiadr2.vari_time_series_statistics as

tsstat INNER JOIN

gaiadr2.vari_rotation_modulation as rotmod ON

tsstat.source_id=rotmod.source_id

Get the final rotation period and amplitude of the first 100 sources in the rotational modulation table:

SELECT top 100 source_id,

best_rotation_period, max_activity_index

FROM gaiadr2.vari_rotation_modulation

Select parallax, magnitudes, period, and modulation amplitude for the first 100 sources in the rotational modulation table:

SELECT TOP 100 gs.source_id, gs.parallax,

gs.phot_g_mean_mag,

gs.phot_bp_mean_mag, gs.phot_rp_mean_mag,

rm.best_rotation_period,

rm.max_activity_index FROM gaiadr2.gaia_source

AS gs INNER JOIN

gaiadr2.vari_rotation_modulation AS rm ON

gs.source_id=rm.source_id

References

All Figures

	Fig. 1 Region of the HRD in which rotational modulation was found. The region is highlighted with colours (rainbow). The colour code indicates the relative density of data points (red for lower density, and purple for higher density). The HRD shown is built using the AGIS02.2 solution for the parallax (left panel, see Lindegren et al. 2018, for details) uncorrected for interstellar extinction.
In the text

	Fig. 2 Simulated light curve for a solar-like star with an effective temperature Teff = 5000 K, a rotation period P = 5 d, and an average Tspot = 4200 K.
In the text

	Fig. 3 Period-folded light curves for the Gaia DR2 source id 894812322514469120. Three similar significantperiods were found in two segments and in the whole time series with similar Afit and Aper.
In the text

	Fig. 4 Sky distribution of Gaia DR2 BY Dra candidates. P is given in day units.
In the text

Fig. 5 Rotation period distribution of Gaia DR2 data for BY Dra stars (P is given in day units). A 3 d binning P histogram is shown in the left panel, together with the Gaussian kernel density estimate with a 3 d width and 3-sigma truncation. The P – ecliptic latitude scatter plot (right panel) shows the higher number of period measurements at ± 45° ecliptic latitudes, particularly regarding longer periods. The colour code (rainbow) is an indication of the relative density of data points (red for lower density, and blue for higher density). Both panels give an indication of the moderate degree of aliasing, the most evident being at ≈23.5, 25.3, and 31.6 d.

In the text

Fig. 6 Left panel: Period comparison with the Hartman et al. (2010) Pleiades data. Black short dashed lines are the half and double period loci, red long dashed lines show the spacecraft rotation alias loci, blue dot-dashed represent the alias loci associated with Gaia spin axis precession, and green short dashed lines show the alias associated with the time delay between the two FoV (basic angle). Right panel: Period comparison with the K2 observations of the Pleiades (Rebull et al. 2016). P is given in day units.

In the text

	Fig. 7 Modulation amplitude (in mag.) dependence on apparent magnitude G. The period-colour diagram and amplitude distribution reported in this paper are built using only the subsample highlighted in orange.
In the text

	Fig. 8 Amplitude-period diagrams in bins of approximate mass showing a bimodality in modulation amplitude. The colour coding gives an indication of the relative density of the data points (red for lower-density, and purple for higher density). Amplitude units are mag., and P is given in day units.
In the text

	Fig. 9 Period-colour diagrams for the low-amplitude population (maxAfit ≤ 0.04, upper panel, and for the high-amplitude population maxAfit > 0.04, lower panel). P is given in day units.
In the text

Fig. A.1 Upper left panel: sky map of the Hartman et al. (2010) sample (red filled dots) with the Gaia rotational modulation variable sample (blue open dots) overlapped. The HEALPix in which we have stars from both samples are outlined, and the colour code indicates the starnumber ratio in each HEALPix. Upper right and bottom panels: Comparison of the P, G, and (GBP − GRP) distributions in Gaia and the Hartman et al. (2010) sample in the overlapping HEALPix. P is given in day units.

In the text

	Fig. A.2 Comparison of the P, G, and BP − RP distributions between the Gaia DR2 BY Dra and the Hartman et al. (2010) samples in the overlapping HEALPix and in common P, G, and BP − RP ranges. P is given in day units.
In the text