A STUDY OF DIFFERENTIAL ROTATION ON II PEGASI VIA PHOTOMETRIC STARSPOT IMAGING

II Pegasi is an SB1 binary system for which the primary component was determined to be of spectral class K2-3 IV–V by Rucinski (1977). It was classified as an RS CVn system by Vogt (1981a). The first photometric light curves were obtained by Chugainov (1976), who found variability with a period of approximately 6.75 days and interpreted the asymmetric light curve in terms of rotational modulation due to large, cool starspots. In 1986 September, the difference between maximum and minimum light for the V filter reached 0.5 mag, implying a projected spot area coverage of the visible hemisphere at minimum light on the order of 50% (Doyle et al. 1989).

On the basis of high-quality radial velocity measurements, Berdyugina et al. (1998b) determined the revolution period of the binary to be 6.724333 ± 0.000010 days. The same authors performed a detailed model atmosphere analysis of high-resolution and high signal-to-noise CCD spectra, obtaining values for the photospheric temperature and surface gravity of the primary star of T_eff = 4600 K and log g = 3.2, with g expressed in cgs units. These values correspond to a K2 IV star of mass M = 0.8 ± 0.1 M_☉. They estimated the radius of the primary as R = 3.4 ± 0.2 R_☉ and the inclination to be i = 60° ± 10° on the assumption that the rotational axis is perpendicular to the orbital plane. Based on the fact that the secondary star is unseen at all wavelengths and thus has a luminosity at least 100 times smaller than that of the primary, they estimated the secondary to be an M0–M3 red dwarf.

II Peg is among the most active RS CVn systems, and it is one of a small number of binaries in which the Hα line is always seen in emission (Nations & Ramsey 1981). Recently, Frasca et al. (2008) reported on contemporaneous photometric and spectroscopic observations of II Peg, finding that the Hα emission and photometric intensity are strongly anticorrelated, suggesting that regions of high chromospheric activity are physically associated with the spots. This conclusion was corroborated by a rotational modulation of the intensity of the He i D₃ line. Based on an estimated radius of R = 2.76 R_☉ and vsin i = 22.6kms⁻¹, they estimated the inclination between the rotation axis and the line of sight to be i = 60°⁺³⁰₋₁₀. Messina (2008) confirmed via long-term monitoring of V as well as the B − V and U − B colors that II Peg is redder when it is dimmer, as would be expected if the dimming is caused by cool spots.

A number of studies have attempted to determine the spot temperatures. Vogt (1981b) modeled light and color curves obtained in 1977 with a single circular spot, finding a spot temperature of T_spot = 3400 ± 100K. Nations & Ramsey (1981) obtained T_spot = 3600K from observations in the Fall of 1979, Poe & Eaton (1985) obtained T_spot = 3620K for Fall 1980, Rodonò et al. (1986) obtained T_spot = 3300K for Fall 1981, Byrne & Marang (1987) obtained T_spot = 3700K for Fall 1986, and Boyd et al. (1987) obtained T_spot = 3450K for 1986–1987. By modeling the strengths of TiO absorption bands, O'Neal et al. (1998) found evidence for multiple spot temperatures, finding T_spot to vary between 3350 ± 60K and 3550 ± 70K as the star was observed through slightly less than one rotational period. More recently, from a spot model applied to contemporaneous photometry and spectroscopy, Frasca et al. (2008) obtained T_spot ≈ 3600K.

Henry et al. (1995) used a simple analytic two-spot model to fit photoelectric light curves of four chromospherically active binaries: λ And, σ Gem, V711 Tau, and II Peg. The II Peg data were acquired from 1973 to 1992 and subdivided into 37 individual light curves. They plotted "migration curves" obtained via the spot-model curve fits for 12 long-lived spots they identified in the data from the times of minimum light. A migration curve shows the variation in the phase of minimum light with time, where the phase is computed using the orbital ephemeris and represents the fractional part of the number of rotation periods since an arbitrary starting time. Assuming tidal locking and no differential rotation, a given spot would always cross the central meridian of the stellar disk as seen from Earth once per revolution period and thus always at the same phase. However, if the star exhibits latitude-dependent differential rotation, we would expect to see a given spot progressively advanced or retarded in phase relative to the orbital ephemeris. A plot of the phase of minimum light versus time for a given spot should then be a straight line with slope determined by the difference in the rotation period of the latitude of the spot and the revolution period. This was precisely what Henry et al. observed in the data, for II Peg and the other stars. The plots for different spots had different slopes, demonstrating latitude-dependent differential rotation. The degree of differential rotation was specified in terms of the differential rotation coefficient, k, defined for the Sun by fitting the rotation period as a function of latitude with the relation

$$\begin{equation} P(\theta) = \frac{P_\mathrm{eq}}{1 - k\sin ^2\theta }, \end{equation} \tag{ 1 }$$

where P(θ) is the rotation period for latitude θ and P_eq is the rotation period at the equator. For the Sun, k = 0.19. If the differential rotation of other stars has the same functional form as for the Sun, and if the rotation periods for spots sampling a range of latitudes are determined for a star, then the coefficient k is given by

$$\begin{equation} \frac{P_\mathrm{max} - P_\mathrm{min}}{P_\mathrm{avg}} = kf, \end{equation} \tag{ 2 }$$

where P_max, P_min, and P_avg are the maximum, minimum, and average observed periods, and f is a distribution function which relates the total range in rotational period sampled to the number of spots for which the period has been determined (Hall & Henry 1994). The value of f ranges from 0.5 for two spots to over 0.9 when the number of spots exceeds 6. Henry et al. (1995) used 8 of the 12 spots they observed on II Peg (four spots were observed over intervals too short to allow their periods to be obtained reliably) to determine k using Equation (2), with the result k = 0.005 ± 0.001.

Siwak et al. (2010) obtained k = 0.0245^+0.0155_−0.0020 by fitting a model with two differentially rotating circular spots to data obtained by the MOST satellite on 31 consecutive days in 2008 September–October. This is nearly five times greater than the value found by Henry et al., while the results of the present study agree better with those of Henry et al. We defer further discussion of these discrepancies to Section 5.

Rodonò et al. (2000) performed an analysis similar to the present study, inverting light curves acquired between 1974 and 1998 to produce images of the stellar surface. In contrast to the smoothing function used here (see Section 2 and in particular Equation (12)), they used maximum entropy and Tikhonov regularization. They concluded that the distribution of spots on II Peg consists of a component distributed uniformly in longitude which does not rotationally modulate the light curve (but does produce a secular variation in the mean intensity), plus an unevenly distributed component responsible for the rotational modulation. Their analysis indicated that the uniformly distributed component varied in total area with a period of ∼13.5yr. They determined the unevenly distributed component to be concentrated around three active longitudes, one of these having an essentially permanent presence but a cycle in spot area with period ∼9.5yr. They found the activity of the other two active longitudes to switch back and forth, with one active while the other is inactive, with a period of ∼6.8yr. However, there is an interval of ∼1.05yr before the switch in which both longitudes are active. There is thus a period of ∼6.8 − 2(1.05) = 4.7yr during which only one of the two longitudes engaged in the "flip-flop" behavior is active, which agrees with the switching period deduced by Berdyugina & Tuominen (1998) from the times of light minima.

From a periodogram analysis, Henry et al. (1995) found periodicities in the mean magnitudes for the spot-model fits of their 37 light curves of 4.4 ± 0.2yr and 11 ± 2yr. They interpreted the 4.4 yr period as reflecting the average lifetime of the spots and the 11 yr period as representing a different timescale. Rodonò et al. interpreted the 4.7 yr periodicity arising in their analysis as corresponding to the 4.4 yr period obtained by Henry et al., while they interpreted the 9.5 yr period they saw in the total area of the spot component which is unevenly distributed in longitude as corresponding to the 11 yr period found by Henry et al.

Henry et al. noted that their two-spot model, which assumed circular spots varying only in radius over time, was not fully adequate to explain the variations with time of the II Peg light curve. In particular, when the amplitude of the rotational modulation due to spots they designated G and H was diminishing, the mean brightness of the star stayed roughly constant. Similar behavior was seen for another pair of spots, which they designated J and K. If the decrease in amplitude were due simply to a decrease in the spot radii, the mean brightness of the star should have increased (assuming no change in the brightness of the photosphere outside the spots). On the other hand, if instead the spots were being drawn out in longitude by differential rotation while maintaining nearly constant area, then the amplitude would decrease while the mean brightness stayed constant, as observed.

The present study is suited to look for evidence for such drawing out (or compression) in longitude of active regions by differential rotation, as we produce images of the active regions and observe changes in them over periods of several months. By simultaneously inverting contemporaneous B and V light curves, we exploit differences in the limb darkening as seen through the two filters to achieve significantly better latitude resolution than is possible when using light curves obtained through only a single filter (Harmon & Crews 2000), thereby allowing us to directly detect differential rotation in our images. It should be noted, however, that we do not claim to obtain accurate spot latitudes from our two-filter inversions; nonetheless, simulations like those detailed in Harmon & Crews show that relative spot latitudes can be obtained with good reliability, i.e., when two spots are present, the one at the lower latitude is rendered as such. In this regard our approach differs from that of Rodonò et al., who used their V-filter surface imagery just to derive quantities that they claim are independent of the regularization criterion, such as the distributions of the spots versus longitude, the changes in the distribution over time, and the variations of the total area covered by spots.

In Section 2, we discuss the method used to invert the light curves. In Section 3, we discuss the division of the over 22 years worth of B and V light curves into separate data sets and the procedure used to process them for inversion. In Section 4, we discuss in detail the results for the six intervals over which we saw good evidence for latitude-dependent differential rotation. Finally, in Section 5, we present our results for the value of the differential rotation coefficient k for II Peg and compare our results to those of Henry et al. (1995) and Siwak et al. (2010).

Light-curve inversion (LI) is a photometric imaging technique which generates a map of a star's surface based on the brightness variations produced as dark (or possibly bright) starspots are carried into and out of view of Earth by the star's rotation. It makes no a priori assumptions regarding the number of spots on the surface or their shapes. The details regarding the implementation of the algorithm are presented in Harmon & Crews (2000), along with the results of extensive tests in which artificial stellar surfaces were used to create light curves, which were then inverted. In Harmon & Crews, the technique is called "Matrix Light-curve Inversion (MLI)," because it evolved from the original formulation described in Wild (1989) and called by that name. However, because the formulation described in Harmon & Crews (2000) and as modified in the present work no longer uses matrices, we shortened the name (while of course fully realizing that many techniques can be categorized as LI algorithms). Here, we outline the method and refer the reader to Harmon & Crews for more details. The reader may note that what we are doing is using an alternative regularization criterion as opposed to the commonly used maximum entropy or Tikhonov regularization described in, e.g., Craig & Brown (1986).

We assume the star to be spherical; no attempt is made to account for ellipsoidal distortion or reflection effects due to a companion star (though it would be possible to modify the method to do so). The stellar surface is subdivided into N bands in latitude of equal angular widths Δθ = π/N. Each latitude band is further subdivided into patches which are all "spherical rectangles" of equal widths in longitude Δϕ = 2π/M_i, where M_i is the number of patches in the ith latitude band. The M_i are chosen to be proportional to the cosine of the latitude (but with the constraint that the M_i must be integers) so that the areas of all the patches are nearly equal. The visible pole is defined to be the north pole, with latitude θ = +90°, while the hidden south pole has latitude θ = −90°. The jth patch in the ith latitude band is designated patch (i, j). The first patch in each latitude band, patch (i, 1), straddles the meridian with longitude ϕ = 0, defined to be the one which intersects the equator on the approaching limb of the star at an arbitrarily chosen reference time t₀. Longitude increases in the direction of the star's rotation, so that the sub-observer longitude at t = t₀ is ϕ = 90°. In the absence of interstellar absorption, at the time t_nk of observation number k through filter n, the intensity I_nk observed at Earth is (in the limit that the number of patches is large)

$$\begin{eqnarray} I_{nk}\index {} &=& \sum _{i=1}^{N_s}\sum _{j=1}^{M_i} \Omega _{nk;ij}L_{nk;ij}J_{n;ij}, \quad n = 1, \ldots, Q,\nonumber\\ k &=& 1, \ldots, P_n, \end{eqnarray} \tag{ 3 }$$

where Q is the number of filters, P_n is the number of observations through filter n, J_n;ij is the specific intensity (W m⁻² sr⁻¹) along the outward normal of patch (i, j) integrated over the passband of filter n, Ω_nk;ij is the solid angle of patch (i, j) as seen from Earth at time t_nk (we set Ω_nk;ij = 0 if the patch is on the far side of the star), L_nk;ij is the factor by which the specific intensity emitted in the direction of Earth is attenuated by limb darkening compared to that emitted along the outward normal (so that L_nk;ijJ_n;ij is the specific intensity emitted along the line of sight to Earth), and N_s is the index of the southernmost latitude band which is visible from Earth.

The goal of LI is to find a set of computed patch intensities $\hat{J}_{n;ij}$ that mimics the actual variations of surface brightness across the stellar surface as closely as possible. (We use a caret over a quantity to indicate that it represents a value as computed by the LI algorithm.) Since we generally do not know the actual radius and distance of the star very precisely, we content ourselves with finding only the relative brightnesses of the patches. To this end, we simply define the radius of the star to be 1 and use the area of a patch projected onto the plane of the sky as a proxy for the solid angle it subtends at Earth.

We use the limb-darkening coefficients published by Van Hamme (1993) to determine the values of the L_nk;ij in Equation (3). The benefit of observing through multiple filters is that we can take advantage of the differences in the degree of limb darkening as seen through different filters in order to significantly increase the latitude resolution of the inversions, as explained in Harmon & Crews (2000). In order to take advantage of this information, we must simultaneously invert all of the filter light curves. This in turn requires that we couple together the $\hat{J}_{n;ij}$ for different values of the filter index, n. To do this, we designate the filter for which the light curve has the lowest noise as the "primary filter" and assign it filter index n = 1. For simplicity we assume that the actual stellar surface can be described via a two-component model in which all the spots have the same temperature T_spot and thus emit the same specific intensity along the outward normal as seen through filter n, which we designate at J_n;spot; similarly, we assume that all points on the surface outside spots are part of a photosphere of uniform temperature T_phot and emitting specific intensity J_n;phot along the outward normal. However, it should be noted that the reconstructed surface created by the inversion of the data does not have the property that the $\hat{J}_{n;ij}$ can have only one of two values; they are continuous variables. We then define the intensities of the patches as viewed through filter n ≠ 1 via the linear scaling

$$\begin{equation} \hat{J}_{n;ij} \equiv \frac{r_n}{1 - s_1}[(s_n - s_1)\hat{J}_{1;\mathrm{avg}} + (1 - s_n)\hat{J}_{1;ij}]. \end{equation} \tag{ 4 }$$

Here $\hat{J}_{1;\mathrm{avg}}$ is the average value of the $\hat{J}_{1;{ij}}$, r_n is the estimated value of J_n;phot/J_1;phot on the actual stellar surface, and s_n is the estimated value of J_n;spot/J_n;phot for the actual stellar surface. We estimate these ratios by calculating the Planck function at the central wavelength of the filter in question at the assumed spot and photosphere temperatures. It would be more accurate to integrate the product of the Planck function and the filter sensitivity functions over their passbands, but typically we do not know the spot and photosphere temperatures with sufficient precision to justify the extra effort. The scaling from $\hat{J}_{1;ij}$ to $\hat{J}_{n;ij}$ given by Equation (4) has the property that when $\hat{J}_{1;ij}/\hat{J}_{1;\mathrm{avg}} = J_{1;\mathrm{spot}}/J_{1;\mathrm{phot}}$ according to our estimate, then $\hat{J}_{n;ij}/\hat{J}_{n;\mathrm{avg}} = J_{n;\mathrm{spot}}/J_{n;\mathrm{phot}}$ as well. We are using $\hat{J}_{1;\mathrm{avg}}$ as a proxy for J_1;phot, which should be a reasonable approximation as long as the stellar photosphere comprises most of the surface.

Since the patch brightnesses through all the filters in the model are entirely determined in terms of their brightnesses $\hat{J}_{1;ij}$ as seen through filter 1, the problem reduces to finding these values, so for notational simplicity we define $\hat{J}_{ij} \equiv \hat{J}_{1;ij}$. As is well known, the problem of determining the $\hat{J}_{ij}$ is extremely sensitive to the presence of even small amounts of noise in the data. This can be seen by considering the effect on the light curve produced by a myriad of small spots distributed all over the surface. As the star rotated, at any given time nearly equal numbers of spots would be rising over the approaching limb and setting over the receding limb. Thus, the total contribution to the star's brightness from the spots would be nearly but not exactly constant, so that the effect of the spots would be to impart a small-amplitude, high-frequency ripple on the light curve, very similar to the effect of random noise in the observations. Conversely, if we attempt to fit noisy data, then unless precautions are taken, the resulting model surface will be covered with spurious small spots introduced in order to "explain" the presence of the noise in the signal.

To avoid this dilemma, rather than simply finding the set of $\hat{J}_{ij}$ that yields the best fit to the light curve data, we determine the $\hat{J}_{ij}$ by finding the set of them which minimizes the objective function (Twomey 1977; Craig & Brown 1986):

$$\begin{equation} E(\mathbf {\hat{J}},\mathbf {I},\lambda,B) = G(\mathbf {\hat{J}},\mathbf {I}) + \lambda S(\mathbf {\hat{J}},B). \end{equation} \tag{ 5 }$$

Here $\mathbf {\hat{J}}$ represents the set of the $\hat{J}_{ij}$, while I represents the set of observed intensities I_nk, i.e., the data light curve. The function $G(\mathbf {\hat{J}},\mathbf {I})$ expresses the goodness of fit of the calculated light curve $\mathbf {\hat{I}}$ (with components $\hat{I}_{nk}$) obtained from $\mathbf {\hat{J}}$ to the data light curve I, such that smaller values of $G(\mathbf {\hat{J}},\mathbf {I})$ imply a better fit to the data. The smoothing function $S(\mathbf {\hat{J}},B)$ is defined such that it takes on smaller values for surfaces that are "smoother" in some appropriately defined sense. Finally, λ is an adjustable Lagrange multiplier called the tradeoff parameter, and B is an adjustable parameter called the bias parameter, which is discussed below. Note that as λ → 0, the first term on the right dominates, so that minimizing E is equivalent to minimizing G, and we obtain the solution which best fits the light curve data but suffers from the spurious noise artifacts discussed above. On the other hand, as λ → ∞, the second term dominates, so that minimizing E produces a very "smooth" surface lacking in noise artifacts, but also producing a very poor fit to the data. Thus, by varying λ, we adjust the tradeoff between goodness of fit and smoothness of the model surface. If we choose λ such that $G(\mathbf {\hat{J}},\mathbf {I})$ is equal to a corresponding estimate of the amount of noise in the data, then in a rough sense we can say that by minimizing the objective function, we find the smoothest solution $\mathbf {\hat{J}}$ for which the corresponding light curve $\mathbf {\hat{I}}$ fits the data light curve I to a degree which is as good as but not better than is justified by the noise in the data. In this way, we obtain a model surface which fits the data well, but not so well that it is dominated by noise artifacts.

For the goodness-of-fit function in this study, we use

$$\begin{equation} G(\mathbf {\hat{J}},\mathbf {I}) = \frac{(2.5 \log _{10}e)^2}{P}\sum _{n=1}^Q \frac{1}{\sigma _n^2} \sum _{k=1}^{P_n}\left(\frac{I_{nk}-\hat{I}_{nk}}{I_{nk}}\right)^2. \end{equation} \tag{ 6 }$$

Here we assume that light curves have been obtained through Q different photometric filters (Q = 2 in the present work since we use B- and V-filter data), and that the magnitudes have been converted to intensities. Since our goal is only to find the relative values of the $\hat{J}_{ij}$, it suffices to use relative rather than absolute intensities for the light curve data in the calculation of $G(\mathbf {\hat{J}},\mathbf {I})$. The number of observations in the light curve obtained through filter n is P_n, while P = ∑_nP_n is the total number of data points in all the light curves. The estimated noise variance in the light curve data for filter n expressed in magnitudes is σ²_n. In Harmon & Crews (2000) it is shown that to a good approximation, the true noise variance $\tilde{\sigma }_n^2$ is given by

$$\begin{equation} \tilde{\sigma }_n^2 \approx \frac{(2.5 \log _{10}e)^2}{P_n} \sum _{k=1}^{P_n}\left(\frac{I_{nk}-\tilde{I}_{nk}}{\tilde{I}_{nk}} \right)^2, \end{equation} \tag{ 7 }$$

where $\tilde{I}_{nk}$ is the true noise-free value of the intensity (which of course is not known unless one is doing a simulation). If we define $\epsilon _{nk} \equiv I_{nk} - \tilde{I}_{nk}$ to be the true error in the measurement I_nk, and $\delta _{nk} \equiv \hat{I}_{nk} - \tilde{I}_{nk}$ to be the deviation between the calculated and true intensities, then with this notation

$$\begin{equation} \frac{I_{nk} - \tilde{I}_{nk}}{\tilde{I}_{nk}} = \frac{\epsilon _{nk}}{\tilde{I}_{nk}}, \end{equation} \tag{ 8 }$$

while

$$\begin{equation} \frac{I_{nk} - \hat{I}_{nk}}{I_{nk}} = \frac{\epsilon _{nk} - \delta _{nk}}{\tilde{I}_{nk} + \epsilon _{nk}}. \end{equation} \tag{ 9 }$$

If $\mathbf {\hat{I}}$ is a good match to the data I, then the _nk and δ_nk are small quantities, and we can expand the right-hand side of Equation (9) as

$$\begin{equation} \frac{I_{nk} - \hat{I}_{nk}}{I_{nk}} = \frac{\epsilon _{nk}-\delta _{nk}}{\tilde{I}_{nk}}\left(1 - \frac{\epsilon _{nk}}{\tilde{I}_{nk}} + \cdots \right) = \frac{\epsilon _{nk} - \delta _{nk}}{\tilde{I}_{nk}} + \cdots. \end{equation} \tag{ 10 }$$

Then

$$\begin{eqnarray} G(\mathbf {\hat{J}},\mathbf {I}) &=& \frac{1}{P} \sum _{n=1}^Q \frac{(2.5 \log _{10}e)^2}{\sigma _n^2}\nonumber\\ &&\times\sum _{k=1}^{P_n} \left[\frac{\epsilon _{nk}^2}{\tilde{I}_{nk}^2} + \frac{\delta _{nk}^2 - 2\epsilon _{nk}\delta _{nk}}{\tilde{I}_{nk}^2} + \cdots \right] \nonumber \\ &=& \frac{1}{P}\sum _{n=1}^{P_n}\Bigg[P_n\frac{\tilde{\sigma }_n^2}{\sigma _n^2} + \frac{(2.5\log _{10}e)^2}{\sigma _n^2}\nonumber\\ &&\times\sum _{k=1}^{P_n} \frac{\delta _{nk}^2 - 2\epsilon _{nk}\delta _{nk}}{\tilde{I}_{nk}^2} + \cdots \Bigg]. \end{eqnarray} \tag{ 11 }$$

If the reconstructed intensities $\hat{I}_{nk}$ perfectly matched the true intensities $\tilde{I}_{nk}$, and in addition the estimated noise variances σ_n were equal to the true noise variance $\tilde{\sigma _n}$, then we would have all δ_nk = 0, and so we would have $G(\mathbf {\hat{J}},\mathbf {I}) = 1$ to lowest order in the _nk. For this reason, given the estimates σ_n and for a given value of the bias parameter B, we vary λ until the stopping criterion $G(\mathbf {\hat{J}},\mathbf {I}) = 1$ is attained to a predetermined precision.

The smoothing function used in this study is

$$\begin{equation} S(\mathbf {\hat{J}},B) = \sum _{i=1}^N \sum _{j=1}^{M_i} c_{ij}(\hat{J}_{ij} - \hat{J}_\mathrm{avg})^2, \end{equation} \tag{ 12 }$$

where c_ij = 1 if $\hat{J}_{ij} \le \hat{J}_\mathrm{avg}$, while c_ij = B if $\hat{J}_{ij} > \hat{J}_\mathrm{avg}$. Thus, patches brighter or darker than average incur a penalty in that they increase the value of $S(\mathbf {\hat{J}},B)$ (and thus of the function $E(\mathbf {\hat{J}},\mathbf {I},\lambda,B)$ which is to be minimized) by an increasing amount as the deviation from the average increases. Note that $S(\mathbf {\hat{J}},B)$ satisfies the criterion that it takes on its minimum possible value of zero for a featureless surface which is perfectly "smooth" in that all the patch brightnesses $\hat{J}_{ij}$ are equal, and that surfaces showing greater deviations about the average are judged as "rougher." For B>0, the penalty for a patch being brighter than average by a given amount is B times larger than for a patch darker than average to the same degree. Thus, B biases the solution toward having most patches just slightly brighter than average to represent the stellar photosphere, which is assumed to be almost uniformly bright like the Sun's, while a smaller number are much darker than average to represent the dark starspots. This is the reason for the name "bias parameter."

The simulations described in Harmon & Crews (2000) show that as B is increased, the ratio $\min (\hat{J}_{ij})/\hat{J}_\mathrm{avg}$ decreases, so that the darkest patch becomes darker relative to the average patch brightness. We use this ratio as a proxy for the assumed ratio of the spot and photosphere brightnesses as seen through filter 1. For a given value of the tradeoff parameter λ, the bias parameter B can be varied until $\min (\hat{J}_{ij})/\hat{J}_\mathrm{avg} = s_1$, the estimated spot-to-photosphere brightness ratio seen through filter 1, to a predetermined precision. The scaling given by Equation (4) ensures that $\min (\hat{J}_{n;ij})/\hat{J}_{n;\mathrm{avg}} = s_n$, the estimated spot-to-photosphere brightness ratio for filter n.

The procedure for inverting a series of light curves obtained through a set of filter passbands is then as follows. The input parameters are the estimated noise variances σ²_n of the light curves, the estimated spot and photosphere temperatures T_spot and T_phot, and the inclination angle i of the rotation axis to the line of sight. The spot and photosphere temperatures, T_spot and T_phot, are used to obtain the values of s_n and r_n (including n = 1) in Equation (4), and the inclination angle i is used in finding Ω_nk;ij and L_nk;ij in Equation (3). As described in Harmon & Crews (2000), given estimates of the σ_n, two copies of a root-finding subroutine are used in concert so as to find the values of λ and B such that $G(\mathbf {\hat{J}},\mathbf {I}) = 1$ and $\min (\hat{J}_{n;ij})/\hat{J}_{n;\mathrm{avg}} = s_n$ to the desired precision. The result is a set of solutions, one for each combination of the input parameters. How we select one of these to represent the "best" solution is described in Section 3.

The raw data consisted of Johnson B and V differential magnitudes for the target and comparison stars (Var−Cmp) paired with corresponding heliocentric Modified Julian Dates (MJD), acquired from heliocentric MJD 47115.8086–55222.6238 (1987 November 16–2010 January 26) with the Tennessee State University T3 0.4 m Automated Photometric Telescope at Fairborn Observatory in Arizona (Henry 1995a, 1995b). The complete B and V data sets are plotted at upper and lower left in Figure 1, while at upper and lower right are plotted the difference in the magnitudes of the check and comparison stars (Chk−Cmp). The 435 day gap in the data around JD 2448500 corresponds to the installation and testing of a new photometer with lower noise, which can be seen in the smaller scatter in the Chk−Cmp data after the gap. The standard deviations of the B data before and after the gap are 0.0090 and 0.0058 mag, respectively, while for the V data they are 0.0086 and 0.0062 mag.

Zoom In Zoom Out Reset image size

The first task was to convert the MJDs to rotational phases. The rotational phase Φ is defined as

$$\begin{equation} \Phi (t) = \frac{t - t_0}{T} - \left\lfloor \frac{t - t_0}{T}\right\rfloor, \end{equation} \tag{ 13 }$$

where t is the time of the observation, t₀ is an arbitrary reference time used for all the observations, T is the rotational period of the star, and ⌊x⌋ is the greatest integer which is less than or equal to x. Thus, Φ(t) represents the fraction of a rotation through which the star has turned relative to the orientation it had at time t₀. On the assumption that the star exhibits differential rotation, there is no such thing as the rotational period, so T here represents a suitable average of the rotation period over all latitudes. On the assumption of tidal locking, it is reasonable to use the orbital period for this average. In the present study, T = 6.724333 days was used based on the orbital period obtained by Berdyugina et al. (1998b), and t₀ was chosen as JD 2443033.47, based on the orbital ephemeris of Vogt (1981a), in which t₀ represents the time of superior conjunction, when the primary is farthest from the observer.

The next task was to subdivide the data into individual data sets comprised of pairs of B and V light curves suitable for inversion. Each light curve needed to contain enough data points so as to provide good phase coverage. Ideally this would be achieved using data acquired during a single rotation of the star, since this would minimize the chance that the spot configuration had evolved significantly during the time interval spanned by the data. However, in practice this was not feasible, because the star's rotation period is too long to allow for continuous monitoring during a single rotation, and because the telescope was not dedicated solely to observing II Peg. The desire for good phase coverage is thus in opposition to the desire to minimize the number of stellar rotations covered in any one data set, so some compromises were necessary. If data for different revolutions showed a systematic shift in the magnitudes, it was clear that the stellar surface features had evolved by an unacceptable amount during the interval in question; otherwise, data from additional rotations could be included if needed so as to improve the phase coverage. Several groups of observations were discarded because they were temporally isolated by many rotation periods from the observations nearest to them in time and contained an insufficient number of observations to produce good phase coverage. In the end, 79 pairs of B and V light curves were created that were subsequently inverted to produce the results reported in this study. The mean number of observations per B light curve was 21.6, the median was 20, and the standard deviation was 6.0, while for the V light curves the mean was 21.4, the median was 20, and the standard deviation was 5.7.

The final task before inverting the light curves was to convert the differential magnitudes to relative intensities. The formula used for both the B and the V passbands was simply

$$\begin{equation} I = 10^{-0.4(m-m_0)}, \end{equation} \tag{ 14 }$$

where m was the differential magnitude for the observation in question, and the reference magnitude m₀ was the smallest value of m in the entire data set for the given filter. An intensity of I ≡ 1 was thus assigned to this observation. For the B filter, m₀ = 1.843, while for the V filter, m₀ = 1.522. No attempt was made to calibrate the V intensities relative to the B intensities in an absolute sense, because the LI algorithm neither requires nor would make use of this information, as mentioned in Section 2.

Berdyugina et al. (1998b) obtained i = 60° ± 10° for the inclination of the rotation axis. We thus performed inversions for i = 60° and for i = 45° as well, considering this to be prudent given the considerable uncertainty in the inclination. As will be seen from the discussion of the individual data sets below, the results for both assumed inclinations were generally consistent, increasing our confidence in their validity. While one might argue that we also should have performed inversions for assumed inclinations greater than 60°, say 70° or 75°, simulations like those reported in Harmon & Crews (2000) show that the method works poorly in such circumstances, so we chose not to do so.

For the photosphere temperature we used T_phot = 4600 K based on the work of Berdyugina et al. (1998b). As the spot temperature we chose T_spot = 3500 K, which is compatible with the estimates by other authors mentioned in Section 1. Given the considerable range of estimates of T_spot ≈ 3300–3700 K, we also performed inversions assuming T_spot = 3250 K and 3750 K for all the data sets in Data Group 2, the earliest for which we saw evidence for differential rotation. As discussed further in Section 4.1, we determined that the locations of the spots on the reconstructed surfaced were not significantly affected. As a further check, we also inverted the first data sets in the other data groups for which we saw differential rotation, again finding no significant effect on our results. Based on this, we decided that it sufficed to perform inversions only for T_spot = 3500 K on the rest of the data sets.

The estimates r_n of J_n;phot/J_1;phot and s_n of J_n;spot/J_n;phot appearing in the scaling given by Equation (4) were obtained by evaluating the Planck function describing blackbody radiation at the filter effective wavelengths λ_B,eff = 440 nm and λ_V,eff = 550 nm. For T_phot = 4600 K we find J_B;phot/J_V;phot = 0.734. For this photosphere temperature and T_spot = 3500 K, we obtain s_B = 0.107 and s_V = 0.167, while for T_spot = 3250 K we obtain s_B = 0.0522 and s_V = 0.0939, and for T_spot = 3750 K we find s_B = 0.200 and s_V = 0.275.

The photosphere temperature T_phot and surface acceleration due to gravity g are the input parameters used by Van Hamme (1993) to calculate limb-darkening coefficients based on the ATLAS stellar atmosphere models of Kurucz (1991). We used log g = 3.0, which is appropriate for a K2 subgiant (Gray 1992) and is the value in Van Hamme's tables which is closest to the result log g = 3.2 of Berdyugina et al. (1998b). Van Hamme gives coefficients in steps of 250 K for temperatures in the range 3500K < T_phot < 10,000K, so there is no entry for our value of T_phot = 4600 K. We simply substituted the values corresponding to T_phot = 4500 K, the nearest listed temperature to ours, since we do not know T_phot or g accurately enough to justify interpolating. From Table 2 in Van Hamme's paper, we find for the B filter that the limb darkening is

$$\begin{equation} L_B(\mu) = 1 - \epsilon (1 - \mu) - \delta \mu \ln \mu, \end{equation} \tag{ 15 }$$

with = 0.852 and δ = −0.158, where μ is the cosine of the angle between the outward normal to the surface and the observer's line of sight. For the V filter

$$\begin{equation} L_V(\mu) = 1 - \epsilon (1 - \mu) - \delta (1 - \sqrt{\mu }), \end{equation} \tag{ 16 }$$

with = 0.780 and δ = 0.039. Given the latitude and longitude of patch (i, j) and the angle of inclination of the rotation axis, it is straightforward to compute μ for the center of the patch at any given observation time t_nk and thus to obtain L_nk;ij in Equation (3).

The next step was to obtain estimates of the noise variances σ²_B and σ²_V for each of the 79 light curves for each filter. It is well known that in practice Twomey's criterion of choosing the tradeoff parameter λ so that the variance between the data and the reconstruction is equal to the noise variance leads to oversmoothing (Turchin 1967). In our case this would lead to a loss in resolution of the reconstructed surfaces. However, the technique we use to determine λ avoids this problem, because we obtain an "effective noise level" rather than using an estimate of the noise variance based on the scatter in the comparison star magnitudes. Harmon & Crews (2000) describe simulations in which an artificial star is used to generate light curves to which random noise of a known variance $\tilde{\sigma }_n^2$ is added. This allowed the effects of using an underestimate $\sigma _n^2 < \tilde{\sigma }_n^2$ of the noise variance in the goodness-of-fit function $G(\mathbf {\hat{J}},\mathbf {I})$ to be determined. It was found that for a given light curve, when the ratio $\sigma _n^2/\tilde{\sigma }_n^2$ falls below a certain value (typically between 0.90 and 0.98), the solution "falls apart" in that it starts to show very obvious noise artifacts. The transition to this behavior is quite sharp in that it takes place over a narrow range of values of this ratio. This gives a practical means to determine how low σ²_n can be pushed while still yielding acceptable solutions. Thus, we can avoid the oversmoothing associated with the Twomey criterion by performing inversions for a range of values of σ²_n and then choosing the lowest value which leads to a solution free of obvious noise artifacts. This value is what we call the "effective noise level" for the light curve. Because the transition is not perfectly sharp, the precise choice of σ²_n is to some extent a judgment call, so we were conservative in our choices.

The first round of inversions was thus a series of single-filter inversions of all 79 B and 79 V light curves so as to assign an effective noise to each one. The next round was to take each pair of contemporaneous B and V light curves and invert the pair in combination so as to produce the finished surface map for that pair. A complication is that when using LI to simultaneously invert light curves obtained through several filters, the effective noise levels as determined from single-filter inversions sometimes (but not always) lead to undersmoothed surface images corrupted by noise artifacts. Thus, we used the effective noise levels obtained from the initial round of single-filter inversions as starting points, and ran a series of inversions for each pair of light curves using nearby values of σ_B and σ_V. The resulting images were inspected to determine the lowest values of σ_B and σ_V that did not result in obvious noise artifacts; again we were conservative in our judgments.

In general, when simultaneously inverting multiple light curves, the deviation for a given filter between the data and reconstructed light curves differs somewhat from the effective noise. This arises because the convergence criterion for the LI algorithm is based on the overall deviation between the data and reconstructed light curves through the various filters taken together as a whole, rather than on the deviations for the individual filters. This is necessary because the scaling given by Equation (4), which defines the patch intensities assigned to the secondary filter(s) in terms of their values as seen through the primary filter, makes it impossible to independently tweak how well the individual reconstructed light curves match the data light curves.

Each of the 79 data sets consisting of paired B and V light curves was inverted according to the procedure outlined in Section 2. Upon careful examination of the resulting images, the data sets were divided into 24 groups. Each group covers a span of time during which we can see evolution of a particular set of spots, except for three cases in which a group consists of a single data set because it is temporally isolated from the data sets immediately before and after it by long gaps in the data. In six of these groups we saw good evidence for differential rotation, and as looking for such evidence was our primary goal in this work, we present in detail here only the results for these six groups, arranged in chronological order. The time spans corresponding to these six groups are indicated via horizontal bars in the plots of Figure 1. The remainder of our images are presented in Figure 18 in the online version of this paper.

Table 1 summarizes the properties of each of the data sets and groups. The data sets are numbered in chronological order in the first column, with the group to which each set was assigned indicated in the second column. Boldface entries in the first two columns denote the groups for which we present evidence of differential rotation. The starting and ending heliocentric MJDs and corresponding UTC calendar dates for the B and V light curves of each data set are given in the fourth through seventh columns; a blank entry for a V light curve indicates that the value is the same as for the B light curve from the same set. The number of data points in each data set is shown in the N_obs column. The number of rotation periods covered by each data set is listed in the "# Per." column. The last two columns list the "effective noise" in magnitudes for each inversion, as described in section Section 3, for the assumed rotation axis inclinations of i = 45° and i = 60°. The values in the table are 10⁴ times larger than the actual values, e.g., an entry of "154" means that the effective noise was 0.0154 mag.

It may be noted that the values of σ_eff tend to be large for data sets 1–10 in comparison to the rest of the data sets. The explanation for this is the replacement of the photometer with a lower-noise model between the acquisition of data sets 11 and 12, as mentioned above in connection with Figure 1. Fortunately, the amplitudes of the light curves were large during this interval, minimizing the effect of the higher noise on the quality of the light curves.

Berdyugina et al. (1998a, 1999) produced Doppler images of II Peg for 1992–98 using what they call the "Occamian" approach. We comment here on how their images compare qualitatively to our photometric images obtained at nearly the same time.

First, it should be noted that there is a consistent tendency for corresponding spots to appear at higher latitudes in their Doppler images in comparison to our photometric images. This is likely due to our having only two filter passbands available, B and V. As shown in Harmon & Crews (2000), having data through four passbands available (e.g., BVRI) significantly enhances the latitude resolution of the photometric inversions compared to when only two are used. Furthermore, a spectral line profile is effectively a one-dimensional image of the stellar surface (Vogt & Penrod 1983), while the individual data in a light curve are zero-dimensional images in that only the integrated brightness of the surface is indicated by them. Thus, Doppler imaging has intrinsically higher resolution than does photometric imaging, so it is reasonable to presume that the latitudes obtained from the Doppler images are more reliable, and we caution the reader that the latitudes of the spots on our photometric images should not be accepted at face value. However, the extensive simulations reported in Harmon & Crews (2000) show that relative spot latitudes can be reliably determined from LI: when two spots at different latitudes are present, the lower-latitude spot is generally rendered as such, even though the individual spot latitudes may not be accurately rendered, particularly when only two-filter data are available. Thus, we should be able to detect differential rotation in our photometric images, though we cannot go as far as to reliably determine the rotational period as a function of latitude on the stellar surface.

The 1992 August image of Berdyugina et al. (1998a) compares favorably to our Data Set 12 (1992 September 7–October 3) photometric image, though the spot latitudes on the Doppler images are considerably higher. In both images, a large spot has an appendage projecting toward lower longitudes. Their 1993 December image shows two distinct high-latitude spots almost 180° apart, with one moderately larger than the other. Our Set 18 (1993 December 24–1994 February 1) image shows what appear to be two widely separated spots with a "bridge" connecting them. The bridge may be simply an artifact of the inversion, as such bridges are seen in simulations in which two spots close to one another are present on the artificial stellar surface. Their 1994 November and 1995 January images show a pair of spots at higher latitudes than our Group 8 (1994 October 13–1995 February 9) images. Our images show a pair of spots connected by a bridge. Their 1995 July image and our Set 22 (1995 June 5–July 11) image both show a high-latitude feature which extends across a wide range in longitude. Their 1995 October image and our Group 10 (1995 September 22–1996 February 16) images show a pair of isolated spots well-separated in longitude, with one larger than the other.

Their Doppler image for 1996 October shows II Peg shortly before one of the time intervals over which we see good evidence for differential rotation, that of Group 11 (1996 November 4–1998 February 12), and their images for 1997 June, August, and December show the star during the interval covered by our images. Our results for this interval are described in Section 4.3 below. Similarly, their Doppler images for 1998 October and November represent the star during another interval over which we saw evidence of differential rotation, that of Group 12 (1998 September 29–1999 February 15), which we discuss in Section 4.4. We defer the comparison of their Doppler and our photometric images for these intervals to those sections.

4.1. Group 2: MJD 47417.7681–47556.5813

Figure 2 shows plots of the light curves for the four data sets comprising Group 2, which span the interval from heliocentric MJD 47417.7681–47556.5813 (1988 September 11–1989 January 30). Also shown in Figure 2 are the reconstructed light curves, i.e., the light curves computed from the stellar surface maps, for the cases in which the assumed inclination of the rotation axis to the line of sight is i = 45°. The intensities have been normalized such that the brightest datum over the entire data set considered in this study for the given filter is set equal to 1. The quantity σ on each plot represents the rms deviation expressed in magnitudes between the data light curve and the reconstructed light curve.

Zoom In Zoom Out Reset image size

Figure 3 shows the reconstructed surfaces obtained from these light curves via the LI procedure, with the results for an assumed axial inclination of i = 45° shown in the top row and for i = 60° shown in the bottom row. Each column of two images consists of the pair of i = 45° and i = 60° reconstructions for the date ranges indicated at the top of the column. Note that if we let ϕ → −ϕ in the specific intensity distribution J(θ, ϕ), where θ and ϕ are stellar latitude and longitude, and also reverse the direction of stellar rotation, the resulting rotational light curve is unchanged. Thus, based on photometry alone we cannot distinguish between a given brightness distribution and its mirror image. For definiteness, in generating our images we have adopted the convention that the star rotates counterclockwise as viewed from above the visible pole, so that spots would be carried across our view from left to right with the visible pole at the top of the image.

Zoom In Zoom Out Reset image size

We do not show here the other light curves which we inverted. Figure 2 is included so that the reader may get a sense of the comparison between the data and reconstructed light curves for a typical example. However, information about the characteristics of the source light curves along with the displayed longitude of disk center and assumed axial inclination is provided in the form of seven numbers which appear below each of the reconstructed surface images we present. From left to right, these are as follows, using the values for the image at upper left in Figure 3 as examples.

• 1.  
  "190" indicates that the longitude of disk center for the surface as displayed is ϕ = 190°. This is not to be confused with the actual sub-Earth longitude on the stellar surface, which varies from observation to observation. Unless otherwise indicated, the latitude of disk center for all surface images is θ = +30°, as opposed to the actual sub-Earth latitude of 90° − i.
• 2.  
  "45" indicates that i = 45° for the reconstruction.
• 3.  
  "20" and "23" are the number of data points in the B and V light curves used for the reconstruction, respectively.
• 4.  
  "326" is the "effective rms noise" expressed in magnitudes (see Section 3) for the B-filter light curve, multiplied by 10⁴, so that the actual effective rms noise is 0.0326 mag. Similarly, "268" indicates that the effective rms noise for the V-filter light curve is 0.0268 mag. Since the V-filter data were of lower noise in this case, V was designated the "primary filter" n = 1 for this inversion in the sense of the discussion in Section 2. Note that the effective rms noise used in the inversion differs from actual deviation σ = 0.0311 mag reported in Figure 2 for the B light curve of Set 3. The reason for this is discussed at the end of Section 3.
• 5.  
  "2.7" indicates that light curve data used to create the images span 2.7 revolution periods of the binary, which we assume also to be the rotation period of the star because of tidal locking. (Note that if there is differential rotation, this rotation period can, strictly speaking, only apply to some particular latitude.)

Before discussing the images, it should be noted that simulations (Harmon & Crews 2000) show that when two separate circular spots are close together, they may appear in the reconstructions as a single elongated spot because of limitations in the resolution of the method. Indeed, what appears to be a single large spot may in fact be a magnetic active region having complex structure which cannot be resolved. We will thus use "spot" and "active region" somewhat interchangeably, but tending to prefer the latter terminology for a more extensive dark region on an image with a decidedly non-circular appearance. Furthermore, as mentioned above, "bridges" which appear to connect spots together may be artifacts of the inversion procedure.

The Set 3 images for both i = 45° and i = 60° appear to show a lower-latitude spot at left and a higher-latitude spot at right connected by a bridge. The persistence of spottedness at high latitudes near the upper center of each image suggests that there is indeed a high-latitude magnetic active region in this location, and that the bridge is not simply an inversion artifact in this case.

If II Peg exhibits differential rotation in the same sense as the Sun in that lower latitudes have higher angular velocities, then we should see the longitude difference between the lower-latitude spot and the other features diminish over time. This does appear to happen as one peruses the images from left to right. In particular, the lower-latitude spot appears to catch up to and move underneath the high-latitude active region. This could represent two separate spots which appear merged in the images due to the limited resolution of the reconstructions, or it could represent a single large active region extended in latitude, which rotates (roughly speaking) about an axis through its center due to shearing by differential rotation carrying its southern end around the star faster than its northern end.

The interpretation of the spot appearing at far right in the Set 3 images is problematic because it clearly appears at a higher latitude than the spot on the left in this pair of images, but in the Set 4 images the spot on the right appears to be at a comparable latitude to the one on the left. This may simply represent a limitation in the latitude discrimination of the spots in the inversions, which certainly is not perfect, particularly given that only B and V light curves were inverted as opposed to having more filter data available. The spot on the right appears in the i = 45° image for Set 5, but not in the i = 60° image. In simulations it is seen that a small spot close to a much larger spot may be "masked" by the larger spot, causing the smaller spot not to appear in the reconstruction. That this spot appears in the i = 60° image suggests that such masking has occurred for the i = 45° image.

It is important to note that Sets 3 and 4 were acquired over a total span of 7.9 stellar rotations and the corresponding images differ noticeably, while Set 5 by itself spans 8.2 rotations. It is thus likely that there was significant evolution of the spots during the time Set 5 was acquired. In general, all the surface images shown here should be regarded as representing a sort of "average" appearance of the starspots during the interval over which the data were collected. It should also be noted that only 12 B and 13 V observations were available for the Set 6 time span. Thus, some additional caution should be used in interpreting the corresponding images.

Figure 4 shows plots of the total span in longitude Δϕ of the spotted region versus the heliocentric MJD of the midpoint of the time interval over which the data were collected corresponding to each image shown in Figure 3. The plot on the left is for i = 45°, while the one on the right is for i = 60°. To obtain Δϕ, an IDL⁵ widget which displays the reconstructed image for user-adjusted values of the latitude and longitude of disk center was employed. The stellar equator was placed at disk center, and the longitude of disk center varied so as to determine the longitudes of the easternmost and westernmost dark patches within the spotted region. This is admittedly an imprecise determination, since there is no sharp delineation between what constitutes a "spot patch" versus a "photosphere patch." Nevertheless, in practice the edges of the spotted regions in the images are fairly well defined, as can be seen in Figure 3 and all the other surface reconstructions shown in this work. A line was then fitted to the data using a simple unweighted least-squares procedure. The slope of the line represents a measure of the rate at which the longitude span of the spotted region is changing, expressed in degrees per day. There is considerable variation between the results for i = 45° and i = 60°, but both plots behave as would be expected if II Peg exhibits differential rotation in the same sense as the Sun: Δϕ decreases with time, because the trailing region of the active region is at lower latitudes than is the leading region. In Section 5, we show that the fits shown in Figure 4 and similar figures for the other data sets showing evidence of differential rotation are consistent with the conclusions of Henry et al. (1995) regarding the value of the differential rotation coefficient k defined in Equation (1).

Zoom In Zoom Out Reset image size

Because a spot of a given size and location and a lower temperature will modulate the light curve to a greater degree than will a spot of higher temperature, it is conceivable that reducing our assumed value of T_spot would result in smaller spots on the reconstructed surface. This in turn would likely reduce the longitude extent we obtained for the spotted regions and thus alter plots like those in Figure 4. If so, it would be necessary for us to invert for a grid of spot temperatures and compare the results. To examine the effects of variations of T_spot on our models, we inverted all the light curves in Data Group 2 using T_spot = 3250 K and 3750 K in addition to our "canonical" value of 3500 K, a range that more than covers the estimates by other authors of T_spot ≈ 3300–3700 K reported in Section 1. Figure 5 shows the resulting surface images for assumed inclination i = 60°, with the top row corresponding to T_spot = 3250 K, the middle to T_spot = 3500 K, and the bottom to T_spot = 3750 K.

Zoom In Zoom Out Reset image size

Inspection of Figure 5 shows that the geometrical properties of the spots are nearly identical, and indeed with few exceptions the same sets of surface patches are notably darkened in comparison to the surrounding "photosphere" on corresponding reconstructions, as was verified via detailed examination using the IDL display widget. The most significant difference between the images is that the "core" regions of spots containing the darkest patches tend to be smaller when T_spot is smaller and larger when T_spot is larger; the boundaries of the spots are for the most part unchanged. Thus, varying T_spot over the range considered here does not significantly affect the spot shapes, sizes, or locations. The same conclusion was reached upon examining the inversions for i = 45°.

These results gave us confidence that it was not necessary to use a grid of spot temperatures, but as a further check we also inverted the first data sets in all the other data groups for which we saw differential rotation (data sets 7, 27, 34, 38, and 42). In all cases, careful inspection showed that in longitude, spot boundaries varied by, at most, the width of one patch on the surface, though in most cases not even by that. We therefore decided that it was sufficient to invert all the other data sets only for T_spot = 3500 K.

4.2. Group 3: MJD 47779.7480–47926.5752

Figure 6 shows inversions of Group 3, consisting of three data sets obtained from MJD 47779.7480–47926.5752 (1989 September 10–1990 February 4). There are no light curve data for the interval 1989 January 31–September 9. Thus, no information is available regarding the evolution of the surface between the end of the time span covered in Section 4.1 and that covered here.

Zoom In Zoom Out Reset image size

The images in Figure 6 give the distinct impression that the low-latitude spot on the right is moving away from higher-latitude active regions to its left due to differential rotation, though it should be said that the spot at far left in the third column appears to be a new feature which has emerged within the interval covered by these light curves.

An interesting question is whether or not the low-latitude spot at right in the images in Figure 6 is the same as the low-latitude spot at left in the images in Figure 3. If so, then over the combined span MJD 47417.7681–47926.5752 (1988 September 11–1990 February 4) we see the low-latitude spot go from being well behind the high-latitude spot/active region in longitude to being well ahead of it, and thus a distinctive manifestation of differential rotation. The images are consistent with this hypothesis. The midpoint of data set 7 corresponding to the leftmost images in Figure 6 was at MJD 47798.8, while the midpoint of data set 6 corresponding to the rightmost images in Figure 3 was at MJD 47451.1, so that 257.7 days elapsed between the two midpoints. Taking |d(Δϕ)/dt| = 0.345 deg day⁻¹ for i = 45° from Figure 4 as a rough estimate of rate at which the low-latitude spot should separate in longitude from the high-latitude spot/active region near the central meridian of each figure, we would expect the longitude of the low-latitude spot to have increased relative to the high-latitude spot by approximately 90° in this time span. Using |d(Δϕ)/dt| = 0.185 deg day⁻¹ for i = 60° yields a relative increase in longitude of 48° by the low-longitude spot. The actual amount by which the low-latitude spot is ahead of the high-latitude spot in the images for Set 7 is about 60°. If the images for Set 6 in Figure 3 are taken to represent the two spots having approximately equal longitudes, then the low-latitude spot being 60° ahead in Set 7 is within the range of 48°–90° suggested by the plots in Figure 4.

Figure 7 shows how Δϕ for the Group 3 data sets varies in time. In this case we expect d(Δϕ)/dt>0 for solar-sense differential rotation, since the low-latitude spot now starts out ahead in longitude. The magnitudes of the slopes are comparable to those in the plots in Figure 4, strengthening the interpretation that we are seeing the effects of differential rotation throughout this data set and that discussed in the preceding section.

Zoom In Zoom Out Reset image size

4.3. Group 11: MJD 50391.8348–50856.5946

Figure 8 shows the inversions of Group 11, consisting of seven data sets obtained from MJD 50391.8348–50856.5946 (1996 November 4–1998 February 12), all for an assumed inclination of i = 45°. Figure 9 shows inversions of the same light curves for i = 60°. The first four images (Sets 27–30) in both figures appear to show a trailing low-latitude spot at left catching up with and passing under a high-latitude feature near the central meridian, highly suggestive of differential rotation in the same sense as that of the Sun. The last three images do not continue this trend, but also cover a relatively short time interval.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The Set 33 images give the appearance of the sense of the differential rotation having reversed, but it should be noted that both the B and V light curves inverted to create these images had only 13 data points, which calls the reliability of these images into question. It is also possible that new high-latitude spot activity emerged ahead of the older activity.

Berdyugina et al. (1998a) present Doppler images for 1996 October and 1997 June, August, and December, the first of which shortly precedes and the latter three of which overlap in time with the photometric images discussed in this section. We have no data for 1996 October, but their image gives some indication of the high-latitude region seen in between the two initial spots in our later images. As is typically the case, the spots seen in their Doppler images are at higher latitudes than in our photometric images.

Their 1997 June image shows a large high-latitude active region with a large (∼150°) extent in longitude, along with a pair of lower-latitude spots. In their 1997 August image, the extent in longitude of the large active region has diminished considerably, and it appears that this is largely due to the disappearance of activity at the leading edge of the region. A smaller spot persists on the opposite side of the star. In their 1997 December image, the width in longitude of the large active region has further diminished, but this time it appears that either the activity along the trailing edge has diminished, or that the trailing edge has caught up to other parts of the region. There is a southward projection at the trailing edge which is not present in the 1997 August image.

Our Set 29 (1997 May 22–July 13) images differ considerably from their 1997 June image. Ours show a spot close to the equator trailing a higher-latitude active region. Our Set 27 (1996 November 4–December 18) images show the larger low-latitude spot present in their 1997 June image, which appears to fade in our Set 28 (1996 December 19–1997 February 15) images. Our Set 30 (1997 September 23–November 3), Set 31 (1997 November 4–December 13), and Set 32 (1997 December 15–1998 January 24) images are comparable to their 1997 December image in showing a single large active region on one side of the star. The southern protrusion in their image could explain the fact that our images show activity down to near the equator, though in their image the protrusion only extends down to ∼35° latitude.

Figure 10 shows plots of the variation in Δϕ with time for the seven data sets in Group 11. The two plots at top include only the results for the first four data sets for each inclination, while the two at bottom include all seven sets. The fits are better for the former two since the trend toward narrowing of the spotted region in longitude ceases for the latter three data sets.

Zoom In Zoom Out Reset image size

4.4. Group 12: MJD 51085.9791–51224.5942

Figure 11 shows inversions for i = 45° of Group 12, consisting of four data sets obtained from MJD 51085.9791– 51224.5942 (1998 September 29–1999 February 15), while Figure 12 shows the same for i = 60°. In both figures, the images in the top row have latitude θ = 30° at disk center, while for the images in the bottom row the north (visible) rotation pole is at disk center.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

For both assumed inclinations, the images show a leading low-latitude spot or active region at right and a trailing high-latitude spot/active region at left. As time progresses, the separation in longitude between these features increases, just as would be expected if solar-sense differential rotation is carrying the low-latitude feature around the star faster than the high-latitude region. This is most easily seen in the images looking down on the pole. Note that the degree of magnetic activity appears to increase on both ends of the spotted region during this time span.

Berdyugina et al. (1999) present Doppler images for 1998 October and November which overlap the interval under consideration here. Their 1998 October image shows an active region spread extensively in longitude that is not resolved into separate spots; in the November image, this region appears split into three separate spots. The difference is conceivably due to the variation in the noise artifacts between the two images. Our data set 35 (1998 October 30–November 27) image shows a similar elongated feature, but lacks the southward protrusion seen in the middle region of the active region in the Doppler images.

Figure 13 shows the variation of Δϕ with time for this time interval. As expected, Δϕ increases with time, since the low-latitude, faster-moving activity is leading in longitude.

Zoom In Zoom Out Reset image size

4.5. Group 13: MJD 51429.8915–51586.5994

Figure 14 shows inversions of Group 13, composed of four light curves obtained between MJD 51429.8915–51586.5994 (1999 September 8–2000 February 12). For both assumed inclinations, the trailing edge of the spotted region appears to be at lower latitudes than the leading edge. The last pair of images appear to show that the low-latitude trailing spot (or extension of the spotted region) has caught and passed underneath higher-latitude activity.

Zoom In Zoom Out Reset image size

Gu et al. (2003) produced surface maps of II Peg via Doppler imaging for 1999 July–August, 2000 February, and 2001 November–December, the first two of which are relevant to the present discussion. Their 1999 July–August image is qualitatively similar to our Set 38 (1999 September 8–October 24) images, showing a broad region of high-latitude activity with southward projections at its leading and trailing ends, with the trailing end projecting further south. Their 2000 February image primarily gives the appearance that the active region may have simply shrunk in longitude. However, it could also be argued that the trailing edge has caught up to the leading edge.

Figure 15 shows that the extent in longitude Δϕ of the spotted region described here diminishes with time, as expected for solar-sense differential rotation.

Zoom In Zoom Out Reset image size

4.6. Group 14: MJD 51805.8758–51946.6051

The final time interval over which we saw good evidence for differential rotation was from MJD 51805.8758–51946.6051 (2000 September 18–2001 February 6), which we partitioned into three data sets which were assigned to Group 14. Figure 16 shows the corresponding inversions.

Zoom In Zoom Out Reset image size

For both values of i, a protrusion toward lower latitudes can be seen at the bottom of the spotted region. This protrusion appears to drift toward greater longitudes over the time interval covered by these observations, suggesting that it is a spot being carried around the equator faster than the rest of the spotted region because of differential rotation.

However, we cannot quantify this differential rotation by plotting the variation of Δϕ with time as for the data groups discussed above, because in this case the total extent in longitude of the active region did not change significantly. The easternmost and westernmost fringes of the active region appear to lie at similar middle latitudes, so that we would not expect Δϕ to change due to differential rotation here. The longitude of the southward protrusion is within the range of longitudes spanned by the activity to its north, so that its drift in longitude relative to the rest of the spotted region is not manifested in the variation of Δϕ.

Instead, in Figure 17 we show the variation in longitude of the middle of the southern protrusion with time. This is admittedly a vague conception due to the amorphous nature of the protrusion in these images, which is likely just an artifact of the limitations in the resolution of the inversions. The procedure was to use the same IDL widget used to obtain Δϕ for other data sets to estimate the longitude of the midpoint of the southern boundary of the protrusion. The same results were obtained for both i = 45° and i = 60°. As can be seen in Figure 17, the longitude of the protrusion (or spot) increases approximately linearly with time, as would be expected from differential rotation if the protrusion or spot maintains an essentially fixed latitude.

Zoom In Zoom Out Reset image size

The intervals in which we saw evidence for differential rotation are not randomly distributed. Data Groups 2 and 3 span 1988 September to 1990 February, then Data Groups 11–14 span 1996 November to 2001 February, which is suggestive of the 6.8 yr flip-flop period obtained by Rodonò et al. (2000). Indeed, examination of Figure 7 of that paper shows that our Data Groups 2 and 3 overlap a time interval in which their longitude C is ceasing activity while longitude B is resuming activity, these two longitudes being the ones involved in the flip-flop behavior. Similarly, our Data Group 11 closely corresponds to the subsequent end of the activity of their longitude B and resumption of activity of longitude C. Perhaps a common physical mechanism underlies the flip flops and an enhanced ability to detect differential rotation, such as an increase in the spread in latitude of the active regions present on the surface. However, it should be noted that we continued to see differential rotation for Data Groups 12–14, when their longitude B was inactive.

Extrapolating the results of Rodonò et al. (2000) implies that additional flip flops should have occurred in the vicinity of MJD 52500 and 55000, but we saw no evidence of differential rotation at these times. Inspection of Figure 1 shows that the amplitude of the modulation was small then, while it was relatively large during the intervals for which we saw differential rotation. It may well be that it is the large areal coverage by the spots implied by the large amplitude which was primarily responsible for our ability to see differential rotation for these intervals.

In contrast with the relation used by Henry et al. (1995) and given by Equation (2), in the present study it is more natural to characterize the differential rotation in terms of the rate of change in the difference in longitude between spots at different latitudes in our images. If for simplicity we express longitudes in a non-rotating frame of reference (rather than in the co-rotating frame used in the rest of this work), the longitude as a function of time t for a spot having latitude θ is just

$$\begin{equation} \phi (t) = \phi _0 + t\Omega (\theta), \end{equation} \tag{ 17 }$$

where ϕ₀ is the longitude of the spot at time t = 0, and Ω(θ) = 2π/P(θ) is the angular rotation frequency for latitude θ, with P(θ) being the corresponding rotational period. From this the difference in longitude Δϕ = ϕ₂(t) − ϕ₁(t) of two spots at latitudes θ₁ and θ₂ is simply

$$\begin{equation} \Delta \phi = \Delta \phi _0 + t\Delta \Omega, \end{equation} \tag{ 18 }$$

where Δϕ₀ is the separation at t = 0 and ΔΩ = Ω(θ₂) − Ω(θ₁). Thus, a plot of Δϕ versus t will have slope d(Δϕ)/dt = ΔΩ. In terms of the rotational angular frequency, Equation (1) becomes

$$\begin{equation} \Omega (\theta) = \Omega _\mathrm{eq}(1 - k\sin ^2\theta), \end{equation} \tag{ 19 }$$

so that

$$\begin{equation} k = \frac{\Delta \Omega }{\Omega _\mathrm{eq}(\sin ^2\theta _1 - \sin ^2\theta _2)}. \end{equation} \tag{ 20 }$$

In principle, then, if we observe two spots at known latitudes and plot Δϕ versus t to find the slope ΔΩ, we can find k from Equation (20). In practice, this is not possible because we cannot reliably determine precise spot latitudes from our photometric inversions.

However, we can use the slopes d(Δϕ)/dt of the plots in Figures 4, 7, 10, 13, and 15 to estimate the value of k by making reasonable assumptions about the range in spot latitudes. Using R = 3.4 R_☉ and vsin i = 22.6 km s⁻¹ from Berdyugina et al. (1998b), we find Ω_eq = 0.9529day⁻¹. Table 2 shows the values of k obtained using this estimate for Ω_eq, with the assumed minimum and maximum spot latitudes indicated by the notation θ = (θ_min, θ_max).

Table 2. Differential Rotation Coefficient, k

	ΔΩ (deg day−1)		θ = (0°, 90°)		θ = (45°, 80°)		θ = (482, 717)a
Group	i = 45°	i = 60°	i = 45°	i = 60°	i = 45°	i = 60°	i = 45°	i = 60°
2	0.34	0.19	0.0063	0.0034	0.013	0.0072	0.018	0.0098
3	0.28	0.15	0.0051	0.0027	0.011	0.0058	0.015	0.0079
11b	0.24	0.21	0.0044	0.0039	0.009	0.0084	0.013	0.011
12	0.36	0.31	0.0066	0.0057	0.014	0.012	0.019	0.017
13	0.39	0.29	0.0072	0.0052	0.015	0.011	0.021	0.015
14	0.42	0.42	0.0076	0.0076	0.016	0.016	0.022	0.022
		Mean	0.0062	0.0048	0.013	0.010	0.018	0.014
		Standard deviation	0.0012	0.0018	0.003	0.004	0.004	0.005

Notes. ^aValues taken from Siwak et al. (2010). ^bData sets 31–33 excluded.

Download table as:  ASCIITypeset image

Henry et al. assumed θ_min = 0° and θ_max = 90° to obtain k = 0.005 ± 0.001. Table 2 shows that the results of the present study accord well with those of Henry et al. when the same range of spot latitudes is assumed, which lends credence to our assertion that our images represent genuine demonstrations of differential rotation, and to the assertion of Henry et al. that stretching of active regions by differential rotation accounts for features of the light curve of II Peg not explainable by their simple two-spot model.

Of course, the values obtained for k are increased if the span of latitude is not from equator to pole. The Doppler images of Berdyugina et al. (1998a, 1999) show activity extending from mid-to-high latitudes, but not including any polar spots, so it seems reasonable to use θ₁ = 80° and θ₂ = 45° in our estimates. Table 2 shows that these results do not agree as well with the result of Henry et al., but it should be noted that Henry et al. would also have obtained a larger value of k had they assumed a smaller range of latitudes.

Recently, Siwak et al. (2010) analyzed a very high quality light curve of II Peg obtained by the MOST satellite during 31 consecutive days spanning HJD 2454725–2454756 (2008 September 15–October 16), covering more than 4.5 rotations. They fit the light curve with a three-spot model using the starspotz analytic spot modeling program (Croll et al. 2006), all the spots being assumed circular, with photosphere and spot temperatures of T_phot = 4600 K and T_spot = 3600 K. The input parameters for the spots were the latitude, angular diameter, rotation period, and the initial time at which the spot most directly faced the observer. One very large spot centered on the hidden pole (and thus only partly visible from Earth) was used to prevent the radii of the other two spots from being so large that they overlapped, which is not admissible in starspotz. Models were created assuming rigid and differential rotation, and for which proximity effects due to ellipsoidal distortion and reflection arising from the presence of the secondary star were and were not taken into account. Whether or not proximity effects were included, Siwak et al. found that models including differential rotation fit the data better than models assuming rigid rotation, in particular in their ability to account for a steady reduction in the amplitude of the modulation. The fits exhibit reduced weighted χ² values much larger than unity, which Siwak et al. ascribe to the spot shapes not being circular as assumed by starspotz. The observations of Siwak et al. overlap our data sets 73 and 74 of Group 23, for which our inversions are shown in Figure 18 (see images 20 and 21) in the online version of this paper. Careful comparison of the T3 Automated Photometric Telescope observations used in this work and the MOST data shows that the data sets are consistent with one another, but the time resolution of the APT data was not sufficient to reveal the slow decay in the amplitude of the modulation seen by MOST because the peaks were not well sampled. This coupled with the small mean amplitude is likely why we did not see evidence for differential rotation for these sets.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

View figure set (23 panels)

Tarball of figure set images

When proximity effects were included, Siwak et al. obtained k = 0.0245^+0.0155_−0.0020 for the differential rotation coefficient, almost five times greater than that of Henry et al. and considerably larger than all those discussed so far in relation to this study. The nominal value of k corresponded to spots in their model surface having latitudes of 482 and 717. The last two columns in Table 2 show that assuming this range of latitudes in our results leads to values of k more nearly in accord with those of Siwak et al., though ours are still smaller than theirs.

N.V. acknowledges support from the Ohio Wesleyan University Summer Science Research Program. G.W.H. acknowledges support from NASA, NSF, Tennessee State University, and the State of Tennessee through its Centers of Excellence program.

The IDL widget used to interactively display the images was created by Ryan Deskins.

The authors thank the anonymous referee for helpful comments and suggestions for the improvement of this manuscript.