An improved multi-ridge fitting method for ring-diagram helioseismic analysis

Context: There is a wide discrepancy in current estimates of the strength of convection flows in the solar interior obtained using different helioseismic methods applied to observations from SDO/HMI. The cause for these disparities is not known. Aims: As one step in the effort to resolve this discrepancy, we aim to characterize the multi-ridge fitting code for ring-diagram helioseismic analysis that is used to obtain flow estimates from local power spectra of solar oscillations. Methods: We updated the multi-ridge fitting code developed by Greer et al.(2014) to solve several problems we identified through our inspection of the code. In particular, we changed the merit function to account for the smoothing of the power spectra, model for the power spectrum, and noise estimates. We used Monte Carlo simulations to generate synthetic data and to characterize the noise and bias of the updated code by fitting these synthetic data. Results: The bias in the output fit parameters, apart from the parameter describing the amplitude of the p-mode resonances in the power spectrum, is below what can be measured from the Monte-Carlo simulations. The amplitude parameters are underestimated; this is a consequence of choosing to fit the logarithm of the averaged power. We defer fixing this problem as it is well understood and not significant for measuring flows in the solar interior. The scatter in the fit parameters from the Monte-Carlo simulations is well-modeled by the formal error estimates from the code. Conclusions: We document and demonstrate a reliable multi-ridge fitting method for ring-diagram analysis. The differences between the updated fitting results and the original results are less than one order of magnitude and therefore we suspect that the changes will not eliminate the aforementioned orders-of-magnitude discrepancy in the amplitude of convective flows in the solar interior.


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 timedistance (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 in-A&A proofs: manuscript no. Nagashima_arxiv novation: they chose a much denser tile layout than other ringdiagram 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 AT-LAS code from Greer et al. (2014Greer et al. ( , 2015 provides seven parameters for each ridge and five parameters for the background power at each wavenumber k. Greer et al. (2015) applied a threedimensional 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 abovementioned disagreement between the helioseismic measurements of subsurface convection, here we focus on the ringdiagram 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.

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.

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 kR ⊙ ≡ ℓ ∼ 500 (where R ⊙ is the solar radius), which corresponds to k pix ≡ k/h k = 21 (where h k = 3.37 × 10 −2 Mm −1 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.

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: 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 ln O 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.
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 n-th 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 where F is again given by Eq.
(3) and q 0,BG = B 0 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.2, and 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.

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 and q ′ 3,BG ) are usually much smaller than one. As a consequence, the phase (q ′ 6,n or q ′ 4,BG ) does not matter much in the fitting, and, hence, it is unstable. In the particular case of isotropic power spectra with q ′ 5,n = 0 for all n and q ′ 3,BG = 0, the phases q ′ 6,n and q ′ 4,BG are indeterminate.

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): 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,BG = 2. Appourchaux et al. (2002) suggested a generalization of the Harvey (1985) model, where q ′ 2,BG can be the range of 2 to 6. In the ATLAS fittings, however, we typically obtain the index q ′ 2,BG ∼ 1 for HMI observations. Also, the roll-off frequency obtained from the fitting of HMI observations, q ′ 1,BG , is quite low (∼ 1 µHz) and below the frequency resolution for 28.8-hour power spectra (9.7 µHz), which is the typical observation length for 15-degree tiles in HMI pipeline and 16-degree tiles in AT-LAS.
This suggests that the background is not related to either the supergranulation (τ ≈ 10 5 sec, or ν ≈ 10 µHz) or granulation time scales (τ ≈ 4 × 10 2 sec, ν ≈ 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 ν ≫ q ′ 1,BG , the original background model, Eq. (A.3), can be simplified: 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.

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.

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 A&A proofs: manuscript no. Nagashima_arxiv 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.2, and C.3, respectively, in Appendix C. 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 in Appendix C. 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.
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 (σ scat / N 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.
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.

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. 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 −1 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. 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. 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 ) and the width Γ n ′ (q ′ 2,n ′ ) at peaks |n − n ′ | ≤ 1, the roll-off frequency ν bg (q ′ 1,BG ) and the index b (q ′ 2,BG ) of the background, and the index b (q ′ 2,BG ) of the background and A n (q ′ 1,n ) or Γ n (q ′ 2,n ) of higher-n peaks. 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.

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.

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 fourtime 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.
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 −1 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 cos 2θ and sin 2θ as parameters.

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 k 2 x + k 2 y 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 −1 in the range of |u x | ≤ 400 m s −1 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.

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 ringdiagram 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 A&A proofs: manuscript no. Nagashima_arxiv the inversion is done at each tile, Greer et al. (2015) used a 3-D inversion using multiple tiles. This unique step should be be the focus of examination in future works. 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 the peak. . In the middle part of each panel, the deviation of the mean fitting results from the input, δq i,n = q output i,n − q input i,n , 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, N sample (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, k pix , n pix ) = 0.569 as described in Appendix B.2. The error bars on the middle panel of δq 1,n are tiny at this scale; the underestimation is relatively large. However, the errors of δq 1,n are σ scat / N sample and, therefore, available from σ scat in the lower panel and N sample = 500; they are ∼ 0.02 at most. A&A proofs: manuscript no. Nagashima_arxiv 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. Table 1. Background parameters in the input model for the Monte Carlo simulation and the fitting results at ℓ = 492 (k pix = 21). 500 realizations were used. Figure 2 shows the corresponding fitting results for the peak parameters. The standard deviations over the 500 realizations and the scaled errors provided by the updated code are consistent. and for the isotropic model except with u x = 400 m s −1 for the n = 3 peak (on the next page) at ℓ = 492 (k pix = 21). Each box indicates the parameters for n-th peak (n = 0, 1, . . . 5), q i,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. 2 q 0,n = ν n is the frequency of the n-th 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. Background parameters are q 0,BG = B 0 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 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.
Article number, page 9 of 18 A&A proofs: manuscript no. Nagashima_arxiv In the original ATLAS code (Greer et al. 2014(Greer et al. , 2015, the power at each k is modeled with the sum of Lorentzians with seven parameters (q ′ i,n , i = 0, . . . 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,BG , i = 0, . . . 4) at each k. The total number of fit parameters at each k is 7n r + 5.

Appendix 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. (2014Greer et al. ( , 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 chisquare-distribution with two degrees of freedom as PDF is 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 using the mpfit.c code 1 . 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 Appendix 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 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 This is the likelihood function, and the model parameters which maximize this L are the targets in the maximum likelihood method. In practice − ln L is minimized to use standard minimization procedures. To make the computation simpler, ln O i is subtracted from − ln L. The minimization of − ln L − ln O i in terms of q is identical to that of − ln L, because O i is not a function of the model parameter q. In conclusion, is minimized.
Appendix 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 leastsquares 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, σ 2 N ): (B.4) where σ N = 1/ √ 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.
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 Standard deviation of the logarithm of chi-square distributed white noise rebinned on the annulus with the radius of k pix pixels, σ scat , is plotted against the square-root of the pixel number on the annulus after rebinning, √ n pix , at three k pix shown with different symbols in different colors. In the case of n pix ≪ 2πk pix , σ scat = σ α (α = 1.5), which is shown with the dotted lines, is a valid approximation, while in the limit of n pix 2πk pix , σ scat deviates from σ α and approaches a constant of 2/3. 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 σ α (α) ≡ n pix /(2πk pix α). 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. 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 σ 2 x = 2 2 , respectively, and y j ( j = 0, 1, . . . n/a−1) is smoothed x i over each  Left panels show (P − P )/ P , and right panels show ln P − ln P , where indicates the expectation value (model). The input model of the power is isotropized (q i,n = 0 for 3 ≤ i ≤ 6 for all n and q i,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 (dof) for each distribution is determined from the standard deviation of each set of normalized power. See the text for details. a-pixel range: 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 re-sults 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. 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%. Input parameter set at ℓ = 328 (k pix = 14) obtained from the fitting the average power spectrum over one Carrington rotation to the model by the updated code. See text for more details. These parameters are used to construct the limit spectrum for Monte-Carlo simulations in Sect. 3.   Dependence of the peak amplitude parameters at ℓ = 492 (k pix = 21) measured by the updated code on the rebinning strength. Left panel shows the 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). The right panel shows the 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.