The Sun’s Magnetic Power Spectra Over Two Solar Cycles. I. Calibration Between SDO/HMI And SOHO/MDI Magnetograms

In our series of studies, we use two radial synoptic magnetogram data sets covering the 23rd and 24th solar cycles. They are observed by Michelson Doppler Imager (MDI) on board the Solar and Heliospheric Observatory (SOHO) (Scherrer et al., 1995) and the Helioseismic and Magnetic Imager (HMI) on board the Solar Dynamics Observatory (SDO) (Scherrer et al., 2012; Schou et al., 2012), respectively. The two instruments utilize different spectral lines with various resolutions: the Ni I 6768 Å (Fe I 6173 Å) (Norton et al., 2006) is adopted as the spectral line for the MDI (HMI), with coverage of $3600\times 1080\ (3600\times 1440)$ pixels. The data sets begin in Carrington Rotation (CR) 1911 (July 1996) and end in CR 2225 (December 2019, data missing for CRs 1938-1942 and CR 1945). In this paper, we focus on analyzing the synoptic maps for the overlap period, i.e. from CR 2097 to CR 2104.

The initial synoptic maps lack data on the polar magnetic fields, which are crucial for the decomposition of the spherical harmonic function. Therefore, we use polar-corrected data (Sun et al., 2011; Sun, 2018). Specifically, they interpolate and smooth the north (south) polar magnetic fields by using a multi-year series of well-observed polar fields from each September (March). This process enables synoptic maps with polar field correction to be obtained. It is worth noting that the polar correction approach of HMI differs slightly from that of MDI, and the implications of these differences are analyzed in Section 3.2.

For the magnetic field $B(\theta,\phi)$, we can express it by the following Eq. (1)

where $\theta$ is the colatitude, $\phi$ is the longitude, $Y_{l,m}(\theta,\varphi)$ is the spherical harmonic function with degree $l$ and azimuthal order $m$, and $B_{l,m}$ is the corresponding spherical harmonic coefficient.

For each synoptic map, we use the “pyshtools” package in Python to perform a spherical harmonic decomposition. The algorithm is proposed by Driscoll & Healy (1994), who employ a grid suitable for data with different sampling densities at different latitudes and introduce weight coefficients. For example, the longitudinal distances of the pixel points near the equator are larger in synoptic maps than those near the polar region. Therefore, the data need to be weighted when performing the spherical harmonic decomposition. The discretization formula for obtaining the spherical harmonic coefficient for a grid of size $2n\times 2n$, introducing the weighting factor $a_{j}$, is given by:

where $\theta_{j}=\pi j/2n$, $\phi_{k}=\pi k/n$, corresponding to spherical coordinates for any point on the sphere, and the number of grid points in both latitude and longitude is $2n$. The weight coefficient $a_{j}$ is given by the following equation

This algorithm requires the grid to be $n\times 2n$ or $n\times n$. Therefore, as the first step in data processing, we transform the data resolution to $2160\times 1080$ for MDI and $2880\times 1440$ for HMI by a simple linear two-point interpolation algorithm. Afterward, we also convert data from an equal sine-latitude distribution to an equal latitude distribution. At this step, we first calculate the sine of latitude at equal intervals. Then the same two-point interpolation algorithm is used to transform the data based on sine-latitudes. After these processes, the data can be used for “pyshtools”.

Due to the instrumental resolution limitations and sampling theorem, we decompose the data up to $l=539$ and $l=719$ for MDI and HMI synoptic maps, respectively. In the subsequent analysis, we consider only spherical harmonic decomposition results of $l<540$ for both data sets. As an example, Figure 1 presents the magnetic field components at different spatial scales after decomposing the CR 2097 synoptic map, which is shown in Figure 2 (a). Figure 1 (a) is the spherical projection of the original synoptic map. Panels (b)-(d) show maps reconstructed by spherical harmonics up to different maximum degrees. Large-scale structures become progressively more refined from Panels (b) to (d) with $l_{max}$ increasing from 30 to 539. These expected results illustrate the algorithm’s validity for the spherical harmonic decomposition of the magnetograms. We can therefore use the results to obtain the power spectra.

Using the orthogonality between the spherical harmonic functions, Parseval’s theorem in the Cartesian coordinate system can be extended to the spherical geometry. That is,

where $P$ is the power spectrum and it is related to the spherical harmonic coefficients by

It is clear from Eq. (5) that in the power spectrum, information about the scale of azimuthal order $m$ is combined. In our work, the spectral power represents the magnetic energy at the spatial scale corresponding to the spherical harmonic degree $l$. To convert the degree $l$ to the actual spatial scale $\lambda$, we can use the following equation:

where the parameter $R_{\odot}$ is the solar radius.

After the above data processing steps, we obtain the magnetic power spectra of eight CRs during the overlapping period. As examples, Figures 2 (a) and (b) show the synoptic magnetograms for CR 2097 and CR 2101, respectively, while the corresponding HMI and MDI magnetic power spectra are shown in Figures 2 (c) and (d). The same figures but for the remaining six CRs are shown in the appendix.

Figures 2 (c) and (d) show that both magnetic field strength and spatial scales influence the difference between HMI and MDI power spectra, which is consistent with previous studies (Liu et al., 2012; Riley et al., 2014; Pietarila et al., 2013; Virtanen & Mursula, 2017). We next investigate whether the impact of these two factors is minimal enough to use a fixed calibration factor or whether it is necessary to calibrate based on spatial scales.

We present a comparison of the MDI and HMI spectral power for data during the overlapping period in Figure 3. The red line represents the calibration results obtained by Liu et al. (2012). The calibration factor, $1.4$, corresponds to the vertical intercept in the log-log plot. The plot shows that the distribution of points is non-linear, with a significant deviation from the red line. That is, the calibration factor varies enough to make a fixed calibration factor inapplicable. Thus we need to develop a new calibration method applicable to power spectra analysis.

Notably, black dots in Figure 3 (small-scale part of magnetic fields) generally have smaller calibration factors and lower spectral power than grey dots (large-scale part). This may suggest that the calibration factors are smaller for weaker magnetic fields. However, this seems to contradict the findings of Liu et al. (2012), which indicate that the calibration factors are smaller for stronger magnetic fields. This discrepancy may arise from the stronger influence of spatial scale on the calibration factor than that of magnetic strength, with smaller spatial scales corresponding to smaller factors. These results indicate scale dependency in the difference between MDI and HMI magnetograms. Therefore, it is essential to consider the calibration at various scales in our work.

In the theory of spherical harmonic decomposition, the term $l=0$ represents the average result of the global scale magnetic field. For the magnetic field, the term should be zero theoretically. However, this coefficient is not zero based on observed magnetograms, which may be introduced by noise or during the construction of the synoptic magnetograms (Bertello et al., 2014). We will not consider this term in the analysis that follows. The power spectra in this paper all start from $l=1$.

In Section 2.1, we mention that Sun et al. (2011); Sun (2018) use similar approaches to correct the polar magnetic field of the HMI and MDI synoptic magnetograms. However, the approaches require the use of the previous year’s data for interpolation. But HMI synoptic magnetograms for the overlap period do not have the previous year’s data. So for this part, they use the data of MDI instead. In addition, the polar fields are smoothed and weighted with the non-polar data in the correction. These correction methods increase uncertainties of the polar magnetic field and further influence the magnetic power spectra. Therefore, we need to determine the calibration range according to the range of the power spectra affected by the polar field.

We multiply the data above 75° north and south of MDI magnetograms by a factor of 0.8 to test the effect of the polar field. We then obtain their power spectra, which are shown in Figures 2 (e) and (f). It can be seen that the scaled MDI spectral power corresponding to the dash line has a certain decrease at $l=1,3,5$ and are close to HMI. This just supports the idea that the polar magnetic field has an impact on the power spectra of the large scale field ($l\leq 5$). The results are similar to Virtanen & Mursula (2016): The correction of the polar field affects axial dipole ($l=1$) and quadrupole ($l=2$) but has no effect on the results of a smaller scale. Thus, the effect of the polar field uncertainty for $l>5$ can be ignored. Taking these results together, we will calibrate the power spectra for $l>5$ and the calibration results will be applied in the non-overlap period.

In the first two subsections, we have shown that the calibration factors are scale-dependent in the range $l>5$. To investigate this further, we combine data for the eight CRs with $l>5$ and explore the pattern of variation of the spectral power ratio $B^{2}_{l-MDI}/B^{2}_{l-HMI}$ with the spherical harmonic degree $l$. The result is shown in Figure 4. The plot clearly shows a decreasing trend in the ratio as the degree $l$ increases and this trend is consistent across different CRs. This further suggests that the calibration factor is a function of $l$.

Due to the limited number of CRs with the same $l$, the impact of noise or measurement error is relatively large. To mitigate the impact, we exclude special points by discarding data larger than the variance $\sigma$ for a given $l$. The remaining data points are then averaged to obtain the solid blue line in Figure 4. We use a nonlinear function to fit the averaged data and obtain the following calibration function $r(l)$:

The calibration factor $r$ for $l=5\ (\approx 800$Mm) is 1.4, while for $l=30\ (\approx 140$Mm) and $l=140\ (\approx 30$Mm), the calibration factors are 1.35 and 1.23, respectively. The factors are close to the calibration factor of 1.4 obtained in Liu et al. (2012). These demonstrate that the calibration factors obtained through traditional methods for analyzing full-disk magnetograms are mostly relevant for large-scale structures, such as active regions. Our results suggest that the classical methods may have limited applicability to smaller-scale features. In contrast, Riley et al. (2014) found a calibration factor $r\approx 1.05$ for the CR 2097 synoptic magnetogram. This value is smaller than our calibration factor for the large-scale section but larger than that for the small-scale section. As shown in Figure 2 (c), small-scale features contribute significantly to the total power, which suggests that the overall single calibration factor obtained by averaging our calibration factors may be similar to the value of 1.05.

In addition, Virtanen & Mursula (2017) concluded that the calibration factors for low spherical degrees are smaller than those obtained through pixel-by-pixel comparison or histogram techniques when they performed spherical harmonic decomposition for synoptic maps. However, it should be noted that their work did not involve the calibration of HMI and MDI data sets, thus the applicability of their conclusion to the calibration of these two data sets is uncertain. Hence, we should be careful to directly compare our findings with theirs.

After applying the calibration function Eq. (7) to the power spectra, we first evaluate the similarity between the calibrated HMI and MDI power spectra. Specifically, we consider the calibrated HMI and MDI power spectra to be consistent if the difference between them is less than $20\%$ of the MDI spectral power. We list the consistent interval of $l$ based on this criterion in the second column of Table 1. Our calibration method leads to throughout consistent HMI and MDI power spectra for CRs 2099-2103, based on this criterion. Moreover, the power spectra of CR 2097 and CR 2098 exhibit consistency only on the large-scale part, from $l=6$ to around $l\approx 200$. And the consistent interval for CR 2104 is mainly in the small-scale part.

The exception of CR 2104 is also apparent in the corresponding images. In Figure 3, below the red line, there is a cluster of gray dots, which results from CR 2104. In Figure 4 gray points below the black dots show a significant deviation. This part of the data is also from CR 2104. These results suggest that there might be some issues with the CR 2104 HMI or MDI magnetograms, which could be attributed to instrumental measurements.

To illustrate in more detail the applicability of our calibration method for the power spectra analysis, we additionally use two quantitative comparisons. The residual sum of squares ($RSS$) measures the difference between the MDI power spectra and the uncalibrated ($r=1$) or calibrated HMI power spectra. The formula is as follows.

The results are shown in columns 3 to 5 of Table 1, where $r=1.4$ represents the fixed calibration factor 1.4 and $r(l)$ represents the calibration function Eq. (7). The smaller the RSS value is, the smaller difference between the two power spectra is. While using a fixed calibration factor, the RSS is reduced to only one-fifth to one-half of the uncalibrated case. For CR 2098 and CR 2104, the RSS even shows an increase. However, when using the calibration function $r(l)$, all the RSS show a decrease, with the vast majority reduced to a tenth or even a hundredth of the uncalibrated case. This shows that our calibration method is more suitable for power spectra analysis.

In addition to using the RSS to evaluate the effect of our calibration on the power spectra, we also quantify the effect on the magnetograms. To accomplish this, we perform an inverse spherical harmonic decomposition of the power spectra to reproduce the magnetograms. Using the MDI magnetograms as the reference, we separately obtain the HMI magnetogram using the two calibration methods (fixed calibration factor 1.4 and calibration function $r(l)$). The ratio of the difference between the HMI and MDI magnetograms at a given pixel is obtained using the following formula:

For each CR map, we then conduct a linear regression of the results for all pixels. The results are summarized in the last column of Table 1. Notably, the ratios of magnetogram differences for all CRs are within a narrow range of $0.81\pm 0.01$. These quantitative comparisons indicate the efficacy of our calibration approaches in achieving not only a better analysis of power spectra but also a more consistent between HMI and MDI synoptic maps.

The analysis of Dopplergrams using spherical harmonic functions, as discussed in Introduction, shows the peak around $l=120$ on the kinetic power spectra corresponding to supergranulation. In Section 3.3 we have mentioned that MDI and calibrated HMI power spectra are consistent around this degree. So it is feasible to search for supergranule features in the magnetic power spectra.

Based on Figure 2 and the power spectra images in the appendix, it is apparent that a majority of the magnetic power spectra exhibit peaks or knees around $l=120$. Here we use three algorithms to verify the existence of these structures and the corresponding specific spatial scales. We start with a morphological approach based on the code developed by Duarte & Watanabe (2021) to find the peaks in the spectra. We search and identify peaks by the trend of the curve around the local maximum, the threshold, and the width. By adjusting these parameters, the identification results can be optimised. The corresponding results are in columns 2 and 3 ($l_{sgr1}$) of Table 2. The second method is similar to wavelet analysis, where a wavelet is fitted to the peak and the presence of a peak is discerned based on the result of fit (Du et al., 2006). This method is achieved through the ‘scipy.signal.find_peaks_cwt’ within Python. Since the method does not give the peak, but gives the center of the wavelet, we list the range of the peak rather than the exact value, with the results in columns 4 and 5 ($l_{sgr2}$) of Table 2. The third method we use is the algorithm for finding peaks and knees used in the field of neural power spectra: FOOOF (Donoghue et al., 2020). This method gives information on the peaks of the four CRs around the supergranule scale and also gives the knee of the three CRs. The results are in columns 6 and 7 ($l_{sgr3}$) of Table 2.

The supergranulation identification results are between $l=116\ (36$ Mm) and $l=140\ (30$ Mm), which are consistent with the results in the kinetic power spectra. Notably, for the same CR, we observe slight differences in the results obtained from different methods. For instance, in the case of CR 2099, the first two methods yield a difference of $\Delta l=16$. The third method, meanwhile, can assist in determining which result is more reliable. Therefore, by combining the three methods, we can identify the supergranule scale correlated features of the power spectra of most CRs. However, we are unable to identify significant features in the power spectra of CR 2103, possibly due to the portion of power spectra belonging to active regions masking supergranulation.

Furthermore, in the work of Williams et al. (2014), they observe different sizes of supergranulation in the kinetic power spectra for HMI and MDI, with $l=132$ for HMI and $l=119$ for MDI. They attribute the difference to the different spatial resolutions. They show that by smoothing the resolution of HMI Dopplergrams, a power spectrum similar to MDI data can be obtained. However, our identifications do not indicate any significant differences between the HMI and MDI results, which remain generally consistent.

Regarding the absence of differences between the supergranular scales in our results mentioned above, we first exclude the effect of the calibration. We note that the difference of calibration factors ($\Delta r=0.02$) around the supergranule scale cannot significantly alter the scale of the feature. Furthermore, the identified supergranule sizes in HMI and MDI maps of CR 2104 remain consistent, even though their power spectra do not coincide in $l<162$. We suggest that the difference in the identification results is due to the influence of small-scale structures. Hathaway et al. (2002); Rincon et al. (2017) mention that the kinetic power spectra near the size of supergranulation can be affected by granulation, which are not completely removed. The higher the resolution, the more features of granulation are included. Therefore, smoothing maps with high resolution is essentially a way to obtain similar power spectra by reducing the influence of small-scale structures. Magnetic power spectra in this paper are obtained from synoptic magnetograms, which have a relatively low resolution. Multiple full-disk magnetograms are averaged in the construction of a synoptic map. Strong small-scale magnetic structures and noise are removed. Therefore, we can obtain ideal results in the magnetic power spectra. On the other hand, in addition to small-scale structures, large-scale structures such as active regions may also cause an unphysical shift of the position of supergranulation in the power spectra. We need to be careful about this potential effect in future work.