Abstract

We present an analytic model of the light-curve variation for stars with non-evolving starspots on a differentially rotating surface. The Fourier coefficients of the harmonics of the rotation period are expressed in terms of the latitude of the spot, ℓs, and the observer’s line-of-sight direction, ℓo, including the limb-darkening effect. We generate different realizations of multi-spots according to the model, and perform mock observations of the resulting light-curve modulations. We discuss to what extent one can recover the properties of the spots and the parameters for the differential rotation law from the periodogram analysis. Although our analytical model neglects the evolution of spots on the stellar surface (dynamical motion, creation, and annihilation), it provides a basic framework to interpret the photometric variation of stars, in particular from the existing Kepler data and the future space-born mission. It is also applicable to photometric modulations induced by rotation of various astronomical objects.

1 Introduction

The last two decades have seen the birth and growth of space-borne photometry due to missions like MOST (Walker et al. 2003), CoRoT (Baglin et al. 2006), Kepler (Borucki et al. 2010), and TESS (Ricker et al. 2014). These instruments provided for the first time long, continuous high-quality photometric data, enabling the detection of thousands of transiting planets, but also opened a new window on the dynamics and evolution of stars. For example, one of the most remarkable achievements of the Kepler space instrument is the discovery that older low-mass stars rotate too fast compared to theoretical expectations (van Saders et al. 2016). This could only be established by the combined analysis of the stellar photometric variability due to spots and to the stellar pulsations (see, e.g., Kjeldsen & Bedding 1995; Christensen-Dalsgaard et al. 1996; Chaplin et al. 2010). In general, the study of the rotation–age relation (Skumanich 1972; Kawaler 1988; MacGregor & Brenner 1991) is an important tool to evaluate the age of stars. The rotation also plays an important role on the solar and stellar dynamo (Ossendrijver 2003; Varela et al. 2016), itself believed to be important for sustaining a latitudinal differential rotation.

Observationally, the stellar rotation period can be estimated from a few independent methods. First, one can combine the equatorial rotational velocity from Doppler broadening and the stellar radius. The spectroscopically derived rotation period, however, depends on the assumed model for the turbulence, and also requires values of the stellar radius and inclination that are not well-determined in general (Kamiaka et al. 2018). Secondly, the asteroseismic analysis of the stellar pulsation can estimate the rotation period and the stellar inclination simultaneously. The asteroseismology, however, also requires various model assumptions in the analysis, and is applicable only to a relatively small fraction of stars that exhibit measurable oscillations (e.g., Appourchaux et al. 2008; Huber et al. 2013; Benomar et al. 2014; Lund et al. 2017; Kamiaka et al. 2018, 2019). Finally, the photometric variation of light curves is by far the most widely used method to estimate the stellar rotation period, and has been intensively applied to the Kepler data (e.g., McQuillan et al. 2014; Mazeh et al. 2015; Angus et al. 2018).

The photometric variation is induced by starspots corotating with the star. If those spots do not dynamically evolve on the stellar surface, it is relatively easy to estimate the stellar rotation period. In reality, however, the spots have individual lifetimes and even move on the stellar surface, and stars are not necessarily rigid rotators (e.g., Donati & Collier Cameron 1997; Barnes et al. 2005; Donati et al. 2010; Roettenbacher et al. 2013; Walkowicz et al. 2013; Brun et al. 2017; Benomar et al. 2018; Basri & Shah 2020). The formation and dissipation of spots on the differentially rotating stars, therefore, complicate the interpretation of the photometrically estimated rotation period Pphoto. Furthermore, we cannot exclude the possibility that that starspots and stellar pulsations may have similar time scales in some stellar types, even if not so likely.

Properties of spots have been extensively studied in past literature for the Sun (e.g., Maunder 1904; Zharkov et al. 2005; Mandal et al. 2021), and also for other stars (Morris 2020). Roettenbacher et al. (2013), for instance, achieved a wonderful light-curve inversion to predict the starspot evolution on Kepler target KIC 5110407. Nevertheless, it is not easy to predict the nature of spots accurately in general. On the other hand, the photometric rotation periods combined with the spectroscopic Doppler broadening have been extensively used to infer the inclination angle of stars hosting planets (Sanchis-Ojeda et al. 2011; Sanchis-Ojeda & Winn 2011; Hirano et al. 2012; Louden et al. 2021; Albrecht et al. 2021), which have profound implications for the spin-orbit architecture of exoplanetary systems (Queloz et al. 2000; Ohta et al. 2005; Sasaki & Suto 2021). Therefore, it is still useful to have parametrized templates for the photometric variation of stellar light curves due to non-evolving starspots. This is the purpose of this paper. We present an analytic model of the photometric light curves induced by starspots on a differentially rotating stellar surface assuming that they do not evolve during the finite observing duration. We compute mock light curves based on the multi-spot model, and address how to interpret the measured distribution of peaks in the Lomb–Scargle periodogram in terms of the stellar differential rotation law.

The rest of the paper is organized as follows. We derive the photometric variation pattern due to a single infinitesimal spot in section 2. The resulting light-curve modulation including the limb-darkening effect is expressed in the Fourier series expansion. In section 3, we apply the analytic model for multispots on a differentially rotating star. Then we generate simulated light curves adopting the statistical distribution model of the Sun spots, perform the Lomb–Scargle analysis, and examine the information content of the resulting power spectra. The final section is devoted to the summary and conclusion of the paper. The Fourier expansion coefficients in our analytic model are given in the appendices.

2 Photometric variation due to a single starspot

As illustrated in figure 1, we consider a spherical star with radius R, and parametrize a position vector on the stellar surface in terms of its latitude ℓ and longitude ϕ:
(1)
where the z-axis is chosen to be the direction of the stellar rotation. In what follows, we assume that the surface angular velocity ω(ℓ) at ℓ is given by the following parametrized model for the latitudinal differential rotation:
(2)
For the Sun, ω0⊙ ≈ 2.972 × 10−6 rad s−1, α2⊙ ≈ 0.163, and α4⊙ ≈ 0.121 (Snodgrass & Ulrich 1990). Thus, the angular velocity at ℓ = 30° is about |$5\%$| smaller than its equatorial value.
Schematic illustration of the observer and a starspot. The stellar rotation axis is chosen to be the z-axis. The location of the spot is specified by the latitude ℓs and the azimuth angle ϕs, and the direction of the observer’s line-of-sight is defined by the latitude ℓo on the x–z plane (ϕ = 0).
Fig. 1.

Schematic illustration of the observer and a starspot. The stellar rotation axis is chosen to be the z-axis. The location of the spot is specified by the latitude ℓs and the azimuth angle ϕs, and the direction of the observer’s line-of-sight is defined by the latitude ℓo on the xz plane (ϕ = 0).

Without loss of generality, we consider a distant observer located at ϕo = 0 and ℓo. Thus the unit vector toward the observer is
(3)
The longitude of the starspot located at the latitude ℓs at epoch t becomes
(4)
due to the stellar surface rotation, where ϕs0 is the longitude at which the spot is located on the stellar surface initially (t = 0), and ωs ≡ ω(ℓs) is the angular velocity of the spot at ℓs defined as equation (2).

2.1 A single infinitesimal starspot without limb darkening

A normalized lightcurve of the stellar surface is
(5)
where the integration is over the stellar surface, A indicates the surface intensity distribution, and K is the weighting kernel of the surface visible for the observer (e.g., Fujii et al. 2010, 2011; Farr et al. 2018; Haggard & Cowan 2018; Nakagawa et al. 2020).
In the case of a single infinitesimal starspot at (ℓs, ϕs) on a homogeneous sphere, we set
(6)
The dimensionless parameter bspot represents the amplitude of the photometric modulation (Dorren 1987; Haggard & Cowan 2018).

Sunspots consist of the central darker part (umbra), and the surrounding lighter part (penumbra). In addition, there is a type of brighter spots (faculae). Note that it is not necessary to specify the temperature and area of spots separately, since the amplitude of the photometric variations of stars depends on their flux (i.e., bspot) alone. As described in subsection 3.2, we consider the spot distribution directly derived from the observed properties using the photometric variation data of the Sun.

Our analytic formulation, however, is general and applicable to various types of spots including umbra, penumbra, and faculae (bspot < 0), if we employ their distribution function properly.

Without loss of generality, we can define the initial phase of the single starspot to be ϕs(t = 0) = ϕs0 = 0. For an isotropically emitting stellar surface, the weighting kernel K is equivalent to the visibility computed from the direction cosine |$\mu \equiv {\boldsymbol e}_{\rm o}\cdot {\boldsymbol e}_{\star }$| between the stellar surface |$\boldsymbol {r}_\star =R_{\star }\boldsymbol {e}_\star$| and the observer |${\boldsymbol e}_{\rm o}$|⁠. The spot is visible (invisible) to the observer if μ > 0 (μ < 0). Thus the weighting kernel is simply computed from equations (1) and (3) as
(7)
Substituting equations (6) and (7) into equation (5), one obtains a normalized light-curve modulation due to a single starspot on an otherwise homogeneous spherical surface:
(8)
For ℓo = ℓs = ϕs(t) = 0, equation (8) reduces to Ls(t) = −bspot/π, and the denominator π indeed corresponds to the visible projected area of the stellar surface |$\pi R_{\star }^2$|⁠. Thus, we note that bspot represents the effective area of the spot Aspot in units of |$R_{\star }^2$|⁠, instead of |$\pi R_{\star }^2$|⁠. Since bspot in equation (6) is defined with respect to the flux, Aspot is equivalent to the geometric area of the spot S only when it is completely black. In general, Aspot should be interpreted to represent a flux-weighted area of the spot.

In what follows, we adopt a parametrized model of the distribution of Aspot that is directly estimated from the observed photometric variations of the Sun (see subsection 3.2). Then we will compute the dimensionless parameter |$b_{\rm spot}\equiv A_{\rm spot}/R_{\star }^2$|⁠. If the blackbody approximation for the stellar surface and the spot is valid, the effective and geometric areas of the spot are related as |$A_{\rm spot}\approx (1-T_{\rm spot}^4/T_\star ^4)S$|⁠, with T and Tspot being the temperatures of the star and the spot.

Equation (8) indicates that the starspot is visible at t if
(9)
For convenience, let us introduce the parameter
(10)
A starspot with Γ > 1 is always visible to the observer, and equation (8) reduces to
(11)
If Γ < −1, on the other hand, the starspot is totally invisible and Ls(t) = 0.
A starspot with |Γ| ≤ 1 becomes visible periodically as the stellar surface rotation. In this case, equation (8) is expanded analytically in the Fourier series. The result is
(12)
where the parameter θc is defined through Γ = tan ℓotan ℓs ≡ −cos θc(0 < θc < π); see appendix 1 for the derivation of equation (12).

The visibility of a single starspot is determined by the parameter Γ or, equivalently, θc. We plot the contours of Γ and θc on the ℓo–ℓs plane in the left- and right-hand panels of figure 2, respectively. For a roughly edge-on-view observer (|ℓo| ≪ 1), spots located near the equatorial plane (|ℓs| ≪ 1) correspond to |Γ| ≈ |ℓos| ≪ 1, and θc ≈ π/2 + ℓos.

Contours of Γ(≡ tan ℓotan ℓs) and θc(≡ −cos −1Γ for |Γ| < 1) on ℓo–ℓs plane. Red and blue solid lines in both panels indicate contours for 0 ≤ Γ ≤ 1 (90° ≤ θc ≤ 180°) and −1 ≤ Γ < 0 (0° ≤ θc < 90°), for which the corresponding starspot becomes visible periodically to the observer as the star rotates. The orange and black dotted lines correspond to those spots that are always visible and invisible to the observer, respectively.
Fig. 2.

Contours of Γ(≡ tan ℓotan ℓs) and θc(≡ −cos −1Γ for |Γ| < 1) on ℓo–ℓs plane. Red and blue solid lines in both panels indicate contours for 0 ≤ Γ ≤ 1 (90° ≤ θc ≤ 180°) and −1 ≤ Γ < 0 (0° ≤ θc < 90°), for which the corresponding starspot becomes visible periodically to the observer as the star rotates. The orange and black dotted lines correspond to those spots that are always visible and invisible to the observer, respectively.

2.2 A single infinitesimal starspot with limb darkening

The stellar limb darkening produces an additional modulation to the photometric variation due to the starspot. Adopting the quadratic limb-darkening law, the normalized stellar surface intensity at |$\boldsymbol {r}_\star$| is characterized by the two limb-darkening parameters u1 and u2 as
(13)
where |$\mu = {\boldsymbol e}_{\rm o}\cdot {\boldsymbol e}_{\star }$| is the direction cosine that we defined before. We adopt the values of u1 and u2 from the Sun (Cox 2000): u1 = 0.47 and u2 = 0.23 at 550 nm (they become 0.42 and 0.23, respectively, at 600 nm).
Including the limb-darkening effect, equation (5) is generalized to be
(14)
Since the denominator of equation (14) is
(15)
equation (8) is now written as
(16)
where μs = cos ℓocos ℓs(cos ωst + Γ).
Similarly to the previous subsection, equation (16) for |Γ| ≤ 1 is expanded analytically in the Fourier series. The derivation is explicitly given in appendix 1, and the normalized light-curve modulation including the limb-darkening effect is summarized in the following expression:
(17)
where
(18)
and the coefficients an, bn, and cn are explicitly given in appendix 1.

For spots with Γ > 1, μs is always positive, and the corresponding lightcurve is written in the same form as equation (16) by replacing An by |$\tilde{A}_n$|⁠, which are given in appendix 2.

Figures 3, 4, and 5 show the trajectories of a single spot on a rotation stellar surface and the corresponding normalized light curves Ls(t)/bspot against t/Pspin(ℓs) for an observer located at ℓo = 0°, 30°, and 60°, respectively. Black, red, blue and orange curves indicate the results for the spot at the latitude of ℓs = 0°, 15°, 45°, and 75°. Solid and dashed lines in the right-hand panels indicate the light curves with and without limb darkening (LD).

Photometric modulation due to a single spot at ℓs viewed from the observer at ℓo = 0°. Left: Trajectories of four spots at ℓs = 0° (black), 15° (red), 45° (blue), and 75° (orange) on the stellar surface. Right: Modulation curves Ls/bspot over the one rotation period of each spot Pspin(ℓs) in the left-hand panel. Solid and dashed lines correspond to those with and without limb darkening (LD) for the differential rotation parameters of α2 = α2⊙ and α4 = α4⊙; see equation (2).
Fig. 3.

Photometric modulation due to a single spot at ℓs viewed from the observer at ℓo = 0°. Left: Trajectories of four spots at ℓs = 0° (black), 15° (red), 45° (blue), and 75° (orange) on the stellar surface. Right: Modulation curves Ls/bspot over the one rotation period of each spot Pspin(ℓs) in the left-hand panel. Solid and dashed lines correspond to those with and without limb darkening (LD) for the differential rotation parameters of α2 = α2⊙ and α4 = α4⊙; see equation (2).

Same as figure 3 but viewed from the observer at ℓo = 30°.
Fig. 4.

Same as figure 3 but viewed from the observer at ℓo = 30°.

Same as figure 3 but viewed from the observer at ℓo = 60°.
Fig. 5.

Same as figure 3 but viewed from the observer at ℓo = 60°.

The right-hand panel of figure 3 shows that the modulation amplitude |Ls(t)/bspot| without the limb-darkening effect becomes 1/π for ℓo = ℓs = 0° (black-dashed curve) at ϕs(t) = 0. Limb darkening decreases the effective visible area of the entire surface by a factor of (1 − u1/3 − u2/6)−1, while that of the starspot by a factor of Is). Depending on the location of the spot, ℓo, ℓs, and ϕs(t), the resulting |Ls(t)/bspot| with limb darkening becomes either smaller or larger than that without limb darkening; see figures 3, 4, and 5.

Figure 6 plots the ratios of Fourier coefficients of the single spot modulation, An/A1 (n = 2, 3, 4). If limb darkening is neglected, they reduce to an/a1 that are a function of Γ ≡ tan ℓotan ℓs (or θc) alone, which are plotted in dotted lines. When the limb-darkening effect is taken into account, An/A1 depends on both ℓo and ℓs. As figure 6 implies, however, the difference among the three curves for ℓo = 10°, 30°, and 60° is small. Thus, An/A1 is still largely determined by the value of Γ (or θc) even with limb darkening.

Ratios of Fourier coefficients of the photometric modulation due to a single spot, An/A1, plotted against θc (left) and Γ (right). Red, blue, and orange lines indicate A2/A1, A3/A1, and A4/A1, respectively. We assume the limb-darkening parameters of α2 = α2⊙ and α4 = α4⊙, and plot those ratios in solid (ℓo = 10°), dashed (ℓo = 30°), and dot–dashed (ℓo = 60°) lines. For reference, the results without limb darkening (w/o LD) are plotted in dotted lines.
Fig. 6.

Ratios of Fourier coefficients of the photometric modulation due to a single spot, An/A1, plotted against θc (left) and Γ (right). Red, blue, and orange lines indicate A2/A1, A3/A1, and A4/A1, respectively. We assume the limb-darkening parameters of α2 = α2⊙ and α4 = α4⊙, and plot those ratios in solid (ℓo = 10°), dashed (ℓo = 30°), and dot–dashed (ℓo = 60°) lines. For reference, the results without limb darkening (w/o LD) are plotted in dotted lines.

This result suggests that An/A1 may be used to examine if the periodic signals detected from the observed photometric lightcurve are due to starspots, instead of other sources. It may be even possible to put a constraint on Γ from An/A1 in principle. Since ℓo is equivalent to the stellar inclination for the observer that can be independently measured from either spectroscopy or asteroseismology (e.g., Kamiaka et al. 2018, 2019; Sasaki & Suto 2021), the constraint on Γ is translated to that on the spot latitude ℓs. In reality, it is feasible to derive a robust constraint on Γ only for a single-spot case. The statistical distribution of An/A1 for multi-spots is more useful to constrain the differential rotation, as discussed below.

3 Multiple starspots: model predictions and mock data analysis

3.1 Superposition of multiple starspots

If more than one starspot is involved, we have to take into account their relative phases, namely ϕs0 in equation (4), as well. In that case, equation (17) can be generalized to
(19)
 
(19)
Thus, the light curve due to multispots becomes the superposition of the following form:
(20)
In the above equation, Ns and |$\tilde{N}_{\rm s}$| denote the number of spots with |Γ| < 1 and Γ > 1, respectively, bspot,i is the amplitude of the photometric variation, ϕs0,i is the initial phase, ωs,i = ω(ℓs,i) is the angular frequency, and An,i and |$\tilde{A}_{n,i}$| are the Fourier components of the ith starspot.

3.2 Mock light curves and the Lomb–Scargle power spectra

In order to examine to what extent one can extract the characteristic signature of starspots from photometric stellar light curves, we create mock light curves in the time domain, and compute the Lomb–Scargle power spectra. Our fiducial set of parameters is listed in table 1.

Table 1.

Fiducial parameters for mock light curves.

SymbolRangeNote
bspotk = 0.54, λ = 2.88 μHemThe Weibull distribution
Ntot30Total number of generated starspots
ϕs[0, 2π]Uniform
s|ℓs| < ℓs,max = 30°P(ℓs) ∝ |sin ℓs|
oStellar inclination relative to the observer’s line-of-sight
α2α2⊙ = 0.163Differential rotation coefficient
α4α4⊙ = 0.121Differential rotation coefficient
u10.47Linear limb-darkening parameter
u20.23Quadratic limb-darkening parameter
2π/ω010 dEquatorial rotation period
Tsamp30 minCadence of the observation
Tobs90 dDuration of the observation
SymbolRangeNote
bspotk = 0.54, λ = 2.88 μHemThe Weibull distribution
Ntot30Total number of generated starspots
ϕs[0, 2π]Uniform
s|ℓs| < ℓs,max = 30°P(ℓs) ∝ |sin ℓs|
oStellar inclination relative to the observer’s line-of-sight
α2α2⊙ = 0.163Differential rotation coefficient
α4α4⊙ = 0.121Differential rotation coefficient
u10.47Linear limb-darkening parameter
u20.23Quadratic limb-darkening parameter
2π/ω010 dEquatorial rotation period
Tsamp30 minCadence of the observation
Tobs90 dDuration of the observation
Table 1.

Fiducial parameters for mock light curves.

SymbolRangeNote
bspotk = 0.54, λ = 2.88 μHemThe Weibull distribution
Ntot30Total number of generated starspots
ϕs[0, 2π]Uniform
s|ℓs| < ℓs,max = 30°P(ℓs) ∝ |sin ℓs|
oStellar inclination relative to the observer’s line-of-sight
α2α2⊙ = 0.163Differential rotation coefficient
α4α4⊙ = 0.121Differential rotation coefficient
u10.47Linear limb-darkening parameter
u20.23Quadratic limb-darkening parameter
2π/ω010 dEquatorial rotation period
Tsamp30 minCadence of the observation
Tobs90 dDuration of the observation
SymbolRangeNote
bspotk = 0.54, λ = 2.88 μHemThe Weibull distribution
Ntot30Total number of generated starspots
ϕs[0, 2π]Uniform
s|ℓs| < ℓs,max = 30°P(ℓs) ∝ |sin ℓs|
oStellar inclination relative to the observer’s line-of-sight
α2α2⊙ = 0.163Differential rotation coefficient
α4α4⊙ = 0.121Differential rotation coefficient
u10.47Linear limb-darkening parameter
u20.23Quadratic limb-darkening parameter
2π/ω010 dEquatorial rotation period
Tsamp30 minCadence of the observation
Tobs90 dDuration of the observation
The key parameter characterizing the spot in our model is bspot. As described in subsection 2.1, bspot is defined as |$A_{\rm spot}/R_{\star }^2$| in our model. Muñoz-Jaramillo et al. (2015) found that the flux-weighted effective area, Aspot, for the solar spot empirically obeys the Weibull distribution,
(21)
from the observed photometric variation over years. The Weibull distribution is written in terms of Aspot/λ, and the amplitude of the resulting spot modulation is simply scaled to the adopted value of λ.
The expectation value of Aspot from equation (21) is
(22)
where Γ(x) denotes the Gamma function and Γ[1 + (1/k)] ≈ 1.75 for the solar value of k = 0.54. Also, the corresponding cumulative number distribution of Aspot exceeding the threshold value Ath is
(23)
For instance, the top 10th percentile of spots have Aspot > 1.47λ.

The best-fitting values of the two parameters, k and λ, vary for different definitions of spots and different datasets (Muñoz-Jaramillo et al. 2015). For definiteness, we adopt the “Sunspot Umbral Area” from the Helio-seismic and Magnetic Imager on the Solar Dynamics Observatory (see their table 1), and adopt k = 0.54 and |$\lambda = 2.88 \, \mu \mbox{Hem}=2.88\times 10^{-6} (2\pi R_{\star }^2)$|⁠. It is likely that different stars may have different values of k and λ. Since our model is fully analytical, however, it is readily applicable for other choices. Thus, we fix their values below, and generate mock data for multi-spots.

Equation (22) suggests that a characteristic amplitude of the dimensionless parameter bspot of our spots is
(24)
where the factor of 2π comes from the fact that λ is given relative to the area of the hemisphere (Hem), |$2\pi R_{\star }^2$|⁠. Thus we can safely neglect the finite size effect of an individual spot, which is consistent with the assumptions of our analytic model, for our adopted value of λ = 2.88 μHem.
In order to understand the meaning of equation (24), let us define the effective radius of the spot rspot through
(25)
Substituting equation (24) into equation (25), one obtains
(26)
or equivalently
(27)
Equations (26) and (27) correspond to the angular and real size corresponding to rspot in terms of bspot.
We generate Ntot spots with bspot following the Weibull distribution, equation (21). We adopted Ntot = 30 for definiteness so as to reproduce the solar spots roughly. The corresponding fraction of spots over the entire stellar surface may be computed from equation (22):
(28)
The value of Ntot is sensitive to the threshold value Ath in identifying a single spot even for the Sun, and moreover is not clear for other stars. Our analytic formulation can be applied to a different choice of Ntot in a straightforward manner.

The latitudes of spots, ℓs, are drawn from the isotropic distribution function (∝|sin ℓs|) but over the restricted range of −ℓs,max < ℓs < ℓs,max. We choose ℓs,max = 30° as our fiducial value, but consider 75° as well to examine its impact. The initial phases are selected randomly for 0 < ϕs < 2π.

For a given value of the observer’s latitude ℓo, we classify each spot according to |Γ| < 1 and Γ > 1, and compute the number of such spots Ns and |$\tilde{N}_{\rm s} (=N_{\rm tot}-N_{\rm s})$|⁠, respectively. Then the light-curve modulation due to those spots is computed from equation (20).

We generate the mock light curves with cadence Tsamp over the duration of Tobs. We set the fiducial values as Tsamp = 30 min and Tobs = 90 d, following the long cadence observation for one single quarter of the Kepler dataset.

Finally, we add the Gaussian noise to the light curves:
(29)
In what follows, we consider two cases, σn = 0 (noiseless) and σn = 35 ppm as a typical value for the Kepler data (cf. Walkowicz et al. 2013; Basri & Shah 2020), for simplicity. Equation (21) implies that the flux modulation induced by a single spot is typically much smaller than the noise:
(30)
Thus, in the case of σn = 35 ppm, the clear periodic signal is visible only for a relatively big spot (Aspot > 5λ, roughly corresponds to the top 10th percentile) or a clustered group of nearby spots.

Figure 7 shows the mock data for two different sets of realizations of starspots (Ntot = 30 ℓs,max = 30°) for an observer at ℓo = 0°. The fractional area covered by spots varies from 0.3 to 1.3 times the expectation value of equation (28).

Light-curve modulations due to starspots (Ntot = 30, ℓs,max = 30°) for an observer at ℓo = 0°. Left-hand, center, and right-hand panels show the spot distribution on the stellar surface, normalized light curve, and the corresponding Lomb–Scargle power spectrum. Results for σn = 0 and 35 ppm are plotted in blue and red, respectively. Upper and lower panels are different realizations of the statistically same model. The vertical dotted lines indicate the range of the differentially rotation periods; Pspin(ℓs = 0°)/2, Pspin(ℓs = ℓs,max)/2, Pspin(ℓs = 0°), and Pspin(ℓs = ℓs,max).
Fig. 7.

Light-curve modulations due to starspots (Ntot = 30, ℓs,max = 30°) for an observer at ℓo = 0°. Left-hand, center, and right-hand panels show the spot distribution on the stellar surface, normalized light curve, and the corresponding Lomb–Scargle power spectrum. Results for σn = 0 and 35 ppm are plotted in blue and red, respectively. Upper and lower panels are different realizations of the statistically same model. The vertical dotted lines indicate the range of the differentially rotation periods; Pspin(ℓs = 0°)/2, Pspin(ℓs = ℓs,max)/2, Pspin(ℓs = 0°), and Pspin(ℓs = ℓs,max).

We search for periodic signals of an angular frequency ω,
(31)
embedded in the mock light curves (center panels) using the Lomb–Scargle (LS) method (Lomb 1976; Scargle 1982). The conventional LS adopts nterms = 1, and we compute the normalized power spectrum
(32)
where χ2(ω) is the residuals of the fit, with |$\chi _{\rm ref}^2$| being the reference value for a constant model.

The left-hand panels of figure 7 show the spot distribution at t = 0 for two different realizations. The area of each circle is plotted in proportion to bspot, and approximately represents the true ratio of the spot area and the entire stellar surface (but neglecting the distortion due to the projection on to the plane). Since those spots span a range of latitudes, the resulting light curves (center panels) are not exactly periodic in the time domain due to the latitudinal surface differential rotation. The right-hand panels of figure 7 show the corresponding LS power. The highest peaks around 10 d are located over a range of rotation periods spanning Pspin(ℓs = 0°) and Pspin(ℓs = ℓs,max) due to the differential rotation. The secondary peaks around 5 d are the second harmonics. The ratios of those amplitudes carry important information on ℓs and ℓo, and will be discussed later (subsection 3.3) using the LS analysis with nterms = 2.

Our mock data completely neglect the dynamics of spots (their creation and dissipation, and motion on the stellar surface) over the duration of the observation Tobs = 90 d, which corresponds to the duration of a single quarter of the Kepler dataset. In order to empirically evaluate the effects of the spot dynamic, we create 10 totally independent realizations drawn from the statistically same spot distribution, and compute each LS power and the average over the 10 realizations. The latter may be interpreted as the average LS power of the entire Kepler observing period, which is made of up to 10 quarters.

The results are shown in figure 8. Possible signatures of differential rotation may be found in the variance among the LS power spectra for different quarters. The width of a peak with a detected period is determined by the entire duration of the observation Tobs, instead of the cadence Tsamp in the present examples. For instance, one can resolve the periods for different spots if they are static over Tobs = 900 d, but cannot do so for Tobs = 90 d. While the non-evolving spots over Tobs = 900 d may not be so realistic in general, a small fraction of stars may have such spots. Therefore our study suggests that it is worthwhile to attempt searching for such signatures in the Kepler archive data.

LS power spectra for different realizations corresponding to figure 7. Upper and lower panels are for σn = 0 and 35 ppm, and left- and right-hand panels are without and with differential rotation. Thin curves (10 in total) indicate the results for different realizations of the spots for Tsamp = 30 min and Tobs = 90 d, with thick red lines being their average. Thick blue lines show the LS power for one realization but observed for Tobs = 900 d.
Fig. 8.

LS power spectra for different realizations corresponding to figure 7. Upper and lower panels are for σn = 0 and 35 ppm, and left- and right-hand panels are without and with differential rotation. Thin curves (10 in total) indicate the results for different realizations of the spots for Tsamp = 30 min and Tobs = 90 d, with thick red lines being their average. Thick blue lines show the LS power for one realization but observed for Tobs = 900 d.

Figures 9 and 10 show the same plots as figures 7 and 8, but for a wider distribution of spots (ℓs,max = 75°) observed from an observer located far outside the stellar equatorial plane (ℓo = 45°). As expected, the effect of differential rotation is more visible than that for ℓs,max = 30° and ℓo = 0°.

Same as figure 7, but for ℓs,max = 75° and ℓo = 45°.
Fig. 9.

Same as figure 7, but for ℓs,max = 75° and ℓo = 45°.

LS power spectra for different realizations corresponding to figure 9.
Fig. 10.

LS power spectra for different realizations corresponding to figure 9.

The visible periodicity of the light-curve modulation in the center panels of figures 7 and 9 seems to be generated by a relatively small number of large spots. To clarify this point, we repeated the analysis by dividing the 30 spots in the two realizations of figures 7 and 9 into two separate groups; the top 10 spots and the remaining 20 spots. The resulting plots are shown in figures 11 and 12. While those small spots still show periodic signals in the noiseless light curve, they are substantially buried in the case of our adopted noise of σn = 35 ppm. In other words, the peaks in the LS power spectra are dominated by a small fraction of spots, and should mostly represent their properties (size, latitude, and rotation velocity), as long as the Weibull distribution is a good approximation for the spot distribution for stars other than the Sun. The above result also implies that our basic conclusion is not so sensitive to the choice of Ntot; see figure 13.

Same as figure 7, but computed for the top 10 spots (upper two panels) and the remaining 20 spots (lower two panels).
Fig. 11.

Same as figure 7, but computed for the top 10 spots (upper two panels) and the remaining 20 spots (lower two panels).

Same as figure 9, but computed for the top 10 spots (upper two panels) and the remaining 20 spots (lower two panels).
Fig. 12.

Same as figure 9, but computed for the top 10 spots (upper two panels) and the remaining 20 spots (lower two panels).

Same as figure 7, but for Ntot = 10.
Fig. 13.

Same as figure 7, but for Ntot = 10.

3.3 Extracting the spot signature from the amplitude ratios of harmonics

We have shown that a single spot leaves a distinctive modulation pattern in the Fourier coefficients of the harmonics (figure 6). In principle, the signature is important to distinguish between the true and false rotation periods from photometry. Nevertheless, it may be weakened for more realistic cases of multispots, in particular under the presence of the stellar differential rotation. We consider this question in detail using mock data analysis.

Consider a single spot case. Figure 14 compares the theoretical model predictions (A2/A1)2 (solid lines) for the stellar rotation period against the measurement from the mock data for a single spot located at a given ℓs viewed from a line-of-sight direction of ℓo. We choose the value of bspot ≈ 85 ppm so that it corresponds to the top |$10\%$| of the whole spot distribution, i.e., Aspot = 4.7λ from equations (24) and (23). Incidentally, the theoretical curve is invariant with respect to the transformation of (ℓo, ℓs) ↔ (ℓs, ℓo), as equation (11) indicates.

Ratios of the second to fundamental Fourier coefficients for a single spot. Theoretical predictions (A2/A1)2 for ℓs = −45° (green), −30° (purple), 0° (black), 30° (gray), and 45° (orange) are plotted against ℓo in solid (ℓsℓo ≥ 0) and dashed (ℓsℓo < 0) lines. Open circles and crosses indicate the ratios estimated from the mock light curves for bspot ≈ 85 ppm with σ = 0, and 35 ppm, respectively. Left- and right-hand panels plot the result based on the Fourier power spectrum and the LS analysis with nterms = 2.
Fig. 14.

Ratios of the second to fundamental Fourier coefficients for a single spot. Theoretical predictions (A2/A1)2 for ℓs = −45° (green), −30° (purple), 0° (black), 30° (gray), and 45° (orange) are plotted against ℓo in solid (ℓso ≥ 0) and dashed (ℓso < 0) lines. Open circles and crosses indicate the ratios estimated from the mock light curves for bspot ≈ 85 ppm with σ = 0, and 35 ppm, respectively. Left- and right-hand panels plot the result based on the Fourier power spectrum and the LS analysis with nterms = 2.

We adopt two different estimators. One is based on the standard Fourier power spectrum, and plots the corresponding amplitude ratio P2/P1 (left-hand panel). The other is based on the LS analysis. In this case, we first identify the best-fitting angular frequency ωfit from the LS power spectra using nterms = 1 in equation (31). We then fit the data to equation (31) with nterms = 2 by setting ω = ωfit, and obtain the Fourier coefficients S1, S2, S3, and S4 simultaneously. The symbols in the right-hand panel plot the ratio |$(S_3^2+S_4^2)/(S_1^2+S_2^2)$|⁠.

In the noiseless case, the Fourier power spectrum recovers the theoretical predictions very well, but the LS analysis seems slightly but systematically to underestimate the theoretical values. We do not understand why, but the fit to equation (31) with nterms = 2 might be too restrictive and thus very sensitive to the best-fitting value of ωfit estimated from that with nterms = 1.

In any case, the ratios estimated for data with σn = 35 ppm are not so accurate, especially when the latitude of the spot is significantly different from the observer’s line-of-sight (with different signs of ℓs and ℓo, for instance). Therefore, figure 14 implies that it is possible to constrain ℓs and ℓo from the harmonic amplitude ratio for a single spot at least for σn = 35 ppm. For multi-spot cases, however, we find that the amplitude ratio varies significantly due to the differential rotation. Thus this methodology seems to be useful to constrain the spot parameter only when the photometric signal is dominated by a single prominent region.

3.4 Photometric rotation period for differentially rotating stars

Time-dependent distribution of multi-spots over a stellar surface leads to complex photometric modulation signals. Combined with the effect of latitudinal differential rotation, the peak of the rotation period would vary at different observing epochs. In turn, the variation of the rotation period among different quarters may constrain the degree of the differential rotation.

In order to examine to what extent such signatures are indeed detectable from the Kepler data, we perform the LS analysis for seven different sets of mock data and plot distribution of the peak rotation period Prot,LS in figure 15. Basically, we adopt the fiducial values for parameters in table 1; the equatorial rotation period of 10 d, Ntot = 30 spots following the Weibull distribution with k = 0.54 and λ = 2.88 μHem, and the cadence of Tsamp = 30 min over an observing period of Tobs = 90 d corresponding to one quarter of the Kepler long-cadence data. The latitude of the observer’s line-of-sight ℓo and the range of the spot latitude ℓs,max are indicated in each panel.

Distribution of the photometric rotation period and the Fourier coefficient ratio estimated from 300 realizations of Ntot = 30 spots with σn = 35 ppm. We adopt the fiducial parameter set of the equatorial rotation period 10 d, Tsamp = 30 min, Tobs = 90 d. We compute the histograms by varying the observer’s line-of-sight, and the range of the spot latitudes −ℓs,max < ℓs < ℓs,max, and the differential rotation parameters, which are indicated in each panel. Vertical dotted lines indicate the range of the differentially rotation periods as in figure 7 (left-hand panels), and the predicted ratio for ℓs,max (right-hand panels). Panels (a), (b) and (c) use the same set of 300 realizations but viewed from the different observer’s latitudes.
Fig. 15.

Distribution of the photometric rotation period and the Fourier coefficient ratio estimated from 300 realizations of Ntot = 30 spots with σn = 35 ppm. We adopt the fiducial parameter set of the equatorial rotation period 10 d, Tsamp = 30 min, Tobs = 90 d. We compute the histograms by varying the observer’s line-of-sight, and the range of the spot latitudes −ℓs,max < ℓs < ℓs,max, and the differential rotation parameters, which are indicated in each panel. Vertical dotted lines indicate the range of the differentially rotation periods as in figure 7 (left-hand panels), and the predicted ratio for ℓs,max (right-hand panels). Panels (a), (b) and (c) use the same set of 300 realizations but viewed from the different observer’s latitudes.

The left-hand panel of figure 15 shows the histograms of the identified rotation period Prot,LS, while the right-hand panel plots histograms of the corresponding harmonic amplitude ratio (A2/A1)2. Each histogram for seven models is computed from 300 realizations. The first three panels (a), (b), and (c) assume a differential rotation law and spot pattern similar to the Sun. They use the same 300 realizations of the spot pattern over −30° ≤ ℓs ≤ 30°, but viewed from ℓo = 0°, 45°, and 75°, respectively. Similarly, panels (d) and (e) share the same set of 300 realizations with −75° ≤ ℓs ≤ 75° but viewed from ℓo = 0° and 45°, respectively.

According to equation (2), the rotation period of the surface is longer than its equatorial value (10 d), and the width of the distribution reflects the observed range of the spot latitudes ℓs and the values of α2 and α4.

Any difference among panels (a), (b), and (c) is simply due to the fact that the observer at higher ℓo preferentially sees the spots located at higher ℓs, as clearly illustrated in figures 3, 4, and 5. Since the rotation periods estimated by observers at high ℓo should be dominated by a small number of big spots around ℓs > 0°, their distribution is shifted towards the larger Prot,LS due to the differential rotation, and the corresponding amplitude ratio becomes smaller as qualitatively expected from figure 14. A fraction of spot patterns may exhibit an approximate symmetry between ϕs and ϕs + π by chance, which would be interpreted as Prot,LS = 5 d. Such symmetric patterns are more likely to be visible from the edge-on view (ℓo = 0°), which explains the fraction of the second peak around Prot,LS = 5 d in panels (a), (b), and (c).

The next two panels (d) and (e) consider the case for the broader spot distribution over −75° ≤ ℓs ≤ 75°. Because of the presence of a few spots located at higher latitudes, the differential rotation becomes more important, and the distribution of Prot, LS becomes even broader towards its larger value.

The last two panels are shown just for comparison purpose; panel (f) is for the stronger differential rotation case (α2 = 3α2⊙ and α4 = 3α4⊙), and panel (g) is for rigid rotation. Given the same spot distribution pattern, comparison among panels (a), (f), and (g) indicates how the differential rotation law affects the distribution of the rotation period at different quarters of the Kepler data, for instance. This is expected to be directly applicable to put statistical constraints on the degree of latitudinal differential rotation of a population of stars, or to estimate the parameter α2 (and even α4) for stars exhibiting clear photometric light-curve modulations.

While the harmonic ratios shown in the right-hand pane reflect the statistical distribution of the spot latitudes to some extent, they are sensitive to the spot area distribution and do not seem to provide quantitatively useful information. Nevertheless, the histograms are qualitatively consistent with the expected range of the ratios plotted as the vertical dotted lines.

4 Summary and conclusion

We have presented an analytic model of the light-curve variation due to starspots on a differentially rotating surface. If the dynamics of the spots over the timescale of the observing period is neglected, the Fourier coefficients of the harmonics of the rotation period are written primarily in terms of the latitude of spots and the observer’s line-of-sight direction angle.

In order to understand the resulting light-curve variations, we generate various realizations of starspots according to the analytic model, and compute the LS power spectra for the mock datasets.

Even though our analytical model neglects the evolution of spots on the stellar surface (dynamical motion, creation and annihilation), its prediction provides a useful framework to interpret the photometric variation of stars, in particular from the existing Kepler data and the future space-borne missions. The conclusion and implications of this paper are summarized below.

  1. If a photometric light curve of a star exhibits a clear single peak in the LS periodogram, the star may be well approximated as a rigid rotator, and the peak should correspond to the rotation period.

  2. For those stars that have multiple peaks in the LS periodogram, the distribution of the peaks estimated in different quarters may be used to put constraints on parameters characterizing the differential rotation law.

  3. In principle, the ratio of harmonics for the rotation period may constrain the spot latitude ℓs and stellar inclination ℓo given a limb-darkening law. The constraint, however, is sensitive to the spot distribution, and seems to be useful only for a single-spot dominated case. Nevertheless, joint analysis with independent constraints on ℓo from spectroscopic and/or asteroseismic measurements may improve the constraint. We have not explored this possibility in the present paper, but it is worthwhile to pursue in future.

  4. The analytical model presented in the paper is based on the distribution of the effective area of spots bspot alone, and does not require the information of the geometric area and temperature simultaneously. Thus it is applicable not only for spots on main-sequence stars, but for other inhomogeneities on rotating systems. For instance, the recent discovery of the fastest-period white dwarf (Kilic et al. 2021) indicates that the interpretation of the photometric modulation of white dwarfs is crucial in extracting their rotation period. Since it is likely to originate from the small hot spot around the polar region, the determination of ℓs and ℓo with respect to our line-of-sight may be more promising for white dwarfs than for stars with many different spots, as long as the modulation signal-to-noise ratio is sufficiently high.

The above findings may have numerous useful applications even in the existing Kepler data that cover a wide variety of stars with different properties of spots on their surfaces. We are currently working on the joint analysis of photometric and asteroseismic measurements of Kepler stars selected by Kamiaka, Benomar, and Suto (2018), and plan to present the results elsewhere in due course (Y. Lu et al. in preparation).

Acknowledgements

We thank an anonymous referee for various constructive comments on the manuscript. Simulations and analyses in this paper made use of a community-developed core Python package for Astronomy, Astropy. This work is supported by Grants-in Aid for Scientific Research by the Japan Society for Promotion of Science (JSPS) No.18H012 and No.19H01947, and from JSPS Core-to-core Program “International Network of Planetary Sciences”.

Appendix 1 Fourier series expansion for spots with |Γ| < 1

The present paper is based on the analytic Fourier series expansion of the photometric light curve due to a single spot. Those expressions are derived in this appendix.

We first compute the Fourier expansion of equation (8) by setting
(A1)
Thus, the coefficients an are simply given by
(A2)
It is convenient to introduce the angle θc for |Γ| < 1 through
(A3)
Then, equation (A2) reduces to
(A4)
The straightforward integration of equation (A4) yields
(A5)
 
(A6)
 
(A7)
which are a set of coefficients shown in equation (12) of the main text.
If the limb-darkening effect is considered, one has to compute two additional expansions, including
(A8)
and
(A9)
Similarly to equation (A4), the coefficients bn and cn are given as
(A10)
and
(A11)
After tedious but straightforward calculations, we obtain
(A12)
 
(A13)
 
(A14)
 
(A15)
and
(A16)
 
(A17)
 
(A18)
 
(A19)
 
(A20)
The above coefficients are combined and form the coefficients An in equation (18) of the main text:
(A21)

Appendix 2 Fourier series expansion for spots with Γ > 1

For those spots with Γ > 1, the coefficients An for spots with |Γ| < 1 should be replaced by |$\tilde{A}_n$|⁠, which are defined through
(A22)
where
(A23)
 
(A24)
Unlike in appendix 1, the left-hand-side of equation (A22) is explicitly written in terms of up to the third-order polynomials of cos ωst. Thus, |$\tilde{A}_n$| can be explicitly given in the following forms:
(A25)
 
(A26)
 
(A27)
 
(A28)
 
(A29)

References

Albrecht
 
S. H.
,
Marcussen
 
M. L.
,
Winn
 
J. N.
,
Dawson
 
R. I.
,
Knudstrup
 
E.
 
2021
,
ApJ
,
916
,
L1

Angus
 
R.
,
Morton
 
T.
,
Aigrain
 
S.
,
Foreman-Mackey
 
D.
,
Rajpaul
 
V.
 
2018
,
MNRAS
,
474
,
2094

Appourchaux
 
T.
 et al.  
2008
,
A&A
,
488
,
705

Baglin
 
A.
,
Auvergne
 
M.
,
Barge
 
P.
,
Deleuil
 
M.
,
Catala
 
C.
,
Michel
 
E.
,
Weiss
 
W.
 
COROT Team
 
2006a
, in
The CoRoT Mission Pre-Launch Status Stellar: Seismology and Planet Finding
, ed.
Fridlund
 
M.
 et al.
,
ESA SP-1306
(
Noordwijk
:
ESA
),
33

Barnes
 
J. R.
,
Collier Cameron
 
A.
,
Donati
 
J.-F.
,
James
 
D. J.
,
Marsden
 
S. C.
,
Petit
 
P.
 
2005
,
MNRAS
,
357
,
L1

Basri
 
G.
,
Shah
 
R.
 
2020
,
ApJ
,
901
,
14

Benomar
 
O.
 et al.  
2018
,
Science
,
361
,
1231

Benomar
 
O.
,
Masuda
 
K.
,
Shibahashi
 
H.
,
Suto
 
Y.
 
2014
,
PASJ
,
66
,
94

Borucki
 
W. J.
 et al.  
2010
,
Science
,
327
,
977

Brun
 
A. S.
 et al.  
2017
,
ApJ
,
836
,
192

Chaplin
 
W. J.
 et al.  
2010
,
ApJ
,
713
,
L169

Christensen-Dalsgaard
 
J.
 et al.  
1996
,
Science
,
272
,
1286

Cox
 
A. N.
ed.
2000
,
Allen’s Astrophysical Quantities
, 4th ed. (
New York
:
Springer
)

Donati
 
J.-F.
 et al.  
2010
,
MNRAS
,
402
,
1426

Donati
 
J.-F.
,
Collier Cameron
 
A.
 
1997
,
MNRAS
,
291
,
1

Dorren
 
J. D.
 
1987
,
ApJ
,
320
,
756

Farr
 
B.
,
Farr
 
W. M.
,
Cowan
 
N. B.
,
Haggard
 
H. M.
,
Robinson
 
T.
 
2018
,
AJ
,
156
,
146

Fujii
 
Y.
,
Kawahara
 
H.
,
Suto
 
Y.
,
Fukuda
 
S.
,
Nakajima
 
T.
,
Livengood
 
T. A.
,
Turner
 
E. L.
 
2011
,
ApJ
,
738
,
184

Fujii
 
Y.
,
Kawahara
 
H.
,
Suto
 
Y.
,
Taruya
 
A.
,
Fukuda
 
S.
,
Nakajima
 
T.
,
Turner
 
E. L.
 
2010
,
ApJ
,
715
,
866

Haggard
 
H. M.
,
Cowan
 
N. B.
 
2018
,
MNRAS
,
478
,
371

Hirano
 
T.
,
Sanchis-Ojeda
 
R.
,
Takeda
 
Y.
,
Narita
 
N.
,
Winn
 
J. N.
,
Taruya
 
A.
,
Suto
 
Y.
 
2012
,
ApJ
,
756
,
66

Huber
 
D.
 et al.  
2013
,
ApJ
,
767
,
127

Kamiaka
 
S.
,
Benomar
 
O.
,
Suto
 
Y.
 
2018
,
MNRAS
,
479
,
391

Kamiaka
 
S.
,
Benomar
 
O.
,
Suto
 
Y.
,
Dai
 
F.
,
Masuda
 
K.
,
Winn
 
J. N.
 
2019
,
AJ
,
157
,
137

Kawaler
 
S. D.
 
1988
,
ApJ
,
333
,
236

Kilic
 
M.
,
Kosakowski
 
A.
,
Moss
 
A. G.
,
Bergeron
 
P.
,
Conly
 
A. A.
 
2021
,
ApJ
,
923
,
L6

Kjeldsen
 
H.
,
Bedding
 
T. R.
 
1995
,
A&A
,
293
,
87

Lomb
 
N. R.
 
1976
,
Astrophys. Space Sci.
,
39
,
447

Louden
 
E. M.
,
Winn
 
J. N.
,
Petigura
 
E. A.
,
Isaacson
 
H.
,
Howard
 
A. W.
,
Masuda
 
K.
,
Albrecht
 
S.
,
Kosiarek
 
M. R.
 
2021
,
AJ
,
161
,
68

Lund
 
M. N.
 et al.  
2017
,
ApJ
,
835
,
172

MacGregor
 
K.
,
Brenner
 
M.
 
1991
,
ApJ
,
376
,
204

McQuillan
 
A.
,
Mazeh
 
T.
,
Aigrain
 
S.
 
2014
,
ApJS
,
211
,
24

Mandal
 
S.
,
Krivova
 
N. A.
,
Cameron
 
R.
,
Solanki
 
S. K.
 
2021
,
A&A
,
652
,
A9

Maunder
 
E. W.
 
1904
,
MNRAS
,
64
,
747

Mazeh
 
T.
,
Perets
 
H. B.
,
McQuillan
 
A.
,
Goldstein
 
E. S.
 
2015
,
ApJ
,
801
,
3

Morris
 
B. M.
 
2020
,
ApJ
,
893
,
67

Muñoz-Jaramillo
 
A.
 et al.  
2015
,
ApJ
,
800
,
48

Nakagawa
 
Y.
 et al.  
2020
,
ApJ
,
898
,
95

Ohta
 
Y.
,
Taruya
 
A.
,
Suto
 
Y.
 
2005
,
ApJ
,
622
,
1118

Ossendrijver
 
M.
 
2003
,
A&A Rev.
,
11
,
287

Queloz
 
D.
,
Eggenberger
 
A.
,
Mayor
 
M.
,
Perrier
 
C.
,
Beuzit
 
J. L.
,
Naef
 
D.
,
Sivan
 
J. P.
,
Udry
 
S.
 
2000
,
A&A
,
359
,
L13

Ricker
 
G. R.
 et al.  
2014
, in
Proc. SPIE 9143, Space Telescopes and Instrumentation 2014: Optical, Infrared, and Millimeter Wave
, ed.
Oschmann
 
J. M.
 Jr.
 et al.
, (
Bellingham, WA
:
SPIE
),
914320

Roettenbacher
 
R. M.
,
Monnier
 
J. D.
,
Harmon
 
R. O.
,
Barclay
 
T.
,
Still
 
M.
 
2013
,
ApJ
,
767
,
60

Sanchis-Ojeda
 
R.
,
Winn
 
J. N.
 
2011
,
ApJ
,
743
,
61

Sanchis-Ojeda
 
R.
,
Winn
 
J. N.
,
Holman
 
M. J.
,
Carter
 
J. A.
,
Osip
 
D. J.
,
Fuentes
 
C. I.
 
2011
,
ApJ
,
733
,
127

Sasaki
 
S.
,
Suto
 
Y.
 
2021
,
PASJ
,
73
,
1656

Scargle
 
J. D.
 
1982
,
ApJ
,
263
,
835

Skumanich
 
A.
 
1972
,
ApJ
,
171
,
565

Snodgrass
 
H. B.
,
Ulrich
 
R. K.
 
1990
,
ApJ
,
351
,
309

van Saders
 
J.
,
Ceillier
 
T.
,
Metcalfe
 
T.
,
Silva Aguirre
 
V.
,
Pinsonneault
 
M.
,
García
 
R.
,
Mathur
 
S.
,
Davies
 
G.
 
2016
,
Nature
,
529
,
181

Varela
 
J.
,
Strugarek
 
A.
,
Brun
 
A.
 
2016
,
Adv. Space Res.
,
58
,
1507

Walker
 
G.
 et al.  
2003
,
PASP
,
115
,
1023

Walkowicz
 
L. M.
,
Basri
 
G.
,
Valenti
 
J. A.
 
2013
,
ApJS
,
205
,
17

Zharkov
 
S.
,
Zharkova
 
V. V.
,
Ipson
 
S. S.
 
2005
,
Sol. Phys.
,
228
,
377

This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://dbpia.nl.go.kr/journals/pages/open_access/funder_policies/chorus/standard_publication_model)