Issue
A&A
Volume 671, March 2023
Solar Orbiter First Results (Nominal Mission Phase)
Article Number A79
Number of page(s) 11
Section The Sun and the Heliosphere
DOI https://doi.org/10.1051/0004-6361/202245293
Published online 08 March 2023

© The Authors 2023

Licence Creative CommonsOpen Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

This article is published in open access under the Subscribe to Open model. Subscribe to A&A to support open access publication.

1. Introduction

Solar Orbiter is a solar and heliospheric mission led by the European Space Agency (ESA) in partnership with NASA. The scientific payload of Solar Orbiter includes four in situ and six remote sensing instruments, including the Spectrometer/Telescope for Imaging X-rays (STIX). STIX is a hard X-ray imaging spectrometer that detects photons with energies in the range of 4−150 keV, with a 1 keV energy resolution (at 6 keV; Krucker et al. 2020). STIX measures bremsstrahlung emission from solar flares and therefore provides diagnostics on the hottest (≳8 MK) flare plasma (Krucker et al. 2020; Battaglia et al. 2021). This means that STIX is equipped to provide information on the accelerated electrons producing such bremsstrahlung emissions upon Coulomb collisions with ambient ions. STIX has a high time resolution, with the ability to sample down to 0.1 s. It also has a stable background and continuously observes the full solar disc. Thanks to these capabilities, STIX is an excellent instrument for measuring time-varying signatures on short timescales (≳1 s) in the X-ray emission of solar flares.

Such time variations are of great interest because they are related to the fundamental timescales occurring in solar flares, such as energy release and particle acceleration processes, as well as magnetohydrodynamic (MHD) waves and oscillations in or around the flare site. It is imperative to understand the origin and nature of the observed time-varying behaviour in order to achieve a unified solar flare model.

Fast time variations, which are sometimes classified as quasi-periodic pulsations (QPPs), have been observed in the emission from solar flares over the past 50 yr, with some of the earliest studies identifying time-varying signatures in the hard X-ray (HXR) energy range (Parks & Winckler 1969). The particular variations of interest in this work are modulations in the flare intensity-time profiles. The ones classified as QPPs typically have periodicities ranging from a few seconds to several minutes (Zimovets et al. 2021). Sub-second spikes in the HXR flare emission have also been detected (Roberts et al. 1983; Qiu et al. 2012; Knuth & Glesener 2020). In some cases, multiple periods of oscillations have been identified alongside amplitude modulation (Kolotkov et al. 2015; Van Doorsselaere et al. 2016; McLaughlin et al. 2018). Furthermore, Hayes et al. (2016) identified these time-varying signatures in both the impulsive and decay phases, thus challenging our understanding of the flare model. Additionally, studies relating active region and flare properties to QPP periodicities have also been performed (Pugh et al. 2017; Hayes et al. 2020). For example, Hayes et al. (2020) found a positive correlation between QPP period and duration. Despite this, QPP period was found to be independent of flare magnitude. Finally, a significant correlation between QPP period and various ribbon properties was found.

Statistical studies have found a wide range of probabilities that a given > M5 GOES class flare will contain a QPP, with estimates ranging from 30−90% (Simões et al. 2015; Inglis et al. 2017; Dominique et al. 2018; Hayes et al. 2020). It has been made clear that time-varying signatures (whether quasi-periodic or not) are commonly observed phenomena which occur in flare emission across the entire electromagnetic spectrum, from radio waves to γ-rays (Nakariakov & Melnikov 2009; Zimovets et al. 2021). These statistics rely on QPP classifications that search for a statistically significant periodic component in the data and thus exclude all fast-time-varying structures that do not display any significant periodicity. Recent reviews of QPP properties include Zimovets et al. (2021), Kupriyanova et al. (2020), McLaughlin et al. (2018), Van Doorsselaere et al. (2016), Nakariakov & Melnikov (2009).

A number of models have been built to attempt to explain the observed time-varying and oscillatory phenomena in solar flare emission, with at least fifteen different models proposed at present (Zimovets et al. 2021). For a recent review of current models and their observational signatures, we refer to Zimovets et al. (2021) and McLaughlin et al. (2018). According to Kupriyanova et al. (2020), these models can generally be divided into three categories: (1) models in which the observed emission is directly modulated by MHD and electromagnetic (EM) waves; (2) models in which the efficiency of energy release is modulated by MHD waves; and (3) models in which the original energy release process is itself quasi-periodic in nature. Despite the many observations of time-varying behaviour in flares, there are few cases where the underlying mechanism(s) behind the variations are unambiguous. This is often due to observational constraints. In addition to this, many mechanisms have non-unique signatures, making the disambiguation more challenging (McLaughlin et al. 2018; Zimovets et al. 2021).

One of the challenges in the study of time-varying signatures and QPPs in solar flares is the identification of the characteristic timescale or period. Broomhall et al. (2019) analysed the strengths and weaknesses of various state-of-the-art techniques for detecting time-varying signatures classified as QPPs, based on a series of tests. The tests involved a simulated dataset of flare time profiles, some of which include a periodic component. From this analysis, a list of eight recommendations for detection methods was set out. One of these recommendations relates to a common pre-processing step, namely: time-series detrending. Broomhall et al. (2019) showed that detrending can be challenging and if it must be done, it should be performed manually, with a time-dependent smoothing window; otherwise, a bias may be introduced. Furthermore, there is the added complexity of the underlying background trend when analysing coronal time series. Auchère et al. (2016) showed that if the power-law dependency of the Fourier power spectra (also known as red noise) is not considered in detection models, this can lead to false detections of a significant periodic component. Inglis et al. (2015) showed the effect of this power-law component in the Fourier power spectrum and took this red noise component into consideration when detecting QPPs with methods such as AFINO (Inglis et al. 2017). Many time-varying signatures in flares also show non-stationary properties, whereby the amplitude and period change as a function of time (Nakariakov et al. 2019). This is a further challenge to detection methods, for which typical Fourier-based approaches do not work well.

The method developed in this work leverages the recent data from Solar Orbiter’s STIX and takes a new approach to analysing fast-time-varying signatures in solar flare HXR time profiles. Since the modulation depth in the HXR signature is large, we assume that detrending is not a necessary step. Instead, the fast-time-varying behaviour is considered here as a linear combination of individual Gaussian contributions to the total observed signal. To achieve this, first the HXR time series are smoothed using Gaussian Process (GP) Regression. Next, the time profiles are fitted with a linear combination of Gaussians. From the Gaussian decomposition, key characteristics such as the waiting time distribution, time evolution of the peak full width at half maximum (FWHM), and amplitude can be derived. Additionally, this method can detect non-stationary (time-evolving) signatures as the method is agnostic to whether the underlying driver is periodic. The derived timing information can be used to spatially resolve sources and to perform time-dependent spectral analyses of the individual peaks. An important distinction between this method and those mentioned in this section is that the method does not rely on the assumption that the observed oscillatory behaviour is periodic. This method can be applied to HXR time profiles which exhibit fast-time variations that may or may not be identified as quasi-periodic by existing analysis techniques. This aspect broadens the scope of our analysis.

The analysis given in this work is based on four M and X GOES-class flares observed by STIX since September 2021. The opportunity to study time variations in flares has greatly improved thanks to the new observations from STIX on Solar Orbiter, which include hundreds of flares that demonstrate a fast-time variability.

2. Observations

The data set used in this work is from the STIX imaging spectrometer onboard Solar Orbiter. STIX creates images by an indirect Fourier based imaging technique (Krucker et al. 2020). It has a movable attenuator which is inserted during large flares to limit exposure to high count rates from low-energy X-ray photons and to avoid instrument saturation (Krucker et al. 2020). STIX is capable of quantifying the location, spectrum, and energy content of flare-accelerated thermal and non-thermal electrons (Krucker et al. 2020).

STIX has an onboard algorithm which performs the dynamic time binning. This means that the time resolution of data taken by STIX can vary throughout the duration of the flare. The highest possible time cadence is 0.1 s. To date, the highest cadence tested is 0.3 s. It is not possible to take data at 0.1 s for a long time (≳1 h) since this fills up the onboard memory. However, it is possible to run at 0.5 s cadence for as long as desired. As such, the time resolution of the data used for analysis in this work is 0.5 s. The overview plots shown in this work are made using lower latency data, with a dynamically-adjusted temporal resolution, which typically ranges from 2 to 20 s. The energy ranges used vary depending on the individual flare profile and are selected such that they contain mainly non-thermal HXR emission.

STIX has a relatively stable non-solar background and is a full Sun imager. It observes the Sun continuously, unlike its predecessor, the Reuven Ramaty High Energy Solar Spectroscopic Imager (RHESSI), which rotated about its own axis with a period of 4 s (Lin et al. 2002). This rotational period sometimes introduced artificial periodicities in the data, which were difficult to disentangle from true flare fluctuations (Inglis et al. 2011). Due to its orbit, RHESSI had a day-night cycle, which meant that certain sections of a flare were often not observed. These constraints limited RHESSI’s capability for analysing time-varying structures on short timescales. In contrast, STIX is an excellent instrument suited to the task of fast-time variations analysis due to its stable background, high temporal resolution, imaging and spectral capabilities. Furthermore, STIX is on board the Solar Orbiter spacecraft, which is in a unique elliptical orbit about the Sun, observing solar dynamics from a different vantage point than Earth.

3. Methodology

The motivation for this work is to develop a method that systematically identifies and characterises the modulation observed on short timescales in HXR non-thermal emission from solar flares. This characterisation is useful because it enables us to gain information about the timing, shape, and origin of HXR pulsations in a systematic way. The methodology involves two main steps: (1) data pre-processing by means of time series normalisation and signal smoothing and (2) Fitting a linear combination of Gaussians to the pre-processed time series. Each step and its motivation is explained in more detail in the following subsections1.

3.1. Data pre-processing

In this work, a Gaussian process (GP) regression is used to smooth the HXR time profiles. This is a non-parametric, Bayesian approach to regression that uses machine learning techniques to fit the data.

3.1.1. Normalization

To perform GP regression, the first step is normalising the data. To do so, we re-scale the original time profile using scikit-learn’s StandardScalar class (Pedregosa et al. 2011) by subtracting the mean, μ and dividing by the standard deviation, σ, such that the new time profile is characterised by μ′ = 0 and σ′ = 1. Figures 1a,b demonstrates this for the STIX 25−76 keV time series of the SOL2022-05-04 flare.

thumbnail Fig. 1.

Various steps involved in the application of this method for the example M1.2 GOES class flare of SOL2022-05-04. The times shown are in seconds from 2022-05-04 15:16:12 (Earth time).

3.1.2. Gaussian process (GP) regression

Gaussian processes (GPs) are a powerful supervised machine learning tool that provide a means to make predictions about data by incorporating prior knowledge. The most obvious area of application is regression problems, such as time-series forecasting. Overall, GPs are not limited to regression and can also be extended to classification and clustering tasks. In this work, GP regression is used to smooth STIX time profiles. This is necessary to identify the local maxima and minima for fitting. Since GP regression is a non-parametric approach, there is no initial assumption about the functional form of the time profile. In contrast, when performing polynomial regression, we assume that the functional form is a polynomial and then estimates the exact coefficients. Furthermore, GP regression uses Bayesian statistics, meaning that the model begins with an initial guess of the fit and iteratively updates its fit as more data (in this case, the intensity of HXR emission at a given time) is fed to the algorithm. Unlike traditional regression methods, which learn the exact values for each function parameter, GP regression models estimate a probability distribution over all possible fits.

A common issue with traditional methods for smoothing data in search of QPPs is the choice of window size. This choice is typically ambiguous and can lead to erroneous results if a window size is not selected manually for individual flares (Broomhall et al. 2019). Since GP regression is based on Bayesian statistics, the model parameters can be chosen systematically by optimising a loss function, which estimates how well the model fits to the data. Also, GP regression has the added benefit of giving an uncertainty estimate on the fitted smooth curve and, unlike most methods involving smoothing windows, GPs can be applied to unevenly sampled data. This is particularly important for STIX data, as it is binned onboard in a dynamic way to optimise the amount of down-linked data (Krucker et al. 2020). Of course, time series can be interpolated and re-sampled at even cadence, but this introduces additional errors; thus, it is preferable to have a method that does not require such pre-processing.

3.1.3. GP kernel and hyperparameter optimisation

A GP is a collection of random variables such that any finite collection of those random variables has joint Gaussian distributions (Rasmussen 2004). A GP can be entirely described by its mean and covariance function (Rasmussen 2004). This is similar to a Gaussian distribution, which is defined by a mean vector and covariance matrix; however, in a GP the distribution is taken over function space (Rasmussen 2004). In other words, the final best fit obtained using GP regression is given by the maximum likelihood of a multidimensional probability distribution over all possible fits. This distribution is specified by a covariance matrix (also known as a kernel).

In this work, a simple radial basis function (RBF) kernel was chosen because an RBF kernel is a common choice for smooth time series data. No inference is made based on the kernel parameters themselves; the GPs are solely used as a smoothing technique. Hübner et al. (2022) have discussed assigning a physical meaning to kernel parameter choice in the context of QPPs. For that case, it would be important to consider various types of kernels. However, for our purpose, a simple choice of kernel was deemed sufficient.

The kernel chosen for this work has two components: a constant kernel component and a radial basis function kernel (squared exponential kernel) component. It takes the following form:

k ( x i , x j ) = A e d ( x i , x j ) 2 2 l 2 , $$ \begin{aligned} k(x_i,x_j) = Ae^{\frac{-d(x_i,x_j)^2}{2l^2}}, \end{aligned} $$(1)

where A is a constant value, d(xi, xj) is the Euclidean distance between two feature vectors xi, xj, and l is the length scale.

For this kernel choice, there are two hyperparameters to optimise, A and l; here, l roughly describes the length of the (time) scale over which data points are related. For example, for a flare profile with subsecond variations, the length scale parameter is expected to be a lot smaller than for a profile with longer (on the order of seconds) variations. Then, A is the amplitude of the kernel. In addition, there is a parameter α, which is added to the diagonal of the kernel matrix to help with fitting and can be considered as Gaussian measurement noise on the training data.

The choice of kernel and its parameters can be optimised by performing an exhaustive search over a wide range of values and computing the performance of the model using cross validation. This is known as hyperparameter optimisation. In cross-validation, the data is split between training and test data. The training data is used to fit model parameters and the test data is used to assess model performance. In this work the split between training and test data is 80:20. In other words, the model is fit using the training data, and the “goodness of fit” is assessed by comparing the model prediction with the known measured value from the test set. The error on the fit is computed based on a chosen loss function – in this case, this is the mean squared error (MSE). An average MSE value on all test data is used to assess overall model performance. Once the model hyperparameters have been optimised, a smooth time profile is obtained (as shown in Fig. 1c).

3.2. Fitting a linear combination of Gaussians to the HXR profile

In many flare cases, the shape of each HXR pulsation can be aptly fit with a Gaussian function, which gives a symmetric rise and fall. The use of a simple functional form such as a Gaussian is preferable because we can easily derive insightful properties from the mean and standard deviation of each individual peak, including the waiting time between peaks, amplitude, and full width at half maximum (FWHM). To a large degree, the form of a Gaussian accurately models the sudden impulse of non-thermal electrons reaching the chromosphere and interacting to produce non-thermal bremsstrahlung emission. Of course, there are cases in which this form may not be a suitable choice; for instance, when there is clear particle trapping in the coronal loop, leading to asymmetric HXR emission profiles and a longer decay time. However, to a large extent, the choice of this functional form describes the observed HXR emission extremely well. Other reasonable choices would include a triangular pulse or a double exponential (a symmetric exponential rise and decay). In order to account for particle trapping, it could also be useful to consider an asymmetric function, such as one with an exponential rise but with a slower decay rate.

In this step, the smoothed signal is fit by a linear combination of Gaussian functions. The gradient of the time profile is computed and the locations where the gradient is zero are identified as peaks and troughs (see Fig. 1d). Each pair of local maxima and minima are considered as a single Gaussian contribution. From these values, initial parameters for the fitting routine are derived. The time of the peak value is considered to be the mean of the Gaussian and the time between peak and trough is taken as a rough estimate of the FWHM. The height of the curve at the peak time is taken as an initial estimate of amplitude. Scipy’s Curve Fit (Virtanen et al. 2020) routine is used with the derived initial values as input and with upper and lower bounds on the possible parameter space to fit the time profile. The range of reasonable parameter values varies with each flare. A large GOES class flare will have a smaller possible range in peak height for fitting, as the counting statistics are better. In general, the range used is quite large to allow for the curve fitting routine to work effectively, yet it excludes wide Gaussian fits. Without putting a restriction on the range of peak widths, the curve fitting routine may return Gaussian contributions with large FWHM. This is demonstrated in Fig. 10.

4. Results

In Sect. 4.1, we present the detailed results of each step of the method for the SOL2022-05-04 flare, for the purposes of demonstration. The process has been tested on four M- and X-class solar flares observed with STIX since September 2021. In Sect. 4.2, we give an overview of all flares analysed and the results obtained. Finally, Sect. 4.3 presents an analysis of QPP identification and characterisation for each event.

4.1. Method results detailed for the SOL2022-05-04 event

4.1.1. Gaussian process regression

GP regression was applied to the SOL2022-05-04 flare of M1.2 GOES class. A grid search was performed over the hyper-parameter space shown in Table 1. For this flare, the optimal solution was found to be A = 1.8, α = 0.007, and l = 0.08. The smoothed GP prediction for the optimised hyper-parameters is shown in Fig. 1c with a 95% confidence interval.

Table 1.

Model hyperparameter grid search domain for the SOL2022-05-04 flare.

4.1.2. Gaussian fitting

The peaks and troughs of the smoothed SOL2022-05-04 time profile are identified (as shown in Fig. 1d) and used as input into the fitting routine. The fitting routine fits a linear combination of individual Gaussians to the smoothed time profile. Constraints are given on the possible range of Gaussian characteristics, and for the case of SOL2022-05-04, the bounds applied are those shown in Table 2. The resulting fit is shown in Fig. 1e. Overall, a good fit to the data is obtained. However, at t ≈ 60 s, several peaks are not well fit. This is because the smoothing step (see Fig. 1c) has removed some of these peaks from the time series as they are very short-lived. One of the short peaks remains in the smooth profile, although it has a smaller amplitude and, thus, it is also not very well fit in the final Gaussian decomposition. This demonstrates an important drawback of the smoothing step.

Table 2.

Bounds applied to the fitting parameters for SOL2022-05-04 flare.

4.2. Flare overview

4.2.1. SOL2021-09-23 M1.9 GOES class flare

The SOL2021-09-23 flare was an early impulsive M1.9 GOES class flare observed by STIX when Solar Orbiter was at 0.60 AU from the Sun and at an angle of ∼33° to the Sun-Earth line. A case study of this flare by Zoë Stiefel et al. (2023) reveals four non-thermal HXR sources observed by STIX, with two inner sources from the traditional flare loop and outer footpoints which may be associated with the early onset of a filament eruption; however, the evidence on the latter is inconclusive. A Gaussian decomposition for this flare was performed and the resulting fit (shown in Fig. 2) has been assessed by calculating the standard deviation of the normalised residual. The normalised residual is calculated as the difference between the fit and the observed STIX time profile normalised by the error on the measured STIX time profile. The error on the STIX observed time profile is a combination of the compression error and counting statistics error. The fit gives σR = 1.83. The residuals are quite high, particularly at the beginning of the impulsive phase because there are fast-time variation on timescales smoothed out by the Gaussian process regression. However, the fit improves in the later phase.

thumbnail Fig. 2.

Overview of the SOL2021-09-23 flare. The 18−28 keV HXR non-thermal Gaussian decomposition is shown with a standard deviation of the residual σR = 1.83. The first few peaks have substructure on timescales smoothed out in the GP regression output. As a result, variations on these timescales are not well fit.

4.2.2. SOL2021-10-09 M1.6 GOES class flare

The SOL2021-10-09 flare was an M1.6 GOES class flare observed when Solar Orbiter was at a distance of 0.68 AU from the Sun with the spacecraft-Sun-Earth angle of ∼15°. The Gaussian decomposition method gives a fit with standard deviation of the residual σR = 0.62, shown in Fig. 3. In this case, the measured data is noisier than other example flares because the flare is smaller. As a result, it is easier for the fit to appear to perform well on global structures, but the variance among the fits is high.

thumbnail Fig. 3.

Overview of the SOL2021-10-09 flare and the Gaussian decomposition of its 20−25 keV HXR non-thermal time profile. The fit gives a standard deviation of the residual σR = 0.62. This is a strong fit. The noise level on the measured values of the counts is higher for this flare since it is smaller. This makes it easier to fit and thus the residuals are smaller than those of a large flare with higher signal-to-noise ratio.

4.2.3. SOL2022-03-30 X1.4 GOES class flare

The X1.4 GOES class flare observed on 30-03-2022 was the first X class flare to be observed by STIX during perihelion. This observation was taken at a distance of 0.33 AU from the Sun. Furthermore, the angle between Solar Orbiter and the Sun-Earth line at this time was ∼95°. Due to the large flux incident on the instrument during this flare, the aluminium alloy attenuator was automatically inserted when a certain trigger threshold was reached. The periods when the attenuator was inserted are shown in Fig. 4. The attenuator mainly blocks low energy X-ray photons and allows high energy photons to be transmitted. The energy dependent response of the attenuator can be found in Krucker et al. (2020). As the counts incident on the detector increase, the live-time decreases and vice versa when the attenuator is inserted. Therefore, to correctly show the time profile of such a flare, we should carefully account for changes in detector live-time as well as the energy and spectrum dependent attenuator transmission. This correction is non-trivial as it requires detailed knowledge of the intricacies of such an instrument. However, for the background detector (BKG), the pixels with large apertures are turned off when flux is high. This means that the attenuator does not cover this detector when inserted. As such, there is no effect exerted by the attenuator motion on the 4−10 keV background detector time profile shown. The background detector was used for the thermal 4−10 keV profile shown in Fig. 4, as it requires no correction for the attenuator motion.

thumbnail Fig. 4.

Overview plot of the SOL2022-03-30 flare and the Gaussian decomposition of the 32−45 keV HXR profile for three different phases: phases 1, 2, and 3 from left to right. The standard deviation of the residual of each phase is σR, 1 = 1.25, σR, 2 = 1.58 and σR, 3 = 1.11.

This flare is an early impulsive flare which shows time-varying structures. In particular, there are variations on the timescale of several seconds. Later in the flare, during the thermal peak, there are further variations with significantly different fluctuation timescales (∼14 s) and, later on, even longer periods of ∼35 s. The pulsation timescales are changing significantly over the course of the flare and so, it is challenging to smooth the time profile using our simple choice of kernel, since the length-scale parameter has no time dependence. To account for this, the time profile was split into three different sections: the early phase (phase 1), middle phase (phase 2), and late phase (phase 3). Figure 4 shows the Gaussian decomposition for each phase with the standard deviation of the residuals being σR, 1 = 1.25, σR, 2 = 1.58 and σR, 3 = 1.11 for phases 1, 2, and 3, respectively. In particular, phases 1 and 3 appear to be well fitted with a linear combination of Gaussians, but the residuals for phase 2 show sinusoidal variation, indicating the data are not well fit by the model.

4.2.4. SOL2022-05-04 M1.2 GOES class flare

Finally, we present the M1.2 GOES class flare on May 4th 2022. At this time, Solar Orbiter was observing the far side of the Sun from the Earth, namely, the spacecraft-Sun-Earth angle was ∼163° and the distance to the Sun was 0.73 AU. As such, the flare was not observed by Earth based observatories. The GOES class of a flare which was not observed from Earth is estimated using a model that fits STIX 4−10 keV counts to the GOES flux2. Similarly to the SOL2022-03-30 flare, the attenuator motion is shown in Fig. 5 and the background detector was used for the thermal 4−10 keV profile. The standard deviation of the residual σR = 1.39. Overall, the data are well fit with a linear combination of Gaussians. The fastest fluctuations in time are not well fitted (e.g. Fig. 1e at t ≈ 60 s) because the time profile output by the GP regression smooths out some very short variations. Therefore, the Gaussian fitting procedure doesn’t fit these shorter peaks. This is a limitation of the method that is discussed further in Sect. 5.2.

thumbnail Fig. 5.

Overview plot of the SOL2022-05-04 flare and the Gaussian decomposition of the impulsive phase 25−76 keV profile. The standard deviation of the residual for the fit is σR = 1.39.

4.3. QPP identification and characterisation

4.3.1. Line Fitting Method

An estimate of signal (quasi-) periodicity can be derived from individual Gaussian components. Figure 6 shows the time of each Gaussian component against the peak number for the SOL2022-05-04 flare. The slope of the line of best fit gives an estimate for the period, ∼5.55 s. This shows that STIX can detect variations on short timescales, which is something that was not feasible with its predecessor RHESSI. This result is then compared with the AFINO analysis method (Inglis et al. 2017). AFINO fits the Fourier power spectrum of a flare signal with four different models: (1) a power law + constant; (2) power law + constant + Gaussian A broken power law + constant; (3) power law + 2 Gaussians. Should a QPP component be present in the signal, the second or fourth model would be favoured, that is, the Gaussian peak fits the enhanced power due to the QPP that is present. The goodness of fit is estimated from a Bayesian information criterion (BIC), given by:

BIC = 2 ln ( L ) + k ln ( n ) , $$ \begin{aligned} \mathrm{BIC} = -2 \ln (L) + k \ln (n), \end{aligned} $$(2)

thumbnail Fig. 6.

Mean time of each Gaussian component against peak number. The slope of the line fit gives an estimate of periodicity in the signal. A stronger fit with multiple lines would indicate that there is non-stationarity in the signal.

where L is the maximum likelihood, k is the number of free parameters, and n is the number of data points in the Fourier power spectrum. Thus, the BIC score penalises over-fitting. A large negative value for a given model indicates a strong fit to the data. One model is said to be strongly preferred over another if |ΔBIC| > 10. In the case of the SOL2022-05-04 flare, the BIC score for model 1 versus model 2 is |ΔBIC12| = 4.3, whereby model 2 gives a larger negative BIC score. Thus, the AFINO method gives a slight preference for the QPP model with single period, P = 4 . 96 0.54 + 0.66 $ P = 4.96^{+0.66}_{-0.54} $ s, over a simple power law model (see Fig. 7). This is largely consistent with the period obtained from the Gaussian decomposition method, which is slightly higher since smoothing the time profile suppresses the fastest time variations. We notice that when assessed based on the standard AFINO criteria (|ΔBIC| > 10), this flare would not be marked as a QPP detection, although there is clearly enhanced power at f0 = 0.202 ± 0.024 Hz ( P = 4 . 96 0.54 + 0.66 $ P = 4.96^{+0.66}_{-0.54} $ s). From the line fitting method, there is the added benefit of obtaining the exact timing information of each pulsation, which can be used to spatially resolve each peak (see Fig. 6).

thumbnail Fig. 7.

Power spectral density of signal shown in Fig. 1a. The AFINO QPP model best fit is shown. AFINO analysis gives a period of 4 . 96 0.54 + 0.66 $ {\sim}4.96^{+0.66}_{-0.54} $ s, with a moderate preference for the QPP model with |ΔBIC| = 4.3.

By fitting two lines, a test can be performed to check for non-stationarity (changes in quasi-periodicity over time) in the flare signal. The derived slopes are 4.48 s and 6.78 s for the early and late stages, respectively. However, the Pearson correlation co-efficients obtained indicate that a single line fit or single periodicity in the signal gives a better fit to the data. This is, of course, biased by the number of data points in each fit. Further analysis was done where the signal was split into two time ranges, corresponding to the time ranges of the two line fits. AFINO analysis was performed on these two time ranges. In this case, the QPP model was only favoured in the first time range. Therefore, the AFINO method favours a single QPP model over two periods in this particular case. This agrees with the apparent preference for a single periodicity given by the line fit. The line fitting method was also applied to the SOL2021-09-23, SOL2021-10-09 and SOL2022-03-30 flares and the results are shown in Figs. 8 and 9.

thumbnail Fig. 8.

Time of each peak from the Gaussian decomposition against peak number for the SOL2021-09-23 and SOL2021-10-09 flare shown in Figs. 2 and 3, respectively. The fit gives an estimate of periodicity in the SOL2021-09-23 non-thermal HXR flare signal of ∼19 s. For the SOL2021-10-09 flare, a periodicity of ∼128 s is derived.

thumbnail Fig. 9.

Line fits to the peak time over peak number for three phases of the SOL2022-03-30 flare shown in Fig. 4. Since the pulsation frequency is increasing significantly in time, each phase was analysed separately. Phase 1 is very well fitted by a line indicating that the timing of peaks is quasi-periodic, phases 2 and 3 are not as well fit by a line. This agrees with the error analysis presented in Table 3, where phases 2 and 3 have high variance. Between phases 1 and 3, the periodicity has increased by a factor of ∼5. This indicates that the driver of these pulsations could be different in the later stages of the flare, compared to the initial impulsive phase.

It has been demonstrated that this new method is capable of accurately detecting and characterising fast-time-varying structures and QPPs in the non-thermal HXR emission from flares. The results from this method are consistent with Fourier based analysis methods such as AFINO and in addition, allow for the extraction of important information, for instance, the time between peaks and pulsation duration. Further, this method can be used to detect frequency drifts and non-stationarity in time series.

4.3.2. Errors on QPP analysis

The Gaussian decomposition fits are non-unique and multiple solutions can be derived for a given time series. To determine the effect of this, error analysis was performed for determining the QPP period. For this, an array of errors was added to each STIX time series, which is some random integer multiple of the calculated measurement error (that includes the uncertainty introduced due to compression and counting statistics). The random integer is sampled from a normal white noise distribution with zero mean and unit standard deviation. The entire fitting procedure is then performed 50 times, with different initial error arrays added to the time profile. This gives many unique Gaussian decomposition fits for a given flare profile. As previously, the periodicity of the signal is estimated from a single line fit and an average period for each flare is obtained, with a standard deviation over 50 iterations. The mean periods and standard deviation are as shown in Table 3. An estimation of the range of periods that reasonably describe the variations observed is also given.

Table 3.

Mean and range of periods derived for each flare over 50 iterations, where random error has been added to the initial time series.

These results indicate that the period estimation for flares SOL2021-09-23, SOL2022-30-03, and SOL2022-05-04 in Figs. 2, 4P1, and 5, respectively, are strong fits, since the estimated periods are within 1.1σ of the means obtained from the error analysis. However, the SOL2021-10-09 flare has high variance when noise is added to the time profile. This is because the HXR counts are low and as a result the measurement error is relatively high; therefore, when additional error is added, the Gaussian decomposition fits vary significantly. The fit obtained in Fig. 3 is thus unreliable, as the derived period is 2.9σ from the mean period over 50 iterations. We notice that for the SOL2022-03-30 flare, phase 1 gives a strong fit; however, phases 2 and 3 have a wide range of periods – thus, they are not as well constrained.

5. Discussion and conclusions

5.1. First analysis of fast time-varying structures and QPP detection with STIX data

In this work, we present the first analysis of fast-time-varying structures in the non-thermal HXR emission from flares using Solar Orbiter’s STIX instrument. We have developed a new method for identifying and quantifying fast-time-varying structures in the HXR emission from flares. This method decomposes the non-thermal HXR emission from flares in the impulsive phase into a linear combination of Gaussian pulses. For a sample of four M- and X-class flares, the standard deviation of the normalised residual is ≤1.8. This indicates that the model is a good fit to the data for this selection of events.

From these four flares, fast-time variations on timescales ranging from 4−128 s have been characterised. Furthermore, the first detection of solar flare QPPs from STIX observations has been made. It has been shown that QPPs with timescales down to the order of ∼4 s can be detected with STIX, which was not possible with it predecessor RHESSI due to the spacecraft rotational period of ∼4 s. This work demonstrates that STIX is an instrument well suited to the detection of fast-time variations in the HXR emission from solar flares in the timescale range of seconds to minutes, due to its high time resolution and relatively constant non-solar background.

5.2. Drawbacks

In this section, we present and discuss the drawbacks and subtleties associated with the method and how they may impact the derived results:

1. One important drawback of this method is that GP regression smooths time variations on timescales that are shorter than the optimal length scale obtained from hyperparameter optimisation. This means that some peaks that are very short-lasting are smoothed out and, hence, they are not fit by the Gaussian fitting routine. Importantly, the timescales that are suppressed are a function of the sampling rate, since with a higher time cadence, shorter variations will likely be more prominent and vice versa. As such, the drawback of this method is that variations that occur at the limit of the instrument’s sampling cadence are smoothed out. In this case with STIX, variations on short timescales such as second or sub-second are smoothed out. While QPPs can occur on short timescales such as second or sub-second periodicities, particularly in the radio band (Tan & Tan 2012; Nakariakov et al. 2018; Carley et al. 2019), this work focuses on QPPs that have timescales on the order of seconds to minutes, which are very commonly reported (see also McLaughlin et al. 2018; Hayes et al. 2020; Zimovets et al. 2021). We have shown here that the method works well for these types of QPPs.

2. It is well known that in the case of significant particle trapping in a flare loop, the HXR time profile becomes asymmetric. In this case, fitting Gaussian curves to the time profiles of a trapped particle population is no longer physically meaningful. We should consider fitting other functional forms with asymmetric profiles.

3. A subtlety of this method is that although the Gaussian decomposition step is good at characterising non-stationarity in a signal, GP regression with our choice of kernel, is unable to model large frequency drifts in a flare since the length scale is not time-dependent. This gives some limitation to the method when identifying non-stationarity. This effect is particularly pertinent to the case of the SOL2022-03-30 flare, whereby in the early impulsive phase we observe fluctuations on the order of ∼7 s and later during the thermal peak there are fluctuations on a ∼14 s and then a ∼35 s timescale, as shown in Fig. 9. One way to rectify this is to split the time series into sections and perform the smoothing on different time ranges separately, as was done for the SOL2022-03-30 flare. Another possibility is to consider a more complex kernel choice which has some time dependency. This method was also suggested by Hübner et al. (2022).

4. Another important note to consider is that the Gaussian decomposition fit is non-unique and influenced by the bounds applied to the curve fitting routine. Furthermore, the routine fits a fixed number of Gaussians based on the number of local maxima and minima. As such, several closely separated peaks in quick succession may not be fit by this method, if the gradient remains non-zero. Figure 10 demonstrates the effect of applying boundary conditions on the fitting parameters for the SOL2022-05-04 flare.

thumbnail Fig. 10.

Effect of applying boundary conditions on the fit parameters of the Gaussian decomposition for the SOL2022-05-04 flare. Panel (a) shows the resulting fit when the boundary conditions shown in Table 2 are applied. Panel (b) shows the resulting fit when no boundary conditions are applied. Both cases give a good fit to the data, however, the unconstrained fit (Fig. 10b) fits long period Gaussians with small amplitude peaks, for instance, at ∼15:17:08. This fit is interesting from a physics point of view, however in this work we wish to derive timing information regarding the global peak trend – therefore a constrained fit is preferable.

5.3. Applications and future work

The potential applications of this method are wide and varied. Obtaining timing information on individual non–thermal HXR pulses in a systematic way enables a wide-scale, in-depth study of fast-time-varying structures. In particular, the timing information obtained in this method can be used to image individual pulse contributions to the HXR flare profile. With such information, we can begin to understand the spatial structure (location, morphology) of the oscillatory source, as suggested by Zimovets et al. (2021), and make a comparison with those expected for various models. For example, should sausage mode oscillations be responsible for fast-time variations in HXR profiles, we would expect the source morphology to change over time and the location to remain fixed. In contrast to this, if repeated reconnection were to be responsible for such modulations, we would expect the HXR source spatial locations to change in time since the reconnecting field lines must change. The time evolution of such structures can be investigated with STIX’s imaging capabilities. Although it should be noted that imaging fast-time-varying HXR structures can be challenging with an indirect imager such as STIX due to its limited dynamic range; in particular, in cases where there is one source that is much brighter than the others. Furthermore, one can perform time dependent spectral analysis with STIX. This can help to identify and investigate the underlying particle population behind such oscillatory signatures. Additionally, time-dependent imaging spectroscopy can be used to localise particle populations.

This method adds to the set of current QPP detection techniques, with the additional ability of characterising non-QPP, fast-time-varying signatures. The quantities and relationships derived by this method are important as they feed back into modeling efforts. For example, large-scale waiting time distributions could be used to assess whether avalanche models can accurately describe the observed fast-time variation phenomenon.

Future work will focus on using the derived timing information and shape to determine over which interval HXR images should be made in order to understand the source origin. As STIX is an indirect Fourier imager, another interesting avenue of investigation will be to analyse the evolution of Fourier components (visibilities) over time. We will focus on combining the analysis with other datasets at various wavelength of emission. This method will be applied to other HXR datasets such as Fermi Gamma-Ray Space Telescope (Meegan et al. 2009) and RHESSI (Lin et al. 2002), and to future instruments including ASO-S/HXI (Zhang et al. 2019) and Aditya-L1/HELIOS (Seetha & Megala 2017). This method will also be applied to other wavelengths of emission such as pulsations observed in radio by the Expanded Owens Valley Solar Array (EOVSA; Gary et al. 2018), which are often correlated with those seen in HXR emission (Aschwanden et al. 1990).

In conclusion, a new method for identifying and characterising fast-time variations in the non-thermal HXR time profiles of solar flares has been developed, in which the signals are decomposed into individual Gaussian contributions. The fits obtained have a standard deviation of the normalised residual of ≤1.8. The first characterisation and detection of fast-time variations in the HXR profile of solar flares with STIX has been made on timescales between 4 and 128 s.

The opportunity to study time variations in flares has greatly improved with new observations from STIX on Solar Orbiter and its observations of numerous flares that demonstrate fast-time variability over a wide range of timescales.


1

The code used here is publicly available at https://github.com/hannahc243/Gaussian_Decomp

Acknowledgments

Solar Orbiter is a space mission of international collaboration between ESA and NASA, operated by ESA. The STIX instrument is an international collaboration between Switzerland, Poland, France, Czech Republic, Germany, Austria, Ireland, and Italy. H.C., A.F.B. and S.K. are supported by the Swiss National Science Foundation Grant 200021L_189180 for STIX. L.A.H. is supported by an ESA Research Fellowship.

References

  1. Aschwanden, M. J., Benz, A. O., & Kane, S. R. 1990, A&A, 229, 206 [NASA ADS] [Google Scholar]
  2. Auchère, F., Froment, C., Bocchialini, K., Buchlin, E., & Solomon, J. 2016, ApJ, 825, 110 [Google Scholar]
  3. Battaglia, A. F., Saqri, J., Massa, P., et al. 2021, A&A, 656, A4 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  4. Broomhall, A.-M., Davenport, J., Hayes, L., et al. 2019, ApJS, 244, 44 [NASA ADS] [CrossRef] [Google Scholar]
  5. Carley, E. P., Hayes, L. A., Murray, S. A., et al. 2019, Nat. Commun., 10, 2276 [NASA ADS] [CrossRef] [Google Scholar]
  6. Dominique, M., Zhukov, A. N., Dolla, L., Inglis, A., & Lapenta, G. 2018, Sol. Phys., 293, 61 [NASA ADS] [CrossRef] [Google Scholar]
  7. Gary, D. E., Chen, B., Dennis, B. R., et al. 2018, ApJ, 863, 83 [CrossRef] [Google Scholar]
  8. Hayes, L. A., Gallagher, P. T., Dennis, B. R., et al. 2016, ApJ, 827, L30 [NASA ADS] [CrossRef] [Google Scholar]
  9. Hayes, L. A., Inglis, A. R., Christe, S., Dennis, B., & Gallagher, P. T. 2020, ApJ, 895, 50 [Google Scholar]
  10. Hübner, M., Huppenkothen, D., Lasky, P. D., et al. 2022, AJ, 936,17 [CrossRef] [Google Scholar]
  11. Inglis, A. R., Zimovets, I. V., Dennis, B. R., et al. 2011, A&A, 530, A47 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  12. Inglis, A. R., Ireland, J., & Dominique, M. 2015, ApJ, 798, 108 [Google Scholar]
  13. Inglis, A., Ireland, J., Dennis, B. R., Hayes, L. A., & Gallagher, P. T. 2017, AAS/Solar Physics Division Meeting, 48, 400.05 [NASA ADS] [Google Scholar]
  14. Knuth, T., & Glesener, L. 2020, ApJ, 903, 63 [NASA ADS] [CrossRef] [Google Scholar]
  15. Kolotkov, D. Y., Nakariakov, V. M., Kupriyanova, E. G., Ratcliffe, H., & Shibasaki, K. 2015, A&A, 574, A53 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  16. Krucker, S., Hurford, G. J., Grimm, O., et al. 2020, A&A, 642, A15 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  17. Kupriyanova, E., Kolotkov, D., Nakariakov, V., & Kaufman, A. 2020, Solar-Terr. Phys., 6, 3 [NASA ADS] [Google Scholar]
  18. Lin, R. P., Dennis, B. R., Hurford, G. J., et al. 2002, Sol. Phys., 210, 3 [Google Scholar]
  19. McLaughlin, J. A., Nakariakov, V. M., Dominique, M., Jelínek, P., & Takasao, S. 2018, Space Sci. Rev., 214, 45 [Google Scholar]
  20. Meegan, C., Lichti, G., Bhat, P. N., et al. 2009, ApJ, 702, 791 [Google Scholar]
  21. Nakariakov, V., & Melnikov, V. 2009, Space Sci. Rev., 149, 119 [NASA ADS] [CrossRef] [Google Scholar]
  22. Nakariakov, V. M., Anfinogentov, S., Storozhenko, A. A., et al. 2018, ApJ, 859, 154 [NASA ADS] [CrossRef] [Google Scholar]
  23. Nakariakov, V. M., Kolotkov, D. Y., Kupriyanova, E. G., et al. 2019, Plasma Phys. Control. Fusion, 61, 014024 [Google Scholar]
  24. Parks, G. K., & Winckler, J. R. 1969, ApJ, 155, L117 [Google Scholar]
  25. Pedregosa, F., Varoquaux, G., Gramfort, A., et al. 2011, J. Mach. Learn. Res., 12, 2825 [Google Scholar]
  26. Pugh, C. E., Nakariakov, V. M., Broomhall, A. M., Bogomolov, A. V., & Myagkova, I. N. 2017, A&A, 608, A101 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  27. Qiu, J., Cheng, J. X., Hurford, G. J., Xu, Y., & Wang, H. 2012, A&A, 547, A72 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  28. Rasmussen, C. E. 2004, Gaussian Processes in Machine Learning (Berlin, Heidelberg: Springer-Verlag), 63 [Google Scholar]
  29. Roberts, B., Edwin, P. M., & Benz, A. O. 1983, Nature, 305, 688 [NASA ADS] [CrossRef] [Google Scholar]
  30. Seetha, S., & Megala, S. 2017, Curr. Sci., 113, 610 [NASA ADS] [CrossRef] [Google Scholar]
  31. Simões, P. J. A., Hudson, H. S., & Fletcher, L. 2015, Sol. Phys., 290, 3625 [CrossRef] [Google Scholar]
  32. Tan, B., & Tan, C. 2012, ApJ, 749, 28 [NASA ADS] [CrossRef] [Google Scholar]
  33. Van Doorsselaere, T., Kupriyanova, E. G., & Yuan, D. 2016, Sol. Phys., 291, 3143 [Google Scholar]
  34. Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nat. Meth., 17, 261 [Google Scholar]
  35. Zhang, Z., Chen, D.-Y., Wu, J., et al. 2019, RAA, 19, 160 [NASA ADS] [Google Scholar]
  36. Zimovets, I. V., McLaughlin, J. A., Srivastava, A. K., et al. 2021, Space Sci. Rev., 217, 66 [NASA ADS] [CrossRef] [Google Scholar]
  37. Zoë Stiefel, M., Battaglia, A. F., Barczynski, K., et al. 2023, A&A, 670, A89 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

All Tables

Table 1.

Model hyperparameter grid search domain for the SOL2022-05-04 flare.

Table 2.

Bounds applied to the fitting parameters for SOL2022-05-04 flare.

Table 3.

Mean and range of periods derived for each flare over 50 iterations, where random error has been added to the initial time series.

All Figures

thumbnail Fig. 1.

Various steps involved in the application of this method for the example M1.2 GOES class flare of SOL2022-05-04. The times shown are in seconds from 2022-05-04 15:16:12 (Earth time).

In the text
thumbnail Fig. 2.

Overview of the SOL2021-09-23 flare. The 18−28 keV HXR non-thermal Gaussian decomposition is shown with a standard deviation of the residual σR = 1.83. The first few peaks have substructure on timescales smoothed out in the GP regression output. As a result, variations on these timescales are not well fit.

In the text
thumbnail Fig. 3.

Overview of the SOL2021-10-09 flare and the Gaussian decomposition of its 20−25 keV HXR non-thermal time profile. The fit gives a standard deviation of the residual σR = 0.62. This is a strong fit. The noise level on the measured values of the counts is higher for this flare since it is smaller. This makes it easier to fit and thus the residuals are smaller than those of a large flare with higher signal-to-noise ratio.

In the text
thumbnail Fig. 4.

Overview plot of the SOL2022-03-30 flare and the Gaussian decomposition of the 32−45 keV HXR profile for three different phases: phases 1, 2, and 3 from left to right. The standard deviation of the residual of each phase is σR, 1 = 1.25, σR, 2 = 1.58 and σR, 3 = 1.11.

In the text
thumbnail Fig. 5.

Overview plot of the SOL2022-05-04 flare and the Gaussian decomposition of the impulsive phase 25−76 keV profile. The standard deviation of the residual for the fit is σR = 1.39.

In the text
thumbnail Fig. 6.

Mean time of each Gaussian component against peak number. The slope of the line fit gives an estimate of periodicity in the signal. A stronger fit with multiple lines would indicate that there is non-stationarity in the signal.

In the text
thumbnail Fig. 7.

Power spectral density of signal shown in Fig. 1a. The AFINO QPP model best fit is shown. AFINO analysis gives a period of 4 . 96 0.54 + 0.66 $ {\sim}4.96^{+0.66}_{-0.54} $ s, with a moderate preference for the QPP model with |ΔBIC| = 4.3.

In the text
thumbnail Fig. 8.

Time of each peak from the Gaussian decomposition against peak number for the SOL2021-09-23 and SOL2021-10-09 flare shown in Figs. 2 and 3, respectively. The fit gives an estimate of periodicity in the SOL2021-09-23 non-thermal HXR flare signal of ∼19 s. For the SOL2021-10-09 flare, a periodicity of ∼128 s is derived.

In the text
thumbnail Fig. 9.

Line fits to the peak time over peak number for three phases of the SOL2022-03-30 flare shown in Fig. 4. Since the pulsation frequency is increasing significantly in time, each phase was analysed separately. Phase 1 is very well fitted by a line indicating that the timing of peaks is quasi-periodic, phases 2 and 3 are not as well fit by a line. This agrees with the error analysis presented in Table 3, where phases 2 and 3 have high variance. Between phases 1 and 3, the periodicity has increased by a factor of ∼5. This indicates that the driver of these pulsations could be different in the later stages of the flare, compared to the initial impulsive phase.

In the text
thumbnail Fig. 10.

Effect of applying boundary conditions on the fit parameters of the Gaussian decomposition for the SOL2022-05-04 flare. Panel (a) shows the resulting fit when the boundary conditions shown in Table 2 are applied. Panel (b) shows the resulting fit when no boundary conditions are applied. Both cases give a good fit to the data, however, the unconstrained fit (Fig. 10b) fits long period Gaussians with small amplitude peaks, for instance, at ∼15:17:08. This fit is interesting from a physics point of view, however in this work we wish to derive timing information regarding the global peak trend – therefore a constrained fit is preferable.

In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.