Probing the rotation rate of solar active regions: the comparison of methods


 Sunspot groups are often used as tracers to probe the differential rotation of the Sun. However, the results on the rotation rate variation obtained by different authors are not always in agreement. The reason for this might be a number of effects. In particular, faster decay of the following part of a sunspot group results in a false apparent shift of the area-weighted center of the group toward the leading part. In this work we analyze how significantly this effect may contribute to the derived rotation rate. For a set of 670 active regions, we compare the rotation rate derived from continuum intensity images to that derived from line-of-sight magnetograms. We found that the difference between the calculated rotation rates is 0.45° day−1 on average. This value is comparable to the difference between the rotation rate of the solar surface near the equator and at 30° latitude. We conclude that the accuracy of the rotation rate measurements using white-light images is not satisfactory. Magnetograms should be used in future research on the differential rotation of the Sun.


Introduction
The differential rotation of the Sun is an important mechanism involved in the generation of solar magnetic fields -the global dynamo process. The differential rotation presumably converts the poloidal global magnetic field to the toroidal one (see, e.g., a review by Charbonneau 2020). The differential rotation was observed for the first time by tracking the sunspots traversing the solar disk back in 1630s (Beck 2000). Recent achievements in the observational instrumentation and mathematical approaches offer sophisticated tools to measure the rotation rate of the solar convective layers at different depths. The description of the methods used to derive the differential rotation of the Sun can be found in an excellent review by Beck (2000).
Sunspots are widely used as tracers to probe the solar differential rotation. The rotation rate is usually estimated by measuring the consecutive positions of an individual sunspot or of a sunspot group observed at different times. Long-term sunspot observations covering four centuries allow one to analyze the rotation rate of sunspots in detail.
Many analyses suggest that the rotation rate of sunspot groups depends on size and age of the tracer. For instance, Beck (2000) reviewed a number of previous results and concluded that larger sunspot groups tend to rotate slower as compared to smaller ones. At the same time, young sunspot groups tend to rotate faster in comparison with older sunspot groups. Faster rotation of young sunspot groups as compared to older ones implies that sunspot groups decelerate with time. This conclusion was supported by Pulkkinen and Tuominen (1998) who performed the most solid measurements of the rotation rate for sunspot groups so far. The authors utilized almost one and a half century long data series of sunspots. Some of their results are presented in Table 1 in the next sections.
An opposite deduction was made in a number of works by Javaraiah and Gokhale (1997);Hiremath (2002); Sivaraman et al. (2003). The authors found that young sunspot groups tend to accelerate. This finding is not in agreement with the results mentioned above. Although the sunspot databases used in most works were the same, there were differences in the approaches applied to analyze the data, for example, in sorting the sunspot groups by age or lifetime.
In addition to different data handling methods, a number of mechanisms (e.g. Petrovay and Christensen 2010) and artifacts may affect the results of the rotation rate measurements. Petrovay (1993) argued that faster decay of the following part of a sunspot group yields apparent motion of the area-weighted center toward the leading part (see Figure 2 in Petrovay 1993). The asymmetry of the emerging magnetic flux loop forming an active region may result in different proper motion of the following and leading part of a sunspot group (e.g. van Driel-Gesztelyi and Petrovay 1990).
The aim of this work is to analyze the influence of the choice of the data on the measured rotation rate of solar active regions. In other words, how dramatically the results would change if we use magnetograms instead of whitelight images. We do this by comparing the rotation rates derived from continuum intensity images and from magnetograms of the same set of active regions. Section 2 describes the data reduction and the methods applied to calculate the rotation rates. The results of the calculations as well as the comparison of the derived values are presented in Section 3. Finally, Section 4 draws the conclusions of this study and discusses the constraints for future research.

Data reduction
In this work the data acquired by the Helioseismic and Magnetic Imager (HMI, Schou et al. 2012) on board the Solar Dynamics Observatory (SDO, Pesnell et al. 2012) were utilized. SDO/HMI is a full-disk filtergraph aiming at measuring solar magnetic fields, Doppler velocities, and other quantities. The instrument uses the Fe I 6173 Å photospheric line. The resolution of SDO/HMI is 1 ′′ with the pixel size of 0.5 ′′ ×0.5 ′′ . The data on magnetic field and continuum intensity are available at 45 s and 720 s cadence.
The data reduction technique in this work is similar to that applied in Kutsenko (2021). Briefly, each active region was manually bounded by a rectangular box at the fulldisk line-of-sight SDO/HMI magnetogram. Then, the active region was tracked at the solar surface for several days as the Sun rotated. The bounded part of the magnetogram was cropped and stored as a magnetic field patch of the active region. The same part of the co-temporary continuum intensity images was also cropped and stored. The tracking stopped when the bounding box center was farther than 60 ∘ from the central meridian. As a result, we obtained two cubes of data containing the maps of magnetic field and continuum intensity of selected active regions. The cadence of the data was set to two hours.
To calculate the heliographic coordinates of the active region at the solar disk we applied the procedure described in detail in Kutsenko (2021). The size of magnetogram and continuum intensity patches was reduced by binning 2×2 pixels. For magnetograms, the magnetic flux centroids of positive and negative magnetic polarities were calculated. The center of the active regions was derived as a geometrical center of these magnetic flux centroids. Note that for unipolar active regions the center of the active region was set as the flux-weighted center of the dominant magnetic polarity. For continuum intensity images, we measured the position of the intensity-weighted center of the active region. The centers of the active region in the pixel coordinates were converted to the heliographic coordinates using IDL World Coordinate System library. As a consequence, at this step we derived two different longitudinal positions of an active region for the entire interval of observations.
The difference between the derived active region centers can be seen in Figure 1. The top panels of the figure show the continuum intensity images of NOAA active region 12489 observed between 2016 January 27 and 2016 February 01. The red crosses mark the derived intensity-weighted center. The co-temporary line-of-sight magnetograms of the active region are shown in the bottom panels. In this case the red crosses mark the geometrical center of the opposite magnetic polarity flux-weighted centroids. One can see an obvious proper motion of the intensity-weighted center toward the leading spot at the top panels. At the same time, the geometrical center at the bottom panels remains nearly  at the same position with respect to the leading and following parts of the active region. Consequently, the geometrical center does not migrate toward the leading part. For each patch of each active region we calculated its total unsigned magnetic flux by summing the pixel magnetic flux density multiplied by the pixel area. The summation was performed over the pixels with the absolute magnetic flux density exceeding 18 Mx cm −2 . Prior to summation the µ-correction (Leka et al. 2017) was applied in order to decrease the influence of the projection effect.
The synodical rotation rate of an active region was evaluated by approximating the longitude versus time curve by a linear fitting. The total unsigned magnetic flux curve was used to manually select the fitting interval. The interval covered the mature state of the active region. The unstable emergence phase as well as late decay phase were excluded from the fitting (see Figure 1 in Kutsenko 2021). Note that an implicit sorting of active regions occurs in some sense: the rotation rate is derived during a certain phase of the active region evolution.
The sidereal rotation rate was calculated from the synodical one by adding a term accounting for the rotation of the Earth with respect to the Sun. The details of the calculation of this term can be found in Skokić et al. (2014) and in Lamb (2017). Through the rest of this paper, the sidereal rotation rates derived from magnetograms and from continuum intensity images will be denoted ω m sid and ω ci sid , respectively.

Results
In all, in this work we analyzed 670 active regions observed between 2010 and 2016. The peak total unsigned magnetic flux of the active regions varied in the range 1.0-12.4×10 22 Mx. Among 670 active regions, 174 ones were identified as unipolar active regions while 496 were bi-or multipolar. The active regions were observed up to heliographic latitudes ϕ of about ±35 ∘ .
The distribution of the difference between the rotation rates ω ci sid − ω m sid is shown in Figure 2. The red curve in the plot shows the histogram for unipolar active regions. As expected, the histogram is centered near the zero value. The blue curve in Figure 2 is related to bi-and multipolar active regions. The distribution is clearly shifted towards positive values. This implies that the continuum intensity images yield in most cases higher rotation rates as compared to the magnetograms. Obviously, the proper motion of the intensity-weighted center toward the leading part in continuum intensity images causes the false acceleration of an active region. The mean value of the ω ci sid − ω m sid difference distribution for bi/multipolar active regions is 0.45 ∘ day −1 . We consider only bi-and multipolar active regions through the rest of this paper.
In order to check whether the tracer rotation rates derived from magnetic field and continuum intensity data are statistically different we performed a simple z-testing. The mean values of the tracer rotation rates are 14.261 ∘ day −1 and 14.706 ∘ day −1 for magnetograms and continuum intensity, respectively. The estimated standard deviations of the distributions are 0.293 ∘ day −1 (magnetograms) and 0.708 ∘ day −1 (continuum intensity). The null hypothesis assumes that the distributions have the same mean value, i.e. the difference between the calculated mean values is due to errors of the measurements. Trivial calculation yields the z-value of 12.9 (the number of data points is N = 496) that is more than the critical z-value ≈ 3.3 for 99.9% confidence interval. Hence, we must reject the null hypothesis and conclude that the tracer rotation rates derived from magnetograms and continuum intensity images form statistically different distributions. The relationship between the rotation rate and the heliographic latitude ϕ for bi/multipolar active region was fitted by a commonly used equation (e.g. Beck 2000) where constant A is the equatorial rotation rate, constants B and C are responsible for the latitudinal dependence of the differential rotation. The derived constants for ω ci sid and ω m sid are listed in Table 1 in rows 5 and 8, respectively. The A,B, and C values derived in a number of other works are presented for comparison. One can see that the constants depend significantly on the data used to probe the rotation rate. Thus, the equatorial rotation rate derived from contin- uum intensity images exhibit systematically higher values compared to that derived from magnetograms. We further divided all the tracers into 2-degree-wide latitudinal bins. For each bin we calculated the mean value and the standard deviation of the tracer rotation rate. The results are presented in Figure 3. Circles represent mean values of the rotation rate of the tracers in each bin. Solid lines in the figure show the best fitting of ω ci sid and ω m sid by equation 1. Note that mean ω ci sid exhibit systematically higher value within the entire latitudinal range. Recall that the same tracers were used to derive the rotation rates shown in blue and red in Figure 3. The difference between the measured rotation rates is exclusively due to the choice of the initial data -either magnetic field maps or continuum intensity images.
We also studied the latitudinal distribution of the tracer rotation rates. To do this in an appropriate way, we could analyze the distribution of the rotation rates of tracers located within 2-degree-wide or, say, 5-degree-wide latitudinal bins. However, we found that in this case the number of points in the bins is too small for some reasonable analysis. The analysis of the entire distribution of the tracer rotation rates without binning might also be incorrect: the expected value of the rotation rates decreases as the latitude increases. In order to mitigate this dependence on the latitude, we applied the following procedure. We subtracted the theoretical values of the rotation rate (equation 1) from the rotation rate of each individual tracer and analyzed the distributions of the differences. These differences represent the deviation of the rotation rate of each tracer from some mean rotation rate at a particular latitude. The theoretical rotation rate was calculated using equation 1 with the constants derived from magnetic field data. Note that we could use the theoretical curve for the tracer rotation rates derived from continuum intensity data. The latter choice would result only in a relative shift of the distributions with respect to zero expected value. The distributions of the calculated differences for the rotation rates derived from magnetic field and continuum intensity data are shown in Figure 4 in orange and light blue colors, respectively. Red and blue curves represent the best fittings of the distributions by a Gaussian. The centers of the Gaussians are located at −0.023±0.004 and (0.391±0.028) ∘ day −1 for red and blue curves, respectively. The width of the Gaussians are 0.206±0.004 (red) and (0.626±0.028) ∘ day −1 (blue). A slight asymmetry in both distributions are clearly seen in Figure 4. The right-hand wings of the distributions are heavy-tailed, i.e. there is an excess of tracers with the rotation rates exceeding mean value. This asymmetry is quite similar to that found by Nagovitsyn et al. (2018) for the rotation rate of sunspot groups (see Figure 1 in that paper). The authors argued that the asymmetry is caused by the existence of two populations of sunspot groups that exhibit different behaviour of the rotation rates. Hence, the distributions in Figure 4 might represent a sum of two Gaussians (see Figures 3 and 5 in Nagovitsyn et al. 2018). We suppose that our data are too scanty to perform a reliable fitting of the distributions by two Gaussians.  Figure 4. The distributions of the differences between the theoretical rotation rate (derived from magnetic field data) and the measured rotation rates derived from magnetic field (orange) and continuum intensity (light blue) data (see text). Red and blue curves represent the fittings of the distributions by a Gaussian. A slight asymmetry of the distributions is clearly seen.

Conclusions and discussion
In this study we used SDO/HMI data to probe the rotation rate of 670 active regions observed between 2010 and 2016. The rotation rates were measured using continuum intensity images and line-of-sight magnetograms of the same active regions. For continuum intensity images, the center of the active region was defined as the intensity-weighted position of the entire active region patch. For magnetograms, the flux-weighted centroids of the opposite magnetic polarity were calculated separately. The center of the tracer was set as the geometrical center of the centroids. The difference in the calculation of the active region center is the only distinction in the applied methods. Continuum intensity images and magnetograms were handled in the same way at all other steps of the data processing.
We found that continuum intensity images yield systematically higher values of the derived rotation rate as compared to magnetograms. The ω ci sid for bi/multipolar active regions exceeds ω m sid by 0.45 ∘ day −1 on average. This value is comparable to the difference in the rotation rate near the solar equator and at 30 ∘ latitude. In other words, the uncertainty in the measured rotation rate of a tracer is comparable to the variation of the solar differential rotation within the ±30 ∘ wide latitudinal zone, i.e. the zone where most of active regions appear during a solar cycle.
We attribute the difference in the derived rotation rates exclusively to the methodological issues. In our opinion, the apparent motion of the intensity-weighted center of a sunspot group toward a more coherent leading part is the main reason for the discrepancy. This issue was pointed out by Petrovay (1993) for the first time. However, surprisingly little attention was paid to that work.
The controversial results on the rotation rate variation of sunspot groups mentioned in the Introduction might be also explained by the data utilized in previous analyses. Indeed, it is obvious that the apparent acceleration or deceleration may be observed during different phases of evolution of a sunspot group. Figure 1 is a good visual explanation of this issue. The intensity-weighted center accelerates toward the West as the following sunspot starts to decay quickly after emergence (see three top left-handed panels in Figure 1). As the following part disappears completely, the intensity-weighted center coincides with the position of the leading sunspot. Obviously, by this time the rotation rate of the intensity-weighted center decelerates down to the rotation rate of the leading sunspot (see two top right-handed panels in Figure 1). Therefore, the variations of the rotation rate measured by an observer depends exclusively on the captured evolution phase of the sunspot group.
As we found in this work, using magnetic field maps results in more robust estimation of the proper rotation rate of an active region. Magnetic field maps are more suitable for analyzing the geometry of an active region as compared to continuum intensity images. On the other hand, the undeniable advantage of the continuum intensity data is the vast time span of the existing databases covering 150 years of observations. These data allow one to analyze time variations of the rotation rate within a dozen of solar cycles (e.g. Pulkkinen and Tuominen 1998). Unfortunately, the databases on full-disk magnetograms are not as large as those of sunspots. However, the increase in the spatial and temporal resolutions in modern magnetic field observation allows one to probe the dynamics of solar active regions in unprecedented detail. As a final remark, we suggest using magnetic field data rather than white-light images in future works on solar differential rotation. ence Team. The study was supported by the Russian Science Foundation, Project 18-12-00131.
Author contributions: All authors have accepted responsibility for the entire content of this manuscript and approved its submission.

Conflict of interest:
The authors state no conflict of interest.