Model fitting of kink waves in the solar atmosphere: Gaussian damping and timedependence
Department of Mathematics & Information Sciences, Northumbria University, Newcastle Upon Tyne, NE1 8ST, UK
email: richard.morton@northumbria.ac.uk
Received: 31 March 2016
Accepted: 17 June 2016
Aims. Observations of the solar atmosphere have shown that magnetohydrodynamic waves are ubiquitous throughout. Improvements in instrumentation and the techniques used for measurement of the waves now enables subtleties of competing theoretical models to be compared with the observed waves behaviour. Some studies have already begun to undertake this process. However, the techniques employed for model comparison have generally been unsuitable and can lead to erroneous conclusions about the best model. The aim here is to introduce some robust statistical techniques for model comparison to the solar waves community, drawing on the experiences from other areas of astrophysics. In the process, we also aim to investigate the physics of coronal loop oscillations.
Methods. The methodology exploits leastsquares fitting to compare models to observational data. We demonstrate that the residuals between the model and observations contain significant information about the ability for the model to describe the observations, and show how they can be assessed using various statistical tests. In particular we discuss the KolmogorovSmirnoff one and two sample tests, as well as the runs test. We also highlight the importance of including any observational trend line in the modelfitting process.
Results. To demonstrate the methodology, an observation of an oscillating coronal loop undergoing standing kink motion is used. The model comparison techniques provide evidence that a Gaussian damping profile provides a better description of the observed wave attenuation than the often used exponential profile. This supports previous analysis from Pascoe et al. (2016, A&A, 585, L6). Further, we use the model comparison to provide evidence of timedependent wave properties of a kink oscillation, attributing the behaviour to the thermodynamic evolution of the local plasma.
Key words: Sun: corona / waves / magnetohydrodynamics (MHD) / Sun: oscillations
© ESO, 2016
1. Introduction
The Sun’s atmosphere is known to be replete with magnetohydrodynamic (MHD) wave phenomena and current instrumentation has enabled the measurement of the waves and their properties. In particular, the transverse kink motions of magnetic structures in the chromosphere and corona has received a great deal of attention, mainly due to their suitability for energy transfer and also their ability for providing diagnostics of the local plasma environment via coronal or solar magnetoseismology.
Current instrumentation has demonstrated the capabilities required to accurately measure the motions of the finescale magnetic structure, with high spatial and temporal resolution and high signaltonoise levels. Further, the almost continuous coverage of Solar Dynamic Observatory (SDO), supplemented with numerous groundbased data sets and Hinode data, has led to the existence of a large catalogue of wave events. The result of these fortuitous circumstances means that there is plenty of high quality data on kink waves, which can be used for probing the physics of the waves and the local plasma. However, these resources have yet to be exploited effectively, with only simple modelling of individual observed wave events. By modelling, we refer to using the observational data to test theoretical ideas of wave behaviour by fitting an expected model.
In a recent study by Pascoe et al. (2016), the authors attempt to exploit the catalogue of events to do just this. Observations of standing kink waves are utilised to test different models of the damping profile of the oscillatory motion. The authors attempt to determine whether the observed oscillatory motion can be best described by an exponential or Gaussian damping profile, after recent analytical studies of the resonant absorption of kink modes suggested that a Gaussian profile can be more appropriate (Hood et al. 2013; Pascoe et al. 2013). However, the method used by the authors to obtain the timeseries data for model fitting and the technique used for the model comparison are unfortunately inadequate. The technique used by the authors ignores known problems with, and associated uncertainties of, the χ^{2} statistic. As such, there is the potential to end up with erroneous values for model parameters, as well as underestimating the associated uncertainties. More importantly, this methodology is not suitable for distinguishing between the two nonlinear models.
The study of Pascoe et al. (2016) is not the first to fit nonexponential profiles to damped kink oscillations. De Moortel et al. (2002) measured the decay in amplitude of a kink oscillation from wavelet analysis of a displacement timeseries. A leastsquares fit of a function of the form exp(−ϵt^{n}) is performed for three cases, finding values of n = { 1.79,2.83,0.42 }. While the fit visually appear to describe the data well, the uncertainties are not displayed. The amplitude spectrum obtained from Fourier techniques is not a consistent statistical estimator of the true amplitude spectrum of the oscillatory signal, largely because of the noise present in the signal, but also the choice of windowing functions and mother wavelets can cause problems owing to for example, spectral leakage. A discussion of the significance of Fourierbased methods takes us away from the central theme of this manuscript, so we do not expand on this any further.
Verwichte et al. (2004) also investigate whether the damping profile of kink waves observed in a coronal loop arcade is exponential. A trend is subtracted from a displacement timeseries before a weighted nonlinear leastsquares fit of damped sinusoids is performed. The damping term has the same form as De Moortel et al. (2002), but the values of n were chosen, namely n = 1,2,3. The authors state they cannot distinguish between the different damping profiles from the fits, although it is not evident how this comparison was made.
In the following, we outline a methodology for comparing models to observations, providing statistical estimators that can be used to assess the suitability of different models. The pitfalls of the methodology used in previous studies that perform leastsquares fitting are highlighted. We demonstrate that, for a particular example of an oscillating coronal loop, the evidence suggests that the Gaussian damping is the better profile to explain the observed attenuation. Further, the model comparison also reveals that the oscillatory signal contains signatures for the dynamic evolution of the background plasma, signified by timedependent behaviour of the wave period.
Fig. 1 Transverse displacement of a quiescent chromospheric fibril. The timeseries of 141 data points represent the location of the centre of the fibrils’ crosssectional flux profile, with the uncertainties in the location given by the error bars (left panel). The overplotted red dashed line shows the best fit of the model. The distribution of the χ^{2} statistic for 500 realisations of the random noise in the timeseries (right panel). 

Open with DEXTER 
2. Observations and measurement
Here, two different data sets are used for the purpose of demonstrating the techniques for modelling. The first is a Hα data set from the 5 May 2013 taken at the Swedish Solar Telescope (SST) with the CRisp Imaging SpectroPolarimeter (CRISP −Scharmer et al. 2008). Details of the data reduction technique are given in Kuridze et al. (2015). The cadence of the observations is 1.34 s and spatial sampling is 0.̋0592 pixel^{1}.
The second data set is 171 Å data from SDO Atmospheric Imaging Assembly (AIA) recorded on 7 January 2013. The loop analysed is Event 43 Loop 4 using the notation from the catalogue of standing kink oscillations identified in Goddard et al. (2016). The timedistance diagram is created from a location that is high along the loop leg, relatively close to the loop apex. The data was prepared using the standard AIA solarsoft packages. The cadence of the observations is 12 s and spatial sampling is 0.̋6 pixel^{1}.
The technique used to measure the waves from these data sets involves the fitting of a Gaussian function to the crosssectional flux profile of the structures in each timeframe, weighting the fits with uncertainties in intensity. The fitting of a Gaussian is an established and common technique for estimating the location of the central position of an oscillating structure (e.g., Aschwanden et al. 2002; Verwichte et al. 2004). Full details of the technique used here can be found in Morton (2014) and the code is freely available (Morton et al. 2016). The estimates for the uncertainties are obtained from: (i) the equations for noise in AIA data (e.g., Yuan & Nakariakov 2012); (ii) image processing techniques to isolate and model the noise in the SST data (Mooroogen et al. 2016).
For both data sets, an estimate of residual motion between frames is required. To assess this, suitable regions of the full field of view (FOV) are identified and crosscorrelation analysis is performed giving displacement vectors. The displacement vectors contain signatures of the longterm, largescale movement of the chromospheric/coronal features due to underlying photospheric flows and/or solar rotation. As such, we detrend the displacement vectors with high order polynomials to remove such motions, leaving essentially a timeseries of residual frametoframe motion that is thought to be present due to either atmospheric seeing (SST) or spacecraft jitter (SDO). This is performed for a number of regions in each data set, and the root mean square (RMS) of each residual motion timeseries is calculated. The RMS values of motion in the x and y directions are summed in quadrature and averaged between the different regions to provide a single value for each data set, which we take to represent the residual motion. This value is added in quadrature to the uncertainties of all Gaussian centroid positions measured from the crosssectional flux profiles.
As we highlight in the next section, this technique probably under estimates the residual motion in the groundbased data. This is likely due to the processing of the data via the MOMFBD and alignment techniques (e.g., van Noort et al. 2005) that ‘destretches and warps’ relatively small subsections of the data to remove the influence of atmospheric seeing conditions, which can vary between different parts of an image. As such, regions of worse seeing may result in increased difficulties in aligning data between frames and lead to larger uncertainties on the relative positions of magnetic features between frames.
3. Model comparison
Here, we utilise leastsquares minimisation to fit a model to the observed data. The technique involves comparing some model f with parameters a to N data values, y(t_{i}), with errors σ_{i} measured at t_{i} by essentially minimising the χ^{2}, given by (1)A standard approach to model comparison is to find a single value of the reduced χ^{2}, i.e., , where d.o.f. is the number of degrees of freedom, and compare the values of for the different models. However, as discussed in Andrae et al. (2010), the use and definition of is subject to various problems if the model is nonlinear and, as such, is not a suitable quantity for model comparison. Further, it is clear from Eq. (1) that the value of the χ^{2} is subject to uncertainty if the data contains random noise, i.e., random noise can lead to variations in measured values of y(t_{i}), hence, the value of χ^{2} is also dependent upon the data noise.
As a demonstration of this, we take an example of a kink wave timeseries measured in the chromosphere. In Fig. 1, the displayed data is the measured displacement of the central axis of a fibril in the quiescent chromosphere. This example is chosen because the wave motion shows a relatively simple displacement profile that is very near sinusoidal and we do not have to concern ourselves (too much) about extra physics. That data is fit with a simple model of the form, (2)where a_{n} are the parameters to be determined. Fitting the data using the Levenberg algorithm (mpfit −Markwardt 2009), with the determined uncertainties, provides a χ^{2} value of 1910.
To demonstrate the variability of the χ^{2}, we use resampling to create new versions of the timeseries. Each data point, y(t_{i}), has an associated uncertainty (σ_{i}), implying the measured value of y(t_{i}) is the “actual” value plus some deviation due to the random noise, with the value of noise coming from a Gaussian distribution centred on the “actual” value with variance . For each data point y(t_{i}), we generate new values by randomly sampling the noise from a Gaussian distribution with the data points given uncertainty. This process is performed to create five hundred representations of the timeseries. The model is fit to each representation and the χ^{2} is calculated. This process mimics the influence of the random error associated with each data point and Fig. 1 shows the resulting distribution for the χ^{2}, which has a mean and standard deviation of 1988 ± 90. Hence, it is clear that quoting a single value for χ^{2} is not adequate for determining the certainty associated with the fitted model.
As mentioned earlier, from examination of the variability of the timeseries and the given uncertainties, the estimated values for some of the σ_{i} do not represent the observed variance between timeframes particularly well (left panel Fig. 1). There are a number of large excursions, lying at least 5σ from the fitted model. This is likely due to an underestimation of the residual motion of the data between frames in this region. As such, we are likely overestimating the magnitude of χ^{2} associated with this data. It is worth noting that underestimation of σ_{i} leads to a wider distribution of the χ^{2}, as well as overestimating its magnitude.
There are a number of well established statistical techniques to overcome the unsuitability of the reduced χ^{2} and enabling model comparison to take place. A brief discussion of such techniques is given in Andrae et al. (2010) and we do not detail the various techniques here. By way of example, in the following we perform model comparison for a standing coronal kink oscillation, whose measured timeseries shows clear indications that detailed physics is required in the model in order to describe the observation.
3.1. Example: exponential vs. Gaussian damping
The first step in model comparison is choosing an appropriate model for the timeseries that is based upon physical considerations. In analysing the loop oscillations, Pascoe et al. (2016) first detrends the oscillations using a spline fit to the maxima and minima of the oscillation. Detrending signals is a popular technique used in timeseries analysis (e.g., Verwichte et al. 2004, 2010; van Doorsselaere et al. 2007) but there are a number of problems that arise with detrending a series before model fitting (see also, e.g., Vaughan 2005; and Gruber et al. 2011, for similar examples with significance testing of Fourier components). First, it removes physics from the data. While the consequences of detrending using low order polynomials may be small on the values of fitted parameters, complicated spline fits of the data are likely to skew results and influence model comparison, since the data can be manipulated in any fashion. On viewing the subtracted “trends” in Pascoe et al. (2016), it appears physical effects are being removed, with the trendlines showing oscillatory behaviour. Second, the fitted trend influences the maximum likelihood values for the model parameters obtained, hence, the residuals and the value of χ^{2}. Since the goal is to derive accurate estimates for the model parameters and their uncertainties, the removal of the trend is counterproductive and leads to an underestimation of the uncertainty with fitted parameters. As such, the trend should not be removed before model fitting, but should form part of the model fitting and the comparison process.
Fig. 2 Model fitting a kink oscillation of a coronal loop. The top panel shows the timedistance diagram for the coronal loop oscillation, where the time shown is in seconds is from 06:29:59 UT. The bottom panel shows the timeseries of 97 data points representing the location of the centre of the coronal loop’s crosssectional flux profile used for the model fitting (solid black curve) with the error bars showing the standard deviation of uncertainties. The time shown is in seconds is from 06:41:23 UT and the corresponding section of the timedistance diagram is highlighted by the dashed white lines in the top panel. The red dashed curve shows the best fit quartic Gaussian model and the green dashdot curve shows the timedependant model, f_{3}. 

Open with DEXTER 
Before performing the model comparison, we select a portion of the timeseries to fit (Fig. 2). We exclude the beginning of the timeseries from 06:29:59 UT to 06:41:23 UT when the wave is excited. If we did not, the model would have to be complicated enough to include the physics related to the excitation of the wave. Since the modelling of the excitation process is not the focus here, the beginning of the timeseries is excluded. We also exclude the last section of the timeseries from 07:01:11 UT. From here, the displacement amplitude of the wave becomes comparable to the size of the uncertainties, hence this section of the timeseries can add no new information to the model fitting.
We begin by choosing a loworder polynomial as the base for the model. Coronal loops are likely subject to long timescale motions from the shuffling and movement of its footpoints by the horizontal flows in the photosphere. An exact physical model for the trend is unknown (e.g., Aschwanden et al. 2002). It is likely that such motions will be nonlinear, hence we opt for a maximum of a fourthorder polynomial to describe any longterm evolution of the loop. In order to demonstrate the influence of the “trend” on the fitted parameters, we also fit a cubic model.
Further, the data shows evidence for damped sinusoidal displacement of the coronal loop, hence a damped sinusoid is required in the model. The damping of kink waves in coronal loops is a well studied phenomenon and resonant absorption is the most promising candidate to explain the observed damping (e.g., Ruderman & Roberts 2002; Goossens et al. 2002; Hood et al. 2013; Pascoe et al. 2013). Finally, we would like to assess two models that are associated with resonant absorption, one with a Gaussian damping profile and another with an exponential profile, hence we multiply the sinusoid by one of these damping terms. The two sets of models for comparison are then: where f_{1} is the exponential profile and f_{2} the Gaussian profile. Here, for the exponential profile a_{3} = 1 /τ_{e} where τ_{e} is the damping time, and for the Gaussian profile where τ_{g} is the damping time.
Fig. 3 Distribution of the χ^{2} statistic for the quartic and cubic models. The solid lines show the distribution for the exponential damping, while the dashed line shows the distribution for the Gaussian damping. 

Open with DEXTER 
Fig. 4 Distributions of fitted parameters from the quartic Gaussian damping (solid black), cubic Gaussian damping (dashed green) and quartic Gaussian damping with timedependent period (dashdot red). 

Open with DEXTER 
These models are fit to the data separately. As demonstrated in the previous section, the χ^{2} is subject to uncertainties. To assess this, we use resampling of the uncertainties to obtain distributions for the χ^{2} values (Fig. 3) and the fitted parameters of the model (some of which are shown in Fig. 4). The mean values of χ^{2} and their standard deviation are given in Table 1. The χ^{2} values appear to suggest that the Gaussian damping profile performs better then exponential damping profile independent of the trend line used, minimising the sum of squared residuals. It is worth highlighting here that if we use a single of χ^{2} to judge the quality of fit, in a small number of cases it is found that the exponential model out performs the Gaussian, i.e., has a smaller value of χ^{2}. If one of these realisations of the random error just happened to be the measured values, then relying on a single value of the χ^{2} would lead us to wrongly conclude the exponential is a better model.
However, it is not clear that either model is a particularly good approximation to the real data. If the model contains all the physics encapsulated in the data, then we should expect that the normalised residuals should be random and distributed normally, i.e., . In Fig. 5 we show the average normalised residuals between the resampled timeseries and models. The residuals are quite clearly not randomly distributed, with apparent periodic structure present. A nonparametric statistical test called the runs test for randomness (e.g., Barlow 1989; Wall & Jenkins 2003) can be implemented for the residuals. The test begins by assigning a 0 or 1 to the residuals if it lies below or above the mean value respectively. Should the residuals be random and normally distributed, the average number of runs of 0 s and 1 s and the variance for the expected number of runs can be estimated. Comparing this to the observed number of runs forms the hypothesis test. Perhaps predictably, there is evidence at the 95% level to suggest that the residuals are not randomly distributed.
In order to quantify the quality of the fitted models, we perform a residuals test that compares the distribution of the normalised residuals to a normal distribution. The model that is favoured by the data is the one whose normalised residuals best match the normal distribution. A powerful test to perform this comparison is the nonparametric KolmogorovSmirnoff (KS) one sample test (e.g., Babu & Feigelson 1996; Wall & Jenkins 2003), which compares the theoretical cumulative distribution function (CDF) to the sample empirical (E)CDF, examples of which are given in Fig. 6. The test uses the null hypothesis that the sample distribution is from the theoretical CDF and calculates the maximum deviation between the sample and theoretical CDF as the test statistic. The value of the maximum deviation can be used to calculate a probability, which essentially provides an estimate of the probability required for the true model to produce a residual equal to the maximum deviation. Hence, small pvalues imply that there is evidence that the null hypothesis is false, i.e., the residuals are not from a normal distribution. In Fig. 7, we display the pvalues from the KS test. As with the χ^{2}, the pvalue is also subject to uncertainty, so we calculate the value for each of the 500 resamples. It is clear from the values in Fig. 7 that non of the models used appear to accurately describes the physics in the data, although the Gaussian damping profile is favoured more strongly over the exponential damping profile. The mean value of the distributions is given in Table 1 and agrees with the impression from the plotted distributions.
Parameter estimates and statistics from the model comparison.
Focusing now on the model parameters obtained from the fitting, we demonstrate the importance of including the trend line in the fitting of the model. It is evident from Fig. 3 that the χ^{2} value depends on the trend, hence, any further tests that may be carried out for goodness of fit using the residuals will also be dependent upon the trend used. Further, it is found that whether fitting with a cubic or quartic spline influences the values of the model parameters. As an example, it is found that the initial amplitude of the wave is greater using the quartic trend rather the cubic (see also Table 1). We can further quantify any differences between the two distributions of amplitudes from the two fitted trends by testing whether they are from the same parent distribution. This involves using the nonparametric twosample KS test to compare the two ECDFs, whose null hypothesis is that the two distributions are from the same parent population. The test results provide evidence at the 95% level to suggest that the two distributions are from different parent distributions.
Fig. 5 Averaged normalised residuals for the different models with the fourthorder polynomial trend. The black solid line shows the exponential damping profile, the blue dashed line is the Gaussian damping profile and the red dashdot is the timedependent model. The normalised error bars are also displayed for every third data point. 

Open with DEXTER 
As an additional test, we performed the leastsquares analysis on trend subtracted data. First, a quartic trend is fitted by weighted leastsquares to the timeseries and the found trend is removed. In line with previous studies that remove the trend first, the uncertainties associated with the fitted trend line are not added to the original errors. Performing a leastsquares fit of a damped sinusoid to the trend subtracted series, weighted with the original errors, gives similar values for parameters but underestimates their uncertainties by 20%.
Hence, it is clear that not taking into account the trend will bias the results of the model fitting and also any further tests of the goodness of fit of the models.
Fig. 6 Examples of the KS one sample test. The ECDF of the residuals for the cubic Gaussian damping (red), timedependent quartic Gaussian damping (blue) and normal CDF (dashed). The solid vertical lines show the location of the maximum distance between the ECDF and CDF. 

Open with DEXTER 
3.2. Example: timedependence
As indicated by examination of the residuals, it is suggested that the physics encapsulated in the fitted model did not accurately describe the observed oscillation. There is potentially much physics we could incorporate into the model, but perhaps the most obvious exclusion is the dynamical behaviour of the background plasma. It is now well documented that a significant proportion of observed coronal loops undergo evolution, continually out of hydrostatic equilibrium, exhibiting cycles of heating and cooling where mass is frequently exchanged with the lower solar atmosphere through evaporation and condensation (e.g., Aschwanden et al. 2000; Winebarger et al. 2003; Aschwanden & Terradas 2008; UgarteUrra et al. 2009; Viall & Klimchuk 2012; Tripathi et al. 2012; McIntosh et al. 2012). Oscillations and waves can exist on this dynamically evolving background plasma, and the variation in plasma properties has been demonstrated to alter both the period and amplitude of waves (e.g., Morton & Erdélyi 2009; Morton et al. 2010; Ruderman 2010; Terradas et al. 2011; Chen et al. 2014).
It is therefore natural to extend our model to include dynamic effects. The exact form of the timedependence for a particular physical mechanism is likely to be cumbersome and we do not know which dynamical effects are occurring. Hence, a simple approximation for the dynamics is employed, assuming that the period is timedependent with the form P_{0}(1 + t/τ_{d}), where P_{0} is the initial period and τ_{d} represents the dynamic timescale of change. So, we now fit a model of the form, (5)The mean values of the parameter distributions from the fitting of the resampled data are found and used to define the best fit model, which is overplotted in Fig. 2. The value obtained for a_{2} = 1.36 ± 0.14 × 10^{4} s^{1}, meaning the period increases from 247 s to 288 s over the course of the timeseries. Close inspection reveals that this model fits the data better than any of the previous models. This is supported by the distribution of the χ^{2} statistic, whose mean occurs at significantly smaller values than the previous models (Table 1). Further, the residual test results have dramatically improved, with the mean value of KS test statistic increasing by a factor of 10 compared to the timeindependent quartic Gaussian profile (Fig. 8). It worth noting that ~58% of the resampled series have KS test statistic values that fall below the critical value at the 95% level, i.e., in these cases there is no evidence to suggest that residuals are not normally distributed.
A better way to visualise the comparison of residuals is to show the distribution of normalised residuals in relation to normal distribution. In Fig. 9 the results from all the models used are shown. It is clear for the timeindependent models that the distribution of residuals is much wider than the normal distribution, with significantly more values greater than 3 standard deviations from the mean than expected. The timedependent model however, has much less extreme values of the residuals and can be seen to be approaching a normal distribution.
Fig. 7 Histogram of pvalues from KS test. The solid and dashed lines show the values for the cubic and quartic polynomials, respectively. 

Open with DEXTER 
4. Discussion
The observation of wave phenomenon throughout the solar atmosphere now occurs with regularity, and techniques have been developed to provide high quality measurements of the observed waves properties. The quality of these measurements can be utilised to probe the local plasma of the wave guide through seismology (e.g., Morton 2014; Wang et al. 2015). Significantly, this enables models of oscillatory phenomenon to be tested against observed cases, with the quality of available data enabling us move beyond very basic models, e.g., a simple sinusoid, and probe increasingly complex physics. However, in order to make decisions about whether the chosen models describe the observed behaviour well, robust techniques for model comparison are required. These techniques give an indication of the probability that the observed behaviour will arise, assuming the chosen model is correct. A range of such techniques have already been developed and are used widely in astrophysics, and many other areas.
Fig. 8 Histogram of pvalues from KS test for the timedependent model. 

Open with DEXTER 
Here, we demonstrate the benefit of using these techniques for the modelling of oscillations of finescale magnetic structure in the corona. The unsuitability of single value χ^{2} and also statistics to compare models is highlighted, as discussed in Andrae et al. (2010), and we demonstrate how greater confidence in the comparison of competing models can be achieved by using robust statistical techniques. In particular, the techniques were applied to the kink oscillations of a coronal loop. It was found that the amplitude profile of the kink wave was best described by a Gaussian damping profile rather than an exponential damping profile. This result had already been suggested in Pascoe et al. (2016), although the methodology used in that study to reach this conclusion is based upon comparison of single values of , and as such, no confidence intervals where associated with the conclusions. Moreover, in a number of cases, Pascoe et al. (2016) suggest the exponential profile is a better fit to the data than a Gaussian profile. These cases could potentially be false results, due to either: (i) the subtraction of a spline profile that distorts the true amplitude envelope; (ii) reliance upon a single value of χ^{2} (or more precisely ); (iii) not including the trend line in the assessment of uncertainties.
Fig. 9 Distribution of normalised residuals for all the fitted models. The black line shows a normal distribution. The legend for the different distributions is as follows: quartic Gaussian (purple); quartic exponential (yellow); cubic Gaussian (green); cubic exponential (blue); timedependent (red). 

Open with DEXTER 
Similar comments can be made regarding previous results that use trend subtraction and then try to estimate other parameters. For example, van Doorsselaere et al. (2007) aim to measure multiple harmonics of a kink oscillation associated with a coronal loop. The displacement timeseries is trend subtracted and fit with a single frequency sinusoid. The residuals between the trend subtracted data and single frequency model are then analysed by an additional leastsquares fit of a secondary sinusoidal component. We suggest that the uncertainties associated with the parameters from the fitting of these residuals would be significantly greater than given by the authors. The subtraction of the trend alone is seen to underestimate the uncertainties (Sect. 3.1), but an additional fit to the residuals will only exacerbate this effect on the parameter estimates of the secondary sinusoidal model. To ensure statistical significance, we suggest a single model incorporating all the physics should be fit to the original timeseries.
It is then unclear whether the secondary harmonic found in van Doorsselaere et al. (2007) is statistically significant. The lack of error bars on the data and residual plots also do not enable a visual assessment of the situation. Applying both the runs test and the KS test would help elucidate whether any additional structure was contained within the residuals above the uncertainties.
Finally, we extended the Gaussian damping model to include the effects of timedependence, namely through a simple modification of the period to permit it to vary as a function of time. Such a model naturally arises when the background plasma that supports the oscillation is subject to thermodynamic evolution, e.g., heating/cooling (cf. Morton & Erdélyi 2009; Ruderman 2010, 2011). We note that the amplitude of the oscillation is also influenced by the timedependent background plasma, indicating that the amplitude envelope described by the Gaussian damping contains the influence of both dynamics and attenuation by e.g., resonant absorption. The analysis performed here suggests that the timedependent model describes the data better than the static model. We believe this is evidence of dynamic coronal loop plasma influencing the properties of the oscillation. A change in periodicity over the observed displacement timeseries has previously been reported by De Moortel et al. (2002) and White et al. (2013).
Further work is need to untangle the information that the fitted model provides about the plasma evolution and the attenuation due to resonant absorption, and should be complemented with analysis of the plasma and magnetic field, e.g., differential emission measure. Moreover, while the timedependent model performs much better than the timeindependent model, it is still clear that something is missing from the current analysis (e.g., as demonstrated by the nonrandom residuals in Fig. 5 and the KS statistic in Fig. 8). This could potentially be due to underestimates of uncertainty, which would imply are model describes all the physics occurring during this event. However, we believe it is more likely that some physics in the data is not captured by the model and this is supported by the lack of randomness in the residuals (Fig. 5). It is unknown what this may be at present. One possibility is the data contains the signature of higher harmonics of the kink wave, which could also have been excited along with the fundamental mode (e.g., Verwichte et al. 2004; van Doorsselaere et al. 2007; De Moortel & Brady 2007; O’Shea et al. 2007; Verth et al. 2008). In order to find answers to some of these questions, an extended study will be required.
Acknowledgments
R.M. is grateful to the Leverhulme Trust for the award of an Early Career Fellowship. K.M. thanks Northumbria University for the award of a Ph.D. studentship. The authors acknowledge IDL support provided by STFC. The authors also thank V. Henriquez for allowing access to the SST data.
References
 Andrae, R., SchulzeHartung, T., & Melchior, P. 2010, ArXiv eprints [arXiv:1012.3754], unpublished [Google Scholar]
 Aschwanden, M. J., & Terradas, J. 2008, ApJ, 686, L127 [NASA ADS] [CrossRef] [Google Scholar]
 Aschwanden, M. J., Nightingale, R. W., & Alexander, D. 2000, ApJ, 541, 1059 [NASA ADS] [CrossRef] [Google Scholar]
 Aschwanden, M. J., de Pontieu, B., Schrijver, C. J., & Title, A. M. 2002, Sol. Phys., 206, 99 [NASA ADS] [CrossRef] [Google Scholar]
 Babu, G. J., & Feigelson, E. D. 1996, Astrostatistics (Chapman and Hall) [Google Scholar]
 Barlow, R. 1989, Statistics, A guide to the use of statistical methods in the physical sciences (J. Wiley & Sons, Inc.) [Google Scholar]
 Chen, S.X., Li, B., Xia, L.D., Chen, Y.J., & Yu, H. 2014, Sol. Phys., 289, 1663 [NASA ADS] [CrossRef] [Google Scholar]
 De Moortel, I., & Brady, C. S. 2007, ApJ, 664, 1210 [NASA ADS] [CrossRef] [Google Scholar]
 De Moortel, I., Hood, A. W., & Ireland, J. 2002, A&A, 381, 311 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Goddard, C. R., Nisticò, G., Nakariakov, V. M., & Zimovets, I. V. 2016, A&A, 585, A137 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Goossens, M., Andries, J., & Aschwanden, M. J. 2002, A&A, 394, L39 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gruber, D., Lachowicz, P., Bissaldi, E., et al. 2011, A&A, 533, A61 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hood, A. W., Ruderman, M., Pascoe, D. J., et al. 2013, A&A, 551, A39 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kuridze, D., Henriques, V., Mathioudakis, M., et al. 2015, ApJ, 802, 26 [NASA ADS] [CrossRef] [Google Scholar]
 Markwardt, C. B. 2009, in Astronomical Data Analysis Software and Systems XVIII, eds. D. A. Bohlender, D. Durand, & P. Dowler, ASP Conf. Ser., 411, 251 [Google Scholar]
 McIntosh, S. W., Tian, H., Sechler, M., & De Pontieu, B. 2012, ApJ, 749, 60 [NASA ADS] [CrossRef] [Google Scholar]
 Morton, R. J. 2014, A&A, 566, A90 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Morton, R. J., & Erdélyi, R. 2009, ApJ, 707, 750 [NASA ADS] [CrossRef] [Google Scholar]
 Morton, R. J., Hood, A. W., & Erdélyi, R. 2010, A&A, 512, A23 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Morton, R. J., Mooroogen, K., & McLaughlin, J. A. 2016, NUWT: Northumbria University Wave Tracking (NUWT) code, http://dx.doi.org/10.5281/zenodo.49563 [Google Scholar]
 O’Shea, E., Srivastava, A. K., Doyle, J. G., & Banerjee, D. 2007, ApJ, 473, L13 [NASA ADS] [CrossRef] [Google Scholar]
 Pascoe, D. J., Hood, A. W., De Moortel, I., & Wright, A. N. 2013, A&A, 551, A40 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Pascoe, D. J., Goddard, C. R., Nisticò, G., Anfinogentov, S., & Nakariakov, V. M. 2016, A&A, 585, L6 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ruderman, M. S. 2010, Sol. Phys., 267, 377 [NASA ADS] [CrossRef] [Google Scholar]
 Ruderman, M. S. 2011, Sol. Phys., 271, 41 [NASA ADS] [CrossRef] [Google Scholar]
 Ruderman, M. S., & Roberts, B. 2002, ApJ, 577, 475 [NASA ADS] [CrossRef] [Google Scholar]
 Scharmer, G. B., Narayan, G., Hillberg, T., et al. 2008, ApJ, 689, L69 [NASA ADS] [CrossRef] [Google Scholar]
 Terradas, J., Arregui, I., Verth, G., & Goossens, M. 2011, ApJ, 729, L22 [NASA ADS] [CrossRef] [Google Scholar]
 Tripathi, D., Mason, H. E., Del Zanna, G., & Bradshaw, S. 2012, ApJ, 754, L4 [NASA ADS] [CrossRef] [Google Scholar]
 UgarteUrra, I., Warren, H. P., & Brooks, D. H. 2009, ApJ, 695, 642 [NASA ADS] [CrossRef] [Google Scholar]
 van Doorsselaere, T., Nakariakov, V. M., & Verwichte, E. 2007, A&A, 473, 959 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 van Noort, M., Rouppe van der Voort, L., & Löfdahl, M. G. 2005, Sol. Phys., 228, 191 [NASA ADS] [CrossRef] [Google Scholar]
 Vaughan, S. 2005, A&A, 431, 391 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Verth, G., Erdélyi, R., & Jess, D. B. 2008, ApJ, 687, L45 [NASA ADS] [CrossRef] [Google Scholar]
 Verwichte, E., Nakariakov, V. M., Ofman, L., & Deluca, E. E. 2004, Sol. Phys., 223, 77 [NASA ADS] [CrossRef] [Google Scholar]
 Verwichte, E., Foullon, C., & Van Doorsselaere, T. 2010, ApJ, 717, 458 [NASA ADS] [CrossRef] [Google Scholar]
 Viall, N. M., & Klimchuk, J. A. 2012, ApJ, 753, 35 [NASA ADS] [CrossRef] [Google Scholar]
 Wall, J. V., & Jenkins, C. R. 2003, in Practical Statistics for Astronomers, Cambridge observing handbooks for research astronomers (Cambridge, UK: Cambridge University Press), 3 [Google Scholar]
 Wang, T., Ofman, L., Sun, X., Provornikova, E., & Davila, J. M. 2015, ApJ, 811, L13 [NASA ADS] [CrossRef] [Google Scholar]
 White, R. S., Verwichte, E., & Foullon, C. 2013, ApJ, 774, 104 [NASA ADS] [CrossRef] [Google Scholar]
 Winebarger, A. R., Warren, H. P., & Seaton, D. B. 2003, ApJ, 593, 1164 [NASA ADS] [CrossRef] [Google Scholar]
 Yuan, D., & Nakariakov, V. M. 2012, A&A, 543, A9 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
All Tables
All Figures
Fig. 1 Transverse displacement of a quiescent chromospheric fibril. The timeseries of 141 data points represent the location of the centre of the fibrils’ crosssectional flux profile, with the uncertainties in the location given by the error bars (left panel). The overplotted red dashed line shows the best fit of the model. The distribution of the χ^{2} statistic for 500 realisations of the random noise in the timeseries (right panel). 

Open with DEXTER  
In the text 
Fig. 2 Model fitting a kink oscillation of a coronal loop. The top panel shows the timedistance diagram for the coronal loop oscillation, where the time shown is in seconds is from 06:29:59 UT. The bottom panel shows the timeseries of 97 data points representing the location of the centre of the coronal loop’s crosssectional flux profile used for the model fitting (solid black curve) with the error bars showing the standard deviation of uncertainties. The time shown is in seconds is from 06:41:23 UT and the corresponding section of the timedistance diagram is highlighted by the dashed white lines in the top panel. The red dashed curve shows the best fit quartic Gaussian model and the green dashdot curve shows the timedependant model, f_{3}. 

Open with DEXTER  
In the text 
Fig. 3 Distribution of the χ^{2} statistic for the quartic and cubic models. The solid lines show the distribution for the exponential damping, while the dashed line shows the distribution for the Gaussian damping. 

Open with DEXTER  
In the text 
Fig. 4 Distributions of fitted parameters from the quartic Gaussian damping (solid black), cubic Gaussian damping (dashed green) and quartic Gaussian damping with timedependent period (dashdot red). 

Open with DEXTER  
In the text 
Fig. 5 Averaged normalised residuals for the different models with the fourthorder polynomial trend. The black solid line shows the exponential damping profile, the blue dashed line is the Gaussian damping profile and the red dashdot is the timedependent model. The normalised error bars are also displayed for every third data point. 

Open with DEXTER  
In the text 
Fig. 6 Examples of the KS one sample test. The ECDF of the residuals for the cubic Gaussian damping (red), timedependent quartic Gaussian damping (blue) and normal CDF (dashed). The solid vertical lines show the location of the maximum distance between the ECDF and CDF. 

Open with DEXTER  
In the text 
Fig. 7 Histogram of pvalues from KS test. The solid and dashed lines show the values for the cubic and quartic polynomials, respectively. 

Open with DEXTER  
In the text 
Fig. 8 Histogram of pvalues from KS test for the timedependent model. 

Open with DEXTER  
In the text 
Fig. 9 Distribution of normalised residuals for all the fitted models. The black line shows a normal distribution. The legend for the different distributions is as follows: quartic Gaussian (purple); quartic exponential (yellow); cubic Gaussian (green); cubic exponential (blue); timedependent (red). 

Open with DEXTER  
In the text 