© K. Nagashima et al. 2020

1. Introduction

Understanding solar interior dynamics is crucial to understanding the mechanisms of the solar dynamo. As one example, convection may play an important role in the formation of magnetic flux tubes, as well as in their rise through the convection zone and their tilts at the solar surface (e.g., Brun & Browning 2017). Helioseismology, which uses observation of oscillations on the solar surface, is an important probe of interior dynamics.

Currently there is a major discrepancy between the time-distance (Hanasoge et al. 2012) and ring-diagram (Greer et al. 2015) estimates of the strength of solar subsurface convection at large spatial scales. The measurements from time-distance helioseismology suggest flows orders-of-magnitude weaker than those seen in convection simulations (e.g., in the anelastic spherical harmonic (ASH) convection simulations by Miesch et al. 2008), while the measurements from ring-diagrams are closer to the expectations from simulations.

Time-distance helioseismology (Duvall et al. 1993; Kosovichev & Duvall 1997) is based on measuring and interpreting the travel times of acoustic and surface-gravity wave packets. These travel times are measured from the temporal cross-covariances between the Doppler observations at pairs of points on the solar surface (see Gizon & Birch 2005, for a review). Ring-diagram analysis (Hill 1988; Antia & Basu 2007) measures the Doppler shift of acoustic and surface-gravity oscillation modes in the local power spectra and uses these Doppler shifts to infer local flows in the solar interior. To compute the local power spectra, the solar surface is divided into a number of spatial tiles. Each tile is tracked at a rate close to the solar rotation rate. The power spectra of the solar oscillations (Doppler observations) are computed for each tile. In the three-dimensional power spectra, roughly concentric rings with high power are present at each frequency and they correspond to the modes of different radial orders. Flow and wave-speed anomalies in the Sun shift and distort the rings and, hence, one can obtain information about the solar interior from the ring parameters. Flow maps from Doppler observations by the Helioseismic and Magnetic Imager (HMI; Schou et al. 2012) onboard the Solar Dynamics Observatory (SDO; Pesnell et al. 2012) are automatically computed by the SDO/HMI ring-diagram pipeline (Bogart et al. 2011a,b) on a daily basis since the SDO launch in 2010.

The HMI ring pipeline codes separately fit each single ridge (single radial order n) in the power in slices at constant horizontal wavenumber or slices in temporal frequency. Greer et al. (2014) developed an alternative approach based on simultaneously fitting multiple ridges (multiple radial orders) at each horizontal wavenumber. Greer et al. (2015) introduced another innovation: they chose a much denser tile layout than other ring-diagram analyses: 16-degree tiles with the separation of 0.25 degree instead of the 7.5 degree spacing (at the equator, the spacing increases at higher latitude to maintain 50% overlap between neighboring tiles) used for 15-degree tiles in the HMI ring pipeline. As described in detail in Appendix A.2, the ATLAS code from Greer et al. (2014, 2015) provides seven parameters for each ridge and five parameters for the background power at each wavenumber k. Greer et al. (2015) applied a three-dimensional flow inversion to these fit results to estimate the three-dimensional flow field in the solar interior. Both the dense packing of tiles and the three-dimensional inversions are unique to ATLAS, however, in this paper we focus only on the fitting component of the code.

As one step toward understanding the causes of the above-mentioned disagreement between the helioseismic measurements of subsurface convection, here we focus on the ring-diagram analysis described by Greer et al. (2015). We revisit the analysis code (Greer et al. 2014; Greer 2015) that was used in that work and identify several issues through a step-by-step examination of the code. In response to these issues, we have developed an updated method and characterize the updated code by applying it to synthetic data generated from Monte-Carlo simulations.

2. Description of the updated code

In this section, we describe our updated ATLAS code. Each of the updates is a response to a problem that we found in our inspection of the original code. In this section we will refer to Appendix A for the details of the original code.

After computing the power spectrum for a particular tile, the processing steps are: (1) remap the power spectrum from Cartesian (k_x, k_y) to polar (k, θ) coordinates, (2) rebin the power in azimuth θ; the number of grid points in θ is reduced from n_pix = 256 to n_pix = 64, (3) fit the logarithm of a Lorentzian model to the logarithm of the smoothed power by least-squares minimization at each k using the Levenberg–Marquardt (Marquardt 1963) technique; the model function has 7n_r + 4 parameters at each horizontal wavenumber k, where n_r is the number of ridges at the particular value of k, and (4) estimate the covariance matrix of the errors of the fitted parameters by computing the inverse of the Hessian matrix of the cost function. In the following subsections we describe the changes that we introduced to each of these steps.

2.1. Re-binning in azimuth

The computed local power spectra O(k_x, k_y, ν) and the interpolated spectra in polar coordinates are non-negative by construction. Following the original code, the interpolated spectra have 256 pixels at each k. At k R_⊙ ≡ ℓ ∼ 500 (where R_⊙ is the solar radius), which corresponds to k_pix ≡ k/h_k = 21 (where h_k = 3.37 × 10⁻² Mm⁻¹ is the grid spacing in k) and which we use in most of the Monte-Carlo test calculations shown later, this is about twice of the number of grid points in θ from the full resolution, which has n_pix ≈ 2πk/h_k. For the sake of computational efficiency, it is desirable to reduce the number of points in azimuth. In the updated code, we use a running box-car smoothing of four-pixel width followed by subsampling by a factor of four to reduce to the number of grid points in θ. This procedure ensures that the resulting smoothed power spectrum O_s(k, θ, ν) will be positive. We expect that as flows produce θ variations that are dominantly at azimuthal wavenumber one, it should be possible retain only very low resolution in θ; Sect. 4.1 discusses this in more detail. We retain 64 points in θ for the examples shown in this paper.

The original code used a low-pass Fourier filter to smooth the power spectrum in the θ direction. This procedure produces occasional points where the smoothed power is negative.

2.2. Least-squares fitting of the logarithm of power

In the updated code, we use a least-squares fitting to fit the logarithm of the model power to the logarithm of the remapped and smoothed power. Following Greer et al. (2014), the fitting is carried out independently at each horizontal wavenumber k. As the amount of smoothing is increased, the probability distribution function (PDF) of the power, as well as its logarithm, approaches a normal distribution (see Appendix B.4 for more details), and it is therefore appropriate to use a least-squares fit. The logarithm of the smoothed power has the convenient property that the variance of the logarithm of power (σ_N in Eq. (B.4)) depends only on the details of the remapping and smoothing and does not depend on θ or ν (see Appendix B.2 for details).

In our approach, the cost function at a single k is:

$$\begin{array}{c}\hfill -lnL(\mathit{q})={\displaystyle \sum _j{[ln{O}_s({\theta }_j,{\nu }_j)-lnP({\theta }_j,{\nu }_j;\mathit{q})]}^2,}\end{array}$$(1)

where O_s(θ, ν) is the observed spectrum at some fixed k after smoothing in θ and P(θ, ν; q) is the model of the spectrum with model parameters q. The summation is taken over all bins j within the fitting range. We note that as the error estimates for lnO_s (σ_N in Eq. (B.4)) are all the same, we set them all to one in writing the cost function. For the sake of readability, throughout the remainder of this paper we do not introduce or carry notation to denote the value of k; the fitting problems at each k are treated as independent and we will not be comparing fit parameters for different values of k. This is an approximation as the interpolation from (k_x, k_y) space to (k, θ) space does imply error correlation between the fit parameters at different values of k.

In the updated code, we use the Levenberg–Marquardt technique to solve the minimization problem. In particular, we use mpfit.c (Markwardt 2009), one of the codes in the MINPACK-1 least-squares fitting library (Moré 1978; Moré & Wright 1993). The covariance matrix of the errors associated with the fitted parameters are estimated at the last step of the Levenberg-Marquardt procedure. As a practical note, the error estimates obtained by this method must be scaled as we have assumed σ = 1 in Eq. (1); see original PDF, Eq. (B.4). In general, calling the code with an incorrect estimate of σ could cause poor performance of the fitting algorithm. In the current case σ is not far from one (σ ∼ 0.5, see Sect. 3.1) and we do not expect that this is a significant issue here. To implement the σ estimation in the code is a task for the future.

The original ATLAS code used a maximum-likelihood method based on the assumption that the power spectrum in a single bin in (k_x, k_y, ν) space follows a chi-squared distribution with two degrees of freedom; see Appendix B.1 for the original likelihood function. This method does not account for the impact of smoothing on the PDF of the power spectrum.

2.3. The functional form of the power spectrum

The model function for the power spectrum at a single k in our updated code is

$$\begin{array}{cc}& P(\theta ,\nu ;{\mathit{q}}_{\mathrm{BG}})=\hfill \\ \hfill & {\displaystyle \sum _{n=0}^{n_{\mathrm{r}}-1}\frac{q_{1,n}(q_{2,n}/2)F(\theta ;q_{5,n},q_{6,n})}{{[\nu -q_{0,n}+k(q_{3,n}cos\theta +q_{4,n}sin\theta )/(2\pi )]}^2+{(q_{2,n}/2)}^2}}\hfill \\ \hfill & +B(\theta ,\nu ;{\mathit{q}}_{\mathrm{BG}}),\hfill \end{array}$$(2)

where

$$\begin{array}{c}\hfill F(\theta ;p_1,p_2)=1+p_{1}cos(2\theta )+p_{2}sin(2\theta ).\end{array}$$(3)

The fitting parameters for the peaks (n = 0, …n_r − 1, where n_r is the number of the ridges) are q_0, n = ν_n is the frequency of the nth peak, q_1, n = A_n is the amplitude, q_2, n = Γ_n is the width, (q_3,n, q_4,n) = u = (u_x,n, u_y,n) is the horizontal velocity, q_5, n = f_c, n and q_6, n = f_s, n are parameters to handle anisotropy. The background is modeled at each k as

$$\begin{array}{c}\hfill B(\theta ,\nu ;{\mathit{q}}_{\mathrm{BG}})=\frac{q_{0,\mathrm{BG}}}{{\nu }^{q_{1,\mathrm{BG}}}}F(\theta ;q_{2,\mathrm{BG}},q_{3,\mathrm{BG}}),\end{array}$$(4)

where F is again given by Eq. (3) and q_0, BG = B₀ is the amplitude, q_1, BG = b is the power-law index, q_2, BG = f_c, bg and q_3, BG = f_s, bg are the parameters to handle the anisotropy. The number of the parameters in total is 7n_r + 4 at each k. For reference, Tables C.1–C.3 show the physical meaning of each of the fitting parameters.

We altered the model function Eq. (2) from the original Eq. (A.1) to make the parameterization more stable. The following subsections describe the motivation for these changes.

2.3.1. Parameters of the anisotropy terms

We replaced F′ defined by Eq. (A.2) with F defined by Eq. (3) and used the same form for the anisotropy terms for the background, although the exact form for the background function is subject to other alterations discussed in Sect. 2.3.2.

Our alteration does not change the space of functions that can be fit with F(θ) but it does not suffer from the issue of indeterminate phase for nearly isotropic power spectra. In the original parameterization, the amplitudes of the anisotropic part of the model function ($q_{5,n}^{\prime }$ and $q_{3,\text{BG}}^{\prime }$) are usually much smaller than one. As a consequence, the phase ($q_{6,n}^{\prime }$ or $q_{4,\text{BG}}^{\prime }$) does not matter much in the fitting, and, hence, it is unstable. In the particular case of isotropic power spectra with $q_{5,n}^{\prime }=0$ for all n and $q_{3,\text{BG}}^{\prime }=0$, the phases $q_{6,n}^{\prime }$ and $q_{4,\text{BG}}^{\prime }$ are indeterminate.

2.3.2. Parameters of the background model

We also changed the background parameterizations for the sake of stability. The background model (Eq. (A.3)) originally contained five parameters. In this section we explain our motivation for reducing this to the four parameters shown in Eq. (4).

The background model in Eq. (A.3) is based on the model of Harvey (1985):

$$\begin{array}{c}\hfill B_{\mathrm{Harvey}}(\nu )=\frac{4{\sigma }_{\mathrm{rms}}^2\tau }{1+{(2\pi \nu \tau )}^2},\end{array}$$(5)

where τ is the characteristic timescale of the velocity field in question and σ_rms is the rms velocity. In this case the index in the original background model (Eq. (A.3)) would be $q_{2,\text{BG}}^{\prime }=2$. Appourchaux et al. (2002) suggested a generalization of the Harvey (1985) model, where $q_{2,\text{BG}}^{\prime }$ can be the range of 2 to 6.

In the ATLAS fittings, however, we typically obtain the index $q_{2,\text{BG}}^{\prime }~1$ for HMI observations. Also, the roll-off frequency obtained from the fitting of HMI observations, $q_{1,\text{BG}}^{\prime }$, is quite low (∼1 μHz) and below the frequency resolution for 28.8 h power spectra (9.7 μHz), which is the typical observation length for 15-degree tiles in HMI pipeline and 16-degree tiles in ATLAS.

This suggests that the background is not related to either the supergranulation (τ ≈ 10⁵ s, or ν ≈ 10 μHz) or granulation time scales (τ ≈ 4 × 10² s, ν ≈ 2.5 mHz), based on Table 1 of Harvey (1985). As these background parameters are not consistent with the original physical model, we might need to reconsider the background model. At this moment, however, this is a task for the future, and we retain this model in the altered form mentioned below.

In the case of $\nu \gg q_{1,\text{BG}}^{\prime }$, the original background model, Eq. (A.3), can be simplified:

$$\begin{array}{cc}& \frac{q_{0,\mathrm{BG}}^{\prime }}{1+{(\nu /q_{1,\mathrm{BG}}^{\prime })}^{q_{2,\mathrm{BG}}^{\prime }}}{F}^{\prime }(\theta ;q_{3,\mathrm{BG}}^{\prime },q_{4,\mathrm{BG}}^{\prime })\hfill \\ \hfill & \sim \frac{q_{0,\mathrm{BG}}^{\prime }{q_{1,\mathrm{BG}}^{\prime }}^{q_{2,\mathrm{BG}}^{\prime }}}{{\nu }^{q_{2,\mathrm{BG}}^{\prime }}}{F}^{\prime }(\theta ;q_{3,\mathrm{BG}}^{\prime },q_{4,\mathrm{BG}}^{\prime }),\hfill \end{array}$$(6)

we therefore redefine the background model B(θ, ν) with four parameters q_i, BG(i = 0, 1, 2, 3) with the altered anisotropy terms in the form of Eq. (4). As the new background model (q_i, BG) has only four parameters instead of the original five, the index i has different meanings in the two models.

3. Performance of the updated code

Monte-Carlo simulations are a powerful tool for testing fitting methods. In this section, we use this approach to characterize the performance of the updated fitting code. In Sect. 3.1 we compare the scatter of the fitting results of Monte-Carlo simulations with the noise estimated by the updated code. In Sect. 3.2 we measure the bias of the flow estimates for some simple cases.

3.1. Error estimates

We use the approach of Gizon & Birch (2004) to generate realizations of the wavefield. The assumption of this approach is that the real and imaginary parts of the wavefield at each point (k_x, k_y, ν) are independent Gaussian random variables. In more detail, the method is as follows: (1) create a Lorentzian model (Eq. (2)) for the limit spectrum using a set of input parameters, (2) pick two standard normally distributed random numbers at each grid point (k_x, k_y, ν) and take sum of the squared numbers divided by two to make a chi-square distribution with two degrees of freedom, and (3) multiply the limit spectrum by these random numbers to obtain one realization of the power spectrum.

After each realization of the power spectrum is generated in (k_x, k_y, ν) space, it is then remapped and rebinned in the same manner as the local power spectra computed from the observations. To obtain the parameters for the input model power spectrum in the first step above, we used the average over Carrington rotation 2211 of the power spectra for the disc center tile from the SDO/HMI pipeline (these average power spectra were obtained using the Data Record Management System (DRMS) specification hmi.rdvavgpspec_fd15[2211][0][0]). We then used the updated code to fit these average power spectra. The parameters resulting from these fits at ℓ = 328, 492, and 984 (k_pix = 14, 21, and 42, respectively) are shown in Tables C.1–C.3, respectively. We later use these parameters to generate Monte-Carlo synthetic data.

Figure 1 shows a few slices of the input limit spectrum, namely, Eq. (2) evaluated for the parameters q listed in Tables C.1–C.3. The original observed power spectrum is also shown. Figure 1 shows that the model power spectrum is reasonable compared to the observables, although the details of the peak shapes and the background slopes are not always reproduced by the model and there is still room for improvement with regard to the model function.

Fig. 1. Slices through the input limit spectrum (black) and observational spectrum averaged over one Carrington rotation (gray), which was used to obtain the input parameters, at θ = 0 and a few different values of k. Vertical dashed lines indicate the fitting ranges and the models are plotted only within these ranges. Lower limits for all k are fixed at 0.4 mHz, while the upper limits depend on the initial guess for each k, and they are the highest peak frequencies plus the widths of the peak. For the cases of ℓ = 328 and 492 (kpix = 14 and 21, respectively), the peaks with the radial order, n, from 0 to 7 are used in the fitting, while for the case of ℓ = 984 (kpix = 42), 0 ≤ n ≤ 3.

We created 500 realizations of the power spectrum using the above procedure. The fitting code occasionally produces outliers and for these computations, we removed them. Specifically, we removed outliers iteratively: in each iteration we computed the standard deviation of the samples between the tenth and 90th percentile and we discarded points at more than five times the standard deviation away from the mean. We repeated this procedure until no further points were removed. After this outlier removal, the number of valid samples are 489 (ℓ = 328, k_pix = 14), 500 (ℓ = 492, k_pix = 21), and 496 (ℓ = 984, k_pix = 42) out of 500 samples.

Figure 2 shows the average and scatter of the output peak parameters associated with the n = 0…5 ridges for ℓ = 492 (k_pix = 21) from fitting 500 Monte-Carlo realizations of the power spectrum. Table 1 shows the average and scatter of the output background parameters from the same set of fitting results. From Fig. 2 and Table 1, we see that most of the averages of the parameters estimated by the fitting are within the expected scatter in the mean (${\sigma }^{\mathrm{scat}}/\sqrt{N_{\mathrm{sample}}}$, where N_sample is the number of Monte-Carlo realizations, hence N_sample = 500 here). The amplitude is an exception; it is always smaller than the input. This is the result of taking the logarithm of the smoothed power; see Appendix D for details. Also, from Fig. 2 and Table 1, we see that the error estimated by the updated code σ^code is consistent with the scatter of the samples σ^scat. This shows that the error estimates produced by the updated code are reasonable.

Fig. 2. Fitting parameters and their errors for the case ℓ = 492 (kpix = 21). Each panel has three parts. In the upper part of each panel, squares indicate the input parameters, $q_{i,n}^{\mathrm{input}}$. The asterisks indicate the fitting results obtained by the updated code, namely, the mean of Monte-Carlo samples, $q_{i,n}^{\text{output}}$. In the middle part of each panel, the deviation of the mean fitting results from the input, $\delta q_{i,n}=q_{i,n}^{\mathrm{output}}-q_{i,n}^{\mathrm{input}}$, are illustrated by the circles and the expected scatter of the mean computed from the square root of variance σscat of the 500 Monte Carlo samples divided by the square root of the sample number, Nsample (here it is 500) are shown as the error bars. The dashed horizontal lines are at δq = 0. In the lower part of each panel the scatter of the 500 samples, σscat and the scaled error estimated by the updated code, σcode, are depicted by the thick gray lines and the thin lines with short horizontal bars on the edge. Errors estimated by the updated code are scaled by σ(α = 1.5, kpix, npix) = 0.569 as described in Appendix B.2. The error bars on the middle panel of δq1, n are tiny at this scale; the underestimation is relatively large. However, the errors of δq1, n are ${\sigma }^{\mathrm{scat}}/\sqrt{N_{\mathrm{sample}}}$ and, therefore, available from σscat in the lower panel and Nsample = 500; they are ∼0.02 at most.

We carried out the analogous Monte-Carlo simulations for the cases of twice larger k and two-thirds k. Although the number of peaks to be fitted is not the same in these cases as the number of prominent peaks is smaller for larger k, the behavior of the fitting results and their error estimates are similar to those in the case shown in Fig. 2 and Table 1.

3.2. Bias of the flow estimates and correlations between the fitting parameters

To measure the bias of the flow estimates and the correlations between the fitting parameters, we make models with a simple flow: isotropic models, except for a specific flow u_x for the n = 3 ridge.

The parameters for the model are identical to those from Sect. 3.1, except for the parameters related to the azimuthal angle θ; namely, the peak parameters q_3, n, q_4, n, q_5, n, and q_6, n (for all n) and the background parameters q_2, b and q_3, b are all set to zero. The flow u_x, n (q_3, n) is non-zero for a single n, n = 3. The limit spectrum is constructed with Eq. (2) and these input parameters, and the realizations of the power spectrum are generated in the same way as those in Sect. 3.1.

In this subsection we compare results from the updated code with results from the original code and also a modified version of the original code. While the original code fits a Lorentzian model to the square root of the power, the modified original code fits a Lorentzian model to the power itself. Although the original and modified original codes assumes the five-parameter background model (Eq. (A.3)), here we created the limit spectrum using the four-parameter background model (Eq. (4)), which we use in our updated code. As described in Sect. 2.3.2, under the current conditions, our updated four-parameter background model can approximate well the original five-parameter background model and the original and the modified original codes are applicable, although we need to keep in mind that the meaning of the background parameters in the results from different codes are different. Comparable tests using input parameters obtained by the original code give similar correlation coefficients maps.

3.2.1. Bias of the flow estimate

Figure 3 shows the fitting results for u_{x, n = 3} at ℓ = 492 (k_pix = 21) in the form of normalized histograms of 500 Monte–Carlo simulations. The input models have q_3, n = u_x, n = 0, 80, …400 m s⁻¹ for the n = 3 ridge. Before plotting, we removed outliers from the output fit parameters. We removed outliers as described in Sect. 3.1. After the outlier removal procedure, the sample numbers for each methods are 394–466 depending on the value of u_{x, n = 3} (79–93%, modified original), 436–488 (87–98%, original), and 499–500 (100%, updated). The updated code provides almost no outliers, while there are some outliers in the fitting results by the modified original and original codes. The number of outliers depends on u_x, n, although there is no clear trend with u_x, n; for example, we cannot say that the stronger flow produces more outliers. We did not further investigate the details of the outliers in the fitting results by the original and modified original codes.

Fig. 3. Fitting results of ux, n = 3 at ℓ = 492 (kpix = 21) in the form of normalized histograms by the original code which fit the model to the square root of power (black), the modified original code which fit the power (blue), and the updated code (red). Horizontal lines indicate the average of 500 Monte–Carlo realizations given by the three codes (in the same colors as the histograms, dash-dotted black, dashed blue, and solid red). The short vertical lines on the horizontal lines indicates the error estimated by the codes (thin lines with short horizontal bars) and the standard deviations of the 500 fitting results (thick lines) in the same color as the histograms centered at the means (horizontal lines) by the three codes. We note that the fitting results with overly large deviation are omitted. See text for details.

Figure 3 shows that the original fitting code underestimates the input flow by about 3%. This trend is consistent with what was reported by Greer et al. (2014). The modified original and updated codes produce less-biased flow estimates. Figure 3 also compares the errors estimated by the code and the scatter in the Monte Carlo simulation. The errors from the updated code are consistent with the scatter of the Monte-Carlo results, while the original code overestimates the errors.

3.2.2. Correlations between the fitting parameters

Figure 4 shows the correlation coefficients between the fitting parameters of 500 Monte Carlo samples. For this computation, outliers were removed as described in Sect. 3.1. The original and modified original codes both produce stronger correlations in some of the output parameters than the updated code. In particular, for the original code, these parameters show the strongest correlations: the amplitude A_n ($q_{1,n}^{\prime }$) and the width Γ_n′ ($q_{2,n^{\prime }}^{\prime }$) at peaks |n − n′| ≤ 1, the roll-off frequency ν_bg ($q_{1,\text{BG}}^{\prime }$) and the index b ($q_{2,\text{BG}}^{\prime }$) of the background, and the index b ($q_{2,\text{BG}}^{\prime }$) of the background and A_n ($q_{1,n}^{\prime }$) or Γ_n ($q_{2,n}^{\prime }$) of higher-n peaks.

Fig. 4. Correlation coefficients between the parameters (500 Monte-Carlo realizations) for the isotropic model (ux, n = 0 for all n) and for the isotropic model except with ux = 400 m s−1 for the n = 3 peak (on the next page) at ℓ = 492 (kpix = 21). Each box indicates the parameters for n-th peak (n = 0, 1, …5), qi, n (i = 0, …6, from left to right and from bottom to top in each box) and the background parameters (BG). As we defined in Sect. 2q0, n = νn is the frequency of the n-th peak, q1, n = An is the amplitude, q2, n = Γn is the width, q4,n) = u = (ux,n, uy,n) is the horizontal velocity, q5, n = fc, n and q6, n = fs, n are parameters to handle anisotropy. Background parameters are q0, BG = B0 is the amplitude, q1, BG = b is the power-law index, q2, BG = fc, bg and q3, BG = fs, bg are the parameters to handle the anisotropy of the background. Color scale is shown by the color bar on the right side. We note that the last columns and rows (i = 4) of BG on the top panels are empty because the number of the background parameters in the model in the updated code is four, instead of five.

Fig. 4. continued.

The modified original code and the updated code did not show such strong biases, except for the width and amplitude of each peak, along with the background index and some weaker peaks (smaller and larger n). In terms of correlations between the parameters, the updated code shows an overall improvement in comparison with the original, but the correlation coefficients between u_x and u_y on the same peaks or on the peaks next to each other are ∼0.1 at most for any u_x and any fitting results by all three codes shown in Fig. 4.

3.3. Summary of the Monte-Carlo test

The Monte-Carlo tests show that the updated code is able to reasonably recover the parameters, apart from the amplitude, that are used to generate the input power spectra. Other than the amplitudes, we were not able to measure a statistically meaningful bias in the fit results using 500 Monte-Carlo simulations. The underestimation of the amplitude is the result of taking the logarithm of the rebinned power. Appendix D discusses this issue in detail. Since current ring-diagram flow inversions are carried out using the mode shifts, (q_3, n, q_4, n) = (u_x, n, u_y, n), and the other fitting parameters including mode amplitudes are not used, we believe underestimated amplitudes will not substantially impact on further analysis.

4. Discussion

4.1. Further rebinning in azimuth

We explored the idea of rebinning the data further in θ. If we rebin further, the PDF of the resulting power spectrum is closer to Gaussian. Additional rebinning has the additional benefit of reducing the underestimation of the amplitude. Another benefit is that if we can reduce the number of pixels involved in the fit without a significant decrease of the fitting quality, it will reduce calculation costs.

Figure 5 shows the fitting results corresponding to the four-time further rebinning. Specifically, we rebin the 64-pixel azimuth grid into 16 pixels. In this case, valid sample numbers for the three codes are 351–429 (70–86%, modified original), 406–487 (81–97%, original), and 500 (100%, updated). While the updated code had no outliers, the numbers of outliers in the results by the original and modified original codes increase compared to the case without further rebinning (see Sect. 3.2). This also confirms the stability of the updated code.

Fig. 5. Fitting results of ux, n = 3 at ℓ = 492 (kpix = 21), similar to Fig. 3, but for the further rebinned data. In this case, the number of azimuthal pixels, npix is reduced to 16, instead of the original 64. Legends and colors are similar to those of Fig. 3. Scale factors for the error estimates by the updated code are adjusted accordingly.

This figure tells us that our updated code is better than the original for this further rebinned case as well, in terms of a more reasonable error estimate and a smaller bias. Moreover, using our updated code, we can rebin further up to 16 pixels at this k (k_pix  =  21) without a significant increase in the noise or the underestimate of the parameters or outliers in comparison with the original 64-pixel case. In the further rebinned case, the correlations between parameters show a trend that is essentially similar to that of Fig. 4. The only exception is in the fitting results obtained by the modified original code; they show less correlation even in the case of u_{x, n = 3}  =  400 m s⁻¹ in the further rebinned case.

In the original code, there is no scaling factor related to the bin number and the error estimate by the code is needed to be scaled properly. For these 64-pixel and 16-pixel cases, the scatters are not significantly changed and it is the case as well for the error estimated with the proper scaling. For these plots, again we note that we need a sufficient number of pixels for a good fitting; in this case 8 or 4 bins were too small, because the functional form of the model has terms that vary as cos θ and sin θ but also as cos2θ and sin2θ as parameters.

4.2. Future work

The modifications described in this work are mainly corrections of problems in the original code. The exception to this is the change of the algorithm from the maximum likelihood method based on the chi-square distribution function to the one based on the normal distribution function, namely the least-squares method, and the fitting not to the square-root of power but to the logarithm of the power. There are several potential improvements of the analysis. While implementing such further alterations are beyond the scope of this paper, we will briefly discuss some potential future improvements.

One of the open issues is the remapping and rebinning. Currently the power is remapped from Cartesian to polar coordinates and then smoothed in the azimuthal direction. However, it is possible to do the analysis in the original Cartesian system as it is done in the HMI pipeline fitsc module (Bogart et al. 2011a), which fits in slices at constant temporal frequency. The fitting approach shown here would only need to be slightly modified to be carried out in a region with $\sqrt{k_x^2+k_y^2}$ near k. We expect that the main issue would be extending the model to allow for small variations in k. As most parameters are presumably smooth in k, we speculate the fitting small range in the k rather than single k would help the stability of the fits.

In the current code, we approximate the PDF of the remapped and smoothed power spectrum as a normal distribution function. But how we rebin, including the topic mentioned in Sect. 4.1, is still open. In the method shown here, we do not rebin the model function but use the model function calculated on the same grid as the rebinned data. This is an approximation in the limit where the rebinning does not significantly change the shape of the limit spectrum and we expect that it will cause a bias in the case of extreme rebinning. Figures 3 and 5 show that the bias in the fitting results of the updated code is less than about 5 m s⁻¹ in the range of |u_x| ≤ 400 m s⁻¹ even in the further rebinned case (Fig. 5).

The choice of the model function (currently Eq. (2)) is also an open issue. As shown in Fig. 1, the current model function does not reproduce the detailed structure of the observed power spectrum. For example, the current model function does not take into account the asymmetry of the ridge shape in terms of frequency. Currently, we use the power spectrum of the tile at the disc center only but that is the simplest case and when we investigate the deep convection, we cannot avoid using the tiles on various locations on the disc. In such a case, in order to construct a model function, we need to take into account further effects, such as the effect of the center-to-limb variations (Zhao et al. 2012), the line-of-sight effect on shape of the power, or the effect of the Postel projection. Also, the effects of differential rotation on eigenfunctions are to be accounted for in future.

5. Conclusions

We identified several problems in the multi-ridge fitting code ATLAS (Greer et al. 2014) and we updated the code in response. We confirmed that flow-estimate biases and error overestimates exist in the fitting results by the original code. The biases that we found are insufficient to resolve the discrepancy presented in Hanasoge et al. (2016).

The updated code is based on a consistent model and an appropriate likelihood function. Monte Carlo tests show the fitting results and their error estimates by the updated code are reasonable and confirm the improvement of the fitting.

The work shown here is limited to the fitting part of the ring-diagram analysis codes. The next step in the ring-diagram analysis is flow inversion using the mode-shift parameters (q_3, n and q_4, n in Eq. (2)). Unlike the HMI ring-diagram pipeline in which the inversion is done at each tile, Greer et al. (2015) used a 3-D inversion using multiple tiles. This unique step should be the focus of examination in future works.

Acknowledgments

We thank Rick Bogart for useful discussions. The German Data Center for SDO (GDC-SDO), funded by the German Aerospace Center (DLR), provided the IT infrastructure to process the data. The HMI data used are courtesy of NASA/SDO and the HMI science teams. BWH acknowledges NASA support through grants 80NSSC17K0008, 80NSSC18K1125, and 80NSSC19K0267. We acknowledge partial support from ERC Synergy grant WHOLE SUN 810218.

References

Appendix A: Original ATLAS code

A.1. Fitting method of the original ATLAS code

In the original ATLAS code (Greer et al. 2014, 2015), the processing steps after computing the power spectrum for each tracked tile are as follows:

1. at each frequency ν, remap the log of the square root of the power spectrum $ln\sqrt{P(\mathit{k},\nu )}$ from the Cartesian (k_x, k_y) to polar (k, θ) coordinates,

2. compute $expln\sqrt{P(k,\theta ,\nu )}$ and carry out Fourier interpolation using a box-car low-pass filter to smooth in the azimuthal (θ) direction. The number of the grid points in azimuth θ is reduced from n_pix = 256 to n_pix = 64,

3. fit a Lorentzian model (see Appendix A.2) to the smoothed square-root power at each k. The fitting is done based on the maximum likelihood method assuming that the PDF of the power is the chi-square distribution with two degrees of freedom (see more details in Appendix B.1), and

4. estimate error of the fitting parameters separately from the fitting using the Fisher information matrix based on the final fitting results.

A.2. Model for the power spectrum in the original ATLAS code

In the original ATLAS code (Greer et al. 2014, 2015), the power at each k is modeled with the sum of Lorentzians with seven parameters ($q_{i,n}^{\prime },i=0,\dots 6$) for each ridge n = 0, …n_r − 1, where n_r is the number of ridges, and a background model with five parameters ($q_{i,\text{BG}}^{\prime },i=0,\dots 4$) at each k. The total number of fit parameters at each k is 7n_r + 5.

The model power spectrum, at a single k, is:

$$\begin{array}{cc}& P^{\prime }(\theta ,\nu ;{\mathit{q}}_{\mathrm{BG}}^{\prime })=\hfill \\ \hfill & {\displaystyle \sum _{n=0}^{n_{\mathrm{r}}-1}\frac{q_{1,n}^{\prime }(q_{2,n}^{\prime }/2)F^{\prime }(\theta ;q_{5,n}^{\prime },q_{6,n}^{\prime })}{{[\nu -q_{0,n}^{\prime }+k(q_{3,n}^{\prime }cos\theta +q_{4,n}^{\prime }sin\theta )/(2\pi )]}^2+{(q_{2,n}^{\prime }/2)}^2}}\hfill \\ \hfill & +B^{\prime }(\theta ,\nu ;{\mathit{q}}_{\mathrm{BG}}^{\prime }),\hfill \end{array}$$(A.1)

where

$$\begin{array}{c}\hfill F^{\prime }(\theta ;p_0,p_1)=1+p_{0}cos[2(\theta -p_1)],\end{array}$$(A.2)

and $q_{0,n}^{\prime }={\nu }_n$ is the frequency of the nth peak, $q_{1,n}^{\prime }$ = A_n is the amplitude, $q_{2,n}^{\prime }$ = Γ_n is the width, $(q_{3,n}^{\prime },q_{4,n}^{\prime })=\mathit{u}=(u_{x,n},u_{y,n})$ is the horizontal velocity, $q_{5,n}^{\prime }=f_n$ and $q_{6,n}^{\prime }={\theta }_n$ are anisotropy terms in Eqs. (5) and (6) in Greer et al. (2014). The background is modeled at each k as

$$\begin{array}{c}\hfill B^{\prime }(\theta ,\nu ;{\mathit{q}}_{\mathrm{BG}}^{\prime })=\frac{q_{0,\mathrm{BG}}^{\prime }}{1+{(\nu /q_{1,\mathrm{BG}}^{\prime })}^{q_{2,\mathrm{BG}}^{\prime }}}{F}^{\prime }(\theta ;q_{3,\mathrm{BG}}^{\prime },q_{4,\mathrm{BG}}^{\prime }),\end{array}$$(A.3)

where $q_{0,\text{BG}}^{\prime }=B_0$ is an amplitude, $q_{1,\text{BG}}^{\prime }$ = ν_bg is a roll-off frequency, $q_{2,\text{BG}}^{\prime }=b$ is the power-law index, $q_{3,\text{BG}}^{\prime }=f_{\text{bg}}$ is the amplitude of the anisotropy term, and $q_{4,\text{BG}}^{\prime }={\theta }_{\text{bg}}$ is the phase of the anisotropy term in Eq. (7) in Greer et al. (2014).

A.3. Problems in the original ATLAS code

First, the Lorentzian model for the power in the original code was fit to the square root of the observed power. This is inconsistent with what is stated in Greer et al. (2014, 2015) and it is also an inconsistency between the model and the observable. Therefore, we made them consistent in the updated code (see Sect. 2.2).

Second, the cost function to be minimized in the original code was not a good approximation. The cost function to be minimized in the maximum likelihood method based on the chi-square-distribution with two degrees of freedom as PDF is

$$\begin{array}{c}\hfill -ln\mathcal{L}(\mathit{q})={\displaystyle \sum _i\{ln\left(\frac{P_i(\mathit{q})}{O_i}\right)+\left(\frac{O_i}{P_i(\mathit{q})}\right)\},}\end{array}$$(A.4)

where O_i is the observed spectrum and P_i is the model, as given in Sect. A.2. The sum is taken over all grid points (θ_i, ν_i) in the fitting range.

Instead, the original code minimizes

$$\begin{array}{c}\hfill {\displaystyle \sum _i{\mathcal{B}}_i^2\equiv \sum _i{\{ln\left(\frac{P_i(\mathit{q})}{O_i}\right)+\left(\frac{O_i}{P_i(\mathit{q})}\right)\}}^2}\end{array}$$(A.5)

using the mpfit.c code¹. This is conceptually inconsistent with the description in Greer (2015), although the likelihood function itself is not explicitly stated there. While it can be shown that the results of minimizing of the exact cost function (Eq. (B.3)) and the wrong one (Eq. (A.5)) are identical in the limit of linear perturbations, there is no reasonable computational or physical reason to take the extra square in the calculations. Moreover, it is not useful for the error estimates. Therefore, we decided to correct this issue; Sect. 2.2 describes our changes.

Third, in the original ATLAS code, error estimates are obtained by an independent calculation of the Fischer information using the final fitting results and using the chi-square distribution with two degrees of freedom as the likelihood function. It is not consistent to compute the error estimates with a likelihood function that is different than what was used in the fitting itself. In Sect. 2.2, we also describe a consistent approach to the computation of error estimates in our updated code.

Fourth, several parameters are unstable. Therefore, we have changed the parameterizations as we discussed in Sects. 2.3.1 and 2.3.2.

Appendix B: Probability distribution function (PDF) and maximum likelihood method based on the PDF

B.1. PDF of the raw power spectrum

Woodard (1984) demonstrated that the PDF of a single observed power spectrum divided by the expectation value of the spectrum is the chi-square distribution with two degrees of freedom. On this basis, Duvall & Harvey (1986) and Anderson et al. (1990) introduced the probability density at a given grid point in the Fourier space (k_x, k_y, ν)_i

$$\begin{array}{c}\hfill \mathcal{P}(O_i)=\frac{1}{P_i(\mathit{q})}exp(-\frac{O_i}{P_i(\mathit{q})}),\end{array}$$(B.1)

where O_i is the observed spectrum, P_i is the model for the limit spectrum, and q is the model parameters. The joint probability density for the experimental outcome at horizontal independent wavevectors and frequencies is given by

$$\begin{array}{c}\hfill \mathcal{L}=\prod _i\mathcal{P}(O_i)=exp[-{\displaystyle \sum _i(ln{P}_i(\mathit{q})+\frac{O_i}{P_i(\mathit{q})})}].\end{array}$$(B.2)

This is the likelihood function, and the model parameters which maximize this ℒ are the targets in the maximum likelihood method. In practice −lnℒ is minimized to use standard minimization procedures. To make the computation simpler, lnO_i is subtracted from −lnℒ. The minimization of −lnℒ − lnO_i in terms of q is identical to that of −lnℒ, because O_i is not a function of the model parameter q. In conclusion,

$$\begin{array}{c}\hfill -ln\mathcal{L}-ln{O}_i={\displaystyle \sum _i(ln\frac{P_i(\mathit{q})}{O_i}+\frac{O_i}{P_i(\mathit{q})})}\end{array}$$(B.3)

is minimized.

B.2. PDF of logarithm of averaged power

The PDF of a well-averaged power spectrum is also a chi-square distribution but with many more than the two degrees of freedom of the original and given the central limit theorem, it can be approximated by a normal distribution. To carry out the least-squares fitting based on the normal distribution function and estimate errors of the fitting results, we need the variance of the normal distribution function. We therefore take advantage of the fact that the logarithm of the averaged spectrum obeys a normal distribution function with a constant variance. In Appendix B.3, we show that if we have a spectrum obtained by averaging over N spectra, each normally distributed with the same expectation value M, the logarithm of averaged spectra, y, obeys the normal distribution function $N(\ln M,{\sigma }_N^2)$:

$$\begin{array}{c}\hfill f(y)=\frac{1}{\sqrt{2\pi }{\sigma }_N}exp[-\frac{1}{2}{\left(\frac{y-lnM}{{\sigma }_N}\right)}^2],\end{array}$$(B.4)

where ${\sigma }_N=1/\sqrt{N}$, therefore, σ_N is a constant over y. We note that the function f(y) is linearized around the mean to derive Eq. (B.4); see Appendix B.3 for details.

In the present case, we fit one k at a time. To do this, the spectra are, at each frequency, linearly interpolated in k_x and k_y to a circle with radius k_pix, where k_pix is k in units of bins, and smoothed to n_pix azimuths. On this basis, the question is if we can still use a least-squares fit with a diagonal covariance and a single σ_N and, if so, which value of σ_N (or equivalently N) in Eq. (B.4) should be used. In particular, it needs to be considered that more than 2πk_pix are used for the interpolation and that the averaged values will be correlated.

In the limit of averaging over the entire circle, namely n_pix = 1, and in the case with k_pix ≫ 1, it might be reasonable to assume that N = 2πk_pixα, with α accounting for the fact that the averaging is essentially over an annulus around k_pix. Indeed, a Monte-Carlo test shows that α ≈ 1.5 gives a very good estimate of the error on the average, which will be shown in Fig. B.1 later in this subsection.

Fig. B.1. Standard deviation of the logarithm of chi-square distributed white noise rebinned on the annulus with the radius of kpix pixels, σscat, is plotted against the square-root of the pixel number on the annulus after rebinning, $\sqrt{n_{\mathrm{pix}}}$, at three kpix shown with different symbols in different colors. In the case of npix ≪ 2πkpix, σscat = σα(α = 1.5), which is shown with the dotted lines, is a valid approximation, while in the limit of npix ≳ 2πkpix, σscat deviates from σα and approaches a constant of 2/3.

For n_pix ≫ 1 the variances, as well as the off-diagonal elements of the covariance matrix of the interpolated datapoints will, in general, depend on azimuth. Assuming that the fitted function may be linearized in the fitted parameters, a linear fit implies that the fitted parameters are given by a linear combination of the observed values. Assuming that the functions are smooth (as in the present case where they are low-order harmonic functions), the coefficients in the linear combination are also smooth. From this it follows that if the variations in the properties of the covariance matrix average on the scale of the variations in the fitted functions, then the same scaling factor may be used and that

$$\begin{array}{c}\hfill N\sim \frac{2\pi k_{\mathrm{pix}}\alpha }{n_{\mathrm{pix}}},\end{array}$$(B.5)

which a Monte-Carlo test again confirms.

To validate the arguments above, a simple Monte-Carlo test is carried out to measure the scatter (standard deviation, σ^scat) of the logarithm of averaged white noise. White-noise fields with the chi-square distribution with two degree of freedom in a three-dimensional Cartesian Fourier space (k_x, k_y, ν) = (384, 384, 1152) [pix] are remapped at several annuli with specific radii k_pix in the same way as the data in the updated analysis code. At first, there are 256 pixels on each annulus after remapping, we then rebin them into n_pix = 4, 8, 16, 32, 64, and 128 pixels. The standard deviation of the logarithm of the rebinned white noise, σ^scat, are computed and plotted against the square root of the number of datapoints after rebin on the annuli, n_pix = 14, 21, and 42 in Fig. B.1. As we mentioned above, in the case of n_pix ≪ 2πk_pix, σ^scat ≃ σ_α(α = 1.5), where ${\sigma }_{\alpha }(\alpha )\equiv \sqrt{n_{\mathrm{pix}}/(2\pi k_{\mathrm{pix}}\alpha )}$. In the case of n_pix ≳ 2πk_pix the deviation of σ^scat from σ_α is not negligible, and σ^scat approaches a constant 2/3; this constant can be derived from the bi-linear interpolations of chi-square distributed random variables in two dimensions. Also we note that we have not included the apodization effect in the data to create the noise field in this test calculation. Without apodization σ^scat might be overestimated; in any case, this does not affect the parameter estimates but only the error estimates.

In our updated code, σ_N is set to 1 (see the cost function Eq. (1)), as we mentioned in Sect. 2.2, and after the fitting is done the error provided by the code is scaled by σ_α. Although calling the minimization code with an incorrect error estimate can lead to worse convergence, the difference in this case should be negligible. Implementing σ_α in the code is a task for the future.

B.3. Derivation of PDF of the logarithm of averaged power

In this section we briefly summarize how to derive the PDF of the logarithm of the power.

Suppose x is the average of N spectra whose averages are M. If N is large enough, then the PDF of this spectrum is the normal distribution

$$\begin{array}{c}\hfill f_x(x)=\frac{1}{\sqrt{2\pi }\sigma }exp[-\frac{1}{2}{\left(\frac{x-M}{\sigma }\right)}^2]\end{array}$$(B.6)

and $\sigma =M/\sqrt{N}$. Here we define

$$\begin{array}{c}\hfill y=g(x)\equiv lnx\end{array}$$(B.7)

and when x is close to M,

$$\begin{array}{c}\hfill y=ln\left[M(1+\frac{x-M}{M})\right]\simeq lnM+\frac{x-M}{M}·\end{array}$$(B.8)

Therefore, the inverse function of g is given by

$$\begin{array}{c}\hfill g^{-1}(y)=x=(y-lnM)M+M\end{array}$$(B.9)

and

$$\begin{array}{c}\hfill \frac{\mathrm{d}}{\mathrm{d}y}(g^{-1}(y))=M.\end{array}$$(B.10)

Hence, the PDF of y = lnx is given by

$$\begin{array}{cc}\hfill f_y(y)& =\left|\frac{\mathrm{d}}{\mathrm{d}y}(g^{-1}(y))\right|f_x(g^{-1}(y))\hfill \\ \hfill & =\frac{1}{\sqrt{2\pi }{\sigma }_y}exp[-\frac{1}{2}{\left(\frac{y-lnM}{{\sigma }_y}\right)}^2],\hfill \end{array}$$(B.11)

where ${\sigma }_y=\sigma /M=1/\sqrt{N}$. Therefore, the logarithm of the averaged power has a constant error, which is independent of the original σ and depends only on the number of the averaged spectra, N. For comparison, approximated PDF and the distribution functions of the Monte-Carlo samples are given in Appendix B.4.

B.4. Comparison of the distribution functions with Gaussian and chi-square distribution functions

In this subsection we show that Gaussian function is a good approximation for the distribution functions of the logarithm of the synthesized power. We also discuss the degrees of freedom of the chi-square distribution functions of the unwrapped and smoothed power here.

Figure B.2 shows the distribution function of the power and the logarithm of the power at k_pix = 21 (ℓ = 492) from the Monte Carlo simulations. The Gaussian function centered at zero whose width is the standard deviation of the samples are overplotted in the figure. These plots show that the Gaussian function is a good approximation for the distribution function of the logarithm of the unwrapped and smoothed power near the mean value, though the tails of the distribution function are longer than those of the Gaussian.

Fig. B.2. Normalized frequency distribution functions of the normalized synthetic power at kpix = 21 (ℓ = 492, black solid lines): top panels are for the subsampled power and the second panels are the unwrapped power without further smoothing, while the lower panels are for the rebinned powers. Left panels: (P − ⟨P⟩)/⟨P⟩, and right panels: lnP − ln⟨P⟩, where ⟨⟩ indicates the expectation value (model). The input model of the power is isotropized (qi, n = 0 for 3 ≤ i ≤ 6 for all n and qi, b = 0 for i = 2, 3) but otherwise it is the same as the model with the parameters given in Table C.2. For this plot, we use 50 realizations for the top two rows and 500 realizations for the rebinned data (lower three rows). The black dashed lines on the lower four panels are the Gaussian function centered at zero whose width is the standard deviation of the samples (shown on the panels). This validates our approximation of the distribution function of the logarithm of the power by the Gaussian function. The thick red dashed curves on the left panels are the chi-square distribution function of various degrees of freedom. The degree of freedom (d.o.f.) for each distribution is determined from the standard deviation of each set of normalized power. See the text for details.

In the left panels of Fig. B.2, the distribution function of the power is also compared with the chi-square distribution functions with various degrees of freedom: the chi-square distribution functions shown here are

$$\begin{array}{c}\hfill f(X;d)={\left(\frac{d}{2}\right)}^{d/2}\frac{{(X+1)}^{d/2-1}}{\mathrm{\Gamma }(d/2)}exp(-\frac{d}{2}(X+1)),\end{array}$$(B.12)

where X = (P − ⟨P⟩)/⟨P⟩, ⟨⟩ indicates the expectation value, and d is the number of degrees of freedom (d.o.f.). Since in this case ⟨X⟩ = 0 and the variance of X is 2/d, we determine d.o.f. from the variance and show the corresponding chi-square distribution function in the figure. As the original synthesized power is created from the model and the two-d.o.f. chi-square distributed random numbers, subsampled power is purely two d.o.f. The unwrapped power has a four-d.o.f. chi-square distribution; this indicates that unwrapped power is effectively the subsampled power rebinned over two pixels. Because in this case the radius of annulus is k_pix = 21 and hence roughly 2πk_pix independent pixels are on the annulus, the unwrapped power with 256 pixels on the annulus is roughly twice oversampled. Therefore, the d.o.f. of the power unwrapped and rebinned over two pixels in azimuth (n_pix = 128, third panel) is five. This is almost unchanged from the original (second panel, d.o.f. = 4). This is also the case when we further rebin the data over four pixels (fourth panel); d.o.f is only eight and this is smaller than 16, which we might have expected for the four-pixel rebinned four-d.o.f. distributed independent variables (i.e., unwrapped power). But we should also note here that 12 is the d.o.f. which we expect on the basis of the effective rebinning numbers we discussed in Appendix B.2 (Eq. (B.5)), N ∼ 3, where α = 1.5. It is also larger than the dof of the distribution function (8); this is partly because the approximation (Eq. (B.5)) is not sufficiently good for n_pix = 64 case as was shown in Fig. B.1. When we see the further rebinned case with n_pix = 16, namely the unwrapped power averaged over 16 pixels (the lower-most panel), the d.o.f. is 25 and close to what we expect from Eq. (B.5), N ∼ 6: in this case the d.o.f. is 6 × 4 = 24. For this case the approximation by Eq. (B.5) with α = 1.5 is good according to Fig. B.1. This closeness of the distribution function of the remapped and smoothed power to the chi-square distribution with some degrees of freedom suggest that one might expect improvement to the fitting method by using the maximum likelihood method based on the chi-square distribution with specific degrees of freedom.

Appendix C: Input parameters for smaller and larger wavenumbers

The parameters at several ℓ obtained from the fitting of the observed spectra averaged over one Carrington rotation 2211 for the disc center tile from SDO/HMI pipeline are shown in Tables C.1 (ℓ  =  328), C.2 (ℓ  =  492), and C.3 (ℓ  =  984). These average power spectra were obtained from DRMS specification hmi.rdvavgpspec_fd15[2211][0][0] and the fitting was done using the updated code. These parameters are used as input parameters to generate Monte-Carlo synthetic data in Sect. 3.

Appendix D: Underestimated amplitude – Logarithm of average and average of logarithm

Figure 2 shows that the amplitudes (q_1, n for all n and q_0, b in Eq. (2)) measured by the updated code are underestimated. This is the result of taking the logarithm of the rebinned power and the effect can be reproduced in the following simple test.

To see why this is the case, we define x_i (i = 0, 1, …n − 1) as a series of random variables whose distribution function is the chi-square distribution with two degrees of freedom. The expectation values and the variance of x are ⟨x⟩ = 2 and ${\sigma }_x^2=2^2$, respectively, and y_j(j = 0, 1, …n/a − 1) is smoothed x_i over each a-pixel range:

$$\begin{array}{c}\hfill y_j=\frac{1}{a}{\displaystyle \sum _{i=ja}^{(j+1)a-1}{x}_i.}\end{array}$$(D.1)

In this case ⟨y⟩=⟨x⟩ = 2, and ln⟨y⟩=ln(2), but ⟨ln(y)⟩ ≤ ln⟨y⟩. The more the array is rebinned (the larger a becomes), the less the scatter of y becomes: ⟨ln(y)⟩ approaches ln⟨y⟩.

This trend is alleviated if the spectra are more strongly averaged. Figure D.1 shows a simple test calculation to illustrate how the average of logarithm is reduced from the logarithm of the average; the ratio of ⟨ln(y)⟩ to ln⟨y⟩ is plotted against the pixel number ratio of the rebinned data to the original. In this simple test calculation, we make 500 realizations of sets of x_i (i = 0, 1, …n − 1, where n = 256), calculate y with a = 1, 2, 4, 8, …256, and compute how much smaller the expectation value of the logarithm of y than the logarithm of the expectation value is. The linear regression of the data points are also shown in the plot. For comparison the amplitude fitting results shown in Fig. D.2 are replotted here with squares. This explains how the logarithm of amplitude is reduced. In the case of no rebinning (X = n_rebin/n_original = 1) or only 2-pixel rebin (X = 1/2) the fitting results are deviated significantly from this simple test calculation, but this comes as no surprise because in these cases, the assumption of the well-averaged power as data is not appropriate.

Fig. D.1. exp⟨ln(y)⟩/⟨y⟩, where y is the a-element average of chi-square distributed random variables, as a function of the pixel number ratio of the rebinned data to the original (asterisks). The average ratio of the amplitude fitting results to the input shown in Fig. D.2 is depicted by squares. The dashed line and the equation on the panel are the linear regression of the data points.

Figure D.2 shows the dependence of the amplitude measurements on the amount of rebinning. The input model is identical to the one used in the calculation in Sect. 3.1. The default rebinning is from 256 pixels to 64 pixels and in this case, the output amplitude is about 88% of the input amplitude but with further rebinning to 16 pixels, it increases up to 96%.

Fig. D.2. Dependence of the peak amplitude parameters at ℓ = 492 (kpix = 21) measured by the updated code on the rebinning strength. Left panel: input (black dotted lines with squares) and measurements with three different rebinning with error bars (scatter of 500 realizations): no rebinning (256 pixels on the azimuthal grid, blue), original rebinning (64 pixels, green), and extra rebinning (16 pixels, red). Right panel: ratio of the amplitude to the input. The numbers are the grid number and the average ratio of six peaks. The input parameters are identical to those in Fig. 2.

All Tables

All Figures

Fig. 1. Slices through the input limit spectrum (black) and observational spectrum averaged over one Carrington rotation (gray), which was used to obtain the input parameters, at θ = 0 and a few different values of k. Vertical dashed lines indicate the fitting ranges and the models are plotted only within these ranges. Lower limits for all k are fixed at 0.4 mHz, while the upper limits depend on the initial guess for each k, and they are the highest peak frequencies plus the widths of the peak. For the cases of ℓ = 328 and 492 (kpix = 14 and 21, respectively), the peaks with the radial order, n, from 0 to 7 are used in the fitting, while for the case of ℓ = 984 (kpix = 42), 0 ≤ n ≤ 3.

In the text

Fig. 2. Fitting parameters and their errors for the case ℓ = 492 (kpix = 21). Each panel has three parts. In the upper part of each panel, squares indicate the input parameters, $q_{i,n}^{\mathrm{input}}$. The asterisks indicate the fitting results obtained by the updated code, namely, the mean of Monte-Carlo samples, $q_{i,n}^{\text{output}}$. In the middle part of each panel, the deviation of the mean fitting results from the input, $\delta q_{i,n}=q_{i,n}^{\mathrm{output}}-q_{i,n}^{\mathrm{input}}$, are illustrated by the circles and the expected scatter of the mean computed from the square root of variance σscat of the 500 Monte Carlo samples divided by the square root of the sample number, Nsample (here it is 500) are shown as the error bars. The dashed horizontal lines are at δq = 0. In the lower part of each panel the scatter of the 500 samples, σscat and the scaled error estimated by the updated code, σcode, are depicted by the thick gray lines and the thin lines with short horizontal bars on the edge. Errors estimated by the updated code are scaled by σ(α = 1.5, kpix, npix) = 0.569 as described in Appendix B.2. The error bars on the middle panel of δq1, n are tiny at this scale; the underestimation is relatively large. However, the errors of δq1, n are ${\sigma }^{\mathrm{scat}}/\sqrt{N_{\mathrm{sample}}}$ and, therefore, available from σscat in the lower panel and Nsample = 500; they are ∼0.02 at most.

In the text

Fig. 3. Fitting results of ux, n = 3 at ℓ = 492 (kpix = 21) in the form of normalized histograms by the original code which fit the model to the square root of power (black), the modified original code which fit the power (blue), and the updated code (red). Horizontal lines indicate the average of 500 Monte–Carlo realizations given by the three codes (in the same colors as the histograms, dash-dotted black, dashed blue, and solid red). The short vertical lines on the horizontal lines indicates the error estimated by the codes (thin lines with short horizontal bars) and the standard deviations of the 500 fitting results (thick lines) in the same color as the histograms centered at the means (horizontal lines) by the three codes. We note that the fitting results with overly large deviation are omitted. See text for details.

In the text

Fig. 4. Correlation coefficients between the parameters (500 Monte-Carlo realizations) for the isotropic model (ux, n = 0 for all n) and for the isotropic model except with ux = 400 m s−1 for the n = 3 peak (on the next page) at ℓ = 492 (kpix = 21). Each box indicates the parameters for n-th peak (n = 0, 1, …5), qi, n (i = 0, …6, from left to right and from bottom to top in each box) and the background parameters (BG). As we defined in Sect. 2q0, n = νn is the frequency of the n-th peak, q1, n = An is the amplitude, q2, n = Γn is the width, q4,n) = u = (ux,n, uy,n) is the horizontal velocity, q5, n = fc, n and q6, n = fs, n are parameters to handle anisotropy. Background parameters are q0, BG = B0 is the amplitude, q1, BG = b is the power-law index, q2, BG = fc, bg and q3, BG = fs, bg are the parameters to handle the anisotropy of the background. Color scale is shown by the color bar on the right side. We note that the last columns and rows (i = 4) of BG on the top panels are empty because the number of the background parameters in the model in the updated code is four, instead of five.

In the text

	Fig. 4. continued.
In the text

	Fig. 5. Fitting results of ux, n = 3 at ℓ = 492 (kpix = 21), similar to Fig. 3, but for the further rebinned data. In this case, the number of azimuthal pixels, npix is reduced to 16, instead of the original 64. Legends and colors are similar to those of Fig. 3. Scale factors for the error estimates by the updated code are adjusted accordingly.
In the text

Fig. B.1. Standard deviation of the logarithm of chi-square distributed white noise rebinned on the annulus with the radius of kpix pixels, σscat, is plotted against the square-root of the pixel number on the annulus after rebinning, $\sqrt{n_{\mathrm{pix}}}$, at three kpix shown with different symbols in different colors. In the case of npix ≪ 2πkpix, σscat = σα(α = 1.5), which is shown with the dotted lines, is a valid approximation, while in the limit of npix ≳ 2πkpix, σscat deviates from σα and approaches a constant of 2/3.

In the text

Fig. B.2. Normalized frequency distribution functions of the normalized synthetic power at kpix = 21 (ℓ = 492, black solid lines): top panels are for the subsampled power and the second panels are the unwrapped power without further smoothing, while the lower panels are for the rebinned powers. Left panels: (P − ⟨P⟩)/⟨P⟩, and right panels: lnP − ln⟨P⟩, where ⟨⟩ indicates the expectation value (model). The input model of the power is isotropized (qi, n = 0 for 3 ≤ i ≤ 6 for all n and qi, b = 0 for i = 2, 3) but otherwise it is the same as the model with the parameters given in Table C.2. For this plot, we use 50 realizations for the top two rows and 500 realizations for the rebinned data (lower three rows). The black dashed lines on the lower four panels are the Gaussian function centered at zero whose width is the standard deviation of the samples (shown on the panels). This validates our approximation of the distribution function of the logarithm of the power by the Gaussian function. The thick red dashed curves on the left panels are the chi-square distribution function of various degrees of freedom. The degree of freedom (d.o.f.) for each distribution is determined from the standard deviation of each set of normalized power. See the text for details.

In the text

	Fig. D.1. exp⟨ln(y)⟩/⟨y⟩, where y is the a-element average of chi-square distributed random variables, as a function of the pixel number ratio of the rebinned data to the original (asterisks). The average ratio of the amplitude fitting results to the input shown in Fig. D.2 is depicted by squares. The dashed line and the equation on the panel are the linear regression of the data points.
In the text

Fig. D.2. Dependence of the peak amplitude parameters at ℓ = 492 (kpix = 21) measured by the updated code on the rebinning strength. Left panel: input (black dotted lines with squares) and measurements with three different rebinning with error bars (scatter of 500 realizations): no rebinning (256 pixels on the azimuthal grid, blue), original rebinning (64 pixels, green), and extra rebinning (16 pixels, red). Right panel: ratio of the amplitude to the input. The numbers are the grid number and the average ratio of six peaks. The input parameters are identical to those in Fig. 2.

In the text