The Dependence of Solar Wind Burst Size on Burst Duration and Its Invariance Across Solar Cycles 23 and 24

Time series of solar wind variables are bursty in nature. Bursts, or excursions, in the time series of solar wind parameters are associated with various transient structures in the solar wind plasma and are often the drivers of increased space weather activity in Earth’s magnetosphere. We define bursts by setting a threshold value of the time series and identifying how often, and for how long, it is exceeded. This allows us to study how the statistical distributions and scaling properties of burst parameters vary over solar cycles 23 and 24. We find that the distributions of burst duration and integrated burst size vary over the solar cycle and between the equivalent phases of consecutive cycles. However, there exists a single power law scaling relation between burst size and duration, with a joint area-duration scaling exponent α that is independent of the solar cycle. This provides a solar cycle invariant constraint between possible sizes and durations of solar wind bursts that can occur. Plain Language Summary The solar wind is an outflow of plasma which streams away from the Sun’s atmosphere. Because it carries a magnetic field, when it reaches the Earth, it can interact with the terrestrial magnetic field and transfer energy. This space weather interaction is responsible for the northern and southern lights; however, energy deposited by the solar wind can also have significant negative consequences, such as disrupting radio communications and navigation and in extreme cases inducing large enough currents in power grids to cause blackouts. The solar wind and its magnetic field are bursty, and these bursts can cause space weather. We study how the likelihood of bursts varies over the two most recent 11-year solar activity cycles. The time duration and size of the bursts are found to be related to each other in a manner that does not vary with solar cycle activity, whereas other detailed properties of the bursts do vary. This constrains the size and duration of potentially space weather-causing bursts that can arrive at the Earth.


Introduction
Time series of space plasma variables, such as the interplanetary magnetic field (IMF) strength or the AE geomagnetic index, characteristically take irregularly spaced large excursions above their average value (Consolini et al., 1996;Hnat et al., 2002;Pagel & Balogh, 2002). One method of exploring these excursions is to set a threshold value of the variable and study for how long and by how much the variable continuously exceeds the threshold (Takalo et al., 1999). By studying the statistics of properties of these bursts, it is possible to probe the dynamical process underlying the time series (Sornette, 2000) and to explore the connection between the solar corona, the solar wind, and the Earth's magnetospheric response (Freeman et al., 2000a;Lepreti et al., 2004;Uritsky et al., 2001). As continuous periods of raised activity in the solar wind lead to an increased occurrence of geomagnetic storms (Gonzalez et al., 1994), understanding the statistics of solar wind bursts and their variation over time is relevant for space weather (Schwenn, 2006a).
Statistical distributions of bursts in space plasma parameters are typically fat-tailed and can often be described by power laws over at least part of their range. These include solar flares (Baiesi et al., 2006;Boffetta et al., 1999), solar wind parameters such as Poynting flux, vB s and (Freeman et al., 2000b;Wanliss & Weygand, 2007), and geomagnetic indices (Consolini & De Michelis, 2002;Freeman et al., 2000a). As well as the probability density function (PDF) observations, Uritsky et al. (2001) found a power law dependence between the size and duration of bursts in AE; similar power laws, with different scaling exponents, have since been found by Wanliss and Uritsky (2010) for SYM-H and by Moloney and Davidsen (2011) for the parameter for extremely burst size S, defined as the integrated area between the time series and the threshold; and waiting time w , denoting the time difference between the start times of consecutive valid bursts. The small excursion at around 60 min is excluded from the burst sample, as it is shorter than the imposed 5-min minimum required burst duration. high burst thresholds. These observations, and theoretical considerations, led to a large body of work exploring self-organized criticality (Bak et al., 1987;Watkins et al., 2016) and turbulence as possible mechanisms to produce the ubiquitous scale-free nature of bursts (Aschwanden, 2011;Boffetta et al., 1999;Chapman et al., 1998;Crosby et al., 1992).
Over the 11-year solar cycle, the topology of the solar corona varies extensively, as the prevalence of different coronal structures changes with solar cycle phase (Luhmann et al., 2002;Schwenn, 2006b). This in turn leads to variation in the rate of flares and coronal mass ejections (Cremades & St.Cyr, 2006;Wheatland, 2000), as well as the smaller-scale streams and transients that constitute the solar wind (Behannon et al., 1989;Wang et al., 2000). The solar cycle variation is then mapped via the solar wind into space weather activity in Earth's magnetosphere, so that the observed statistical distributions of both solar wind (Veselovsky et al., 2010) and geomagnetic (Campbell, 1979) variables vary over the solar cycle. The solar cycle evolution of the statistics of bursts is less well explored; when studying the distribution of waiting times between periods of negative B z , D' Amicis et al. (2006) determined that the PDF decayed as a power law at both solar maximum and minimum; however, at minimum there was an exponential cutoff at high waiting times. Hush et al. (2015) found that the survival distribution of the largest AE burst sizes is exponential, with moments specific to each solar cycle phase. Finally, when looking at bursts in solar wind and B 2 over extremely high thresholds, Moloney and Davidsen (2014) reported that the power law exponents of the PDFs of burst duration and size are constant over the solar cycle and that at these high thresholds burst size scales with duration as a power law with a constant exponent.
In this letter, we explore how the scaling properties and distributions of burst parameters in the solar wind B 2 time series vary over solar cycles 23 and 24, and what constraints this may provide for the expected size and duration of large events in the solar wind. We find that the size and duration of bursts in solar wind B 2 are related via a power law with an exponent that is invariant over solar cycles 23 and 24. However, the statistical distributions of burst size and duration do vary both between solar maximum and minimum and from one cycle to the next. We begin by defining the data set and relevant burst parameters in section 2 before presenting our results in section 3. A discussion of the results in the context of established solar wind results, such as the lognormal distribution of the IMF (Burlaga & Ness, 1998) and the 1/f region in the IMF power spectrum (Matthaeus & Goldstein, 1986), is given in section 4. Finally, we conclude in section 5.

Data and Constructing the Samples of Bursts
We use magnetic field and plasma data from the Magnetic Field Investigation (MFI) (Lepping et al., 1995) and Solar Wind Experiment (SWE) instruments (Ogilvie et al., 1995), respectively, on the NASA's Wind satellite (Russell, 1995), sourced from the OMNI data collection. We take data in 1-year intervals surrounding the minima and maxima of solar cycles 23 and 24, using dates December 1995to November 1996, October 1999to September 2000, August 2007to July 2008, and November 2013to October 2014(see Hush et al., 2015. We take data only at times where Wind was in the upstream solar wind. The data has 1-min resolution and, prior to our access, has been shifted to the bow shock nose of the Earth's magnetosphere (King & Papitashvili, 2005).
The key parameters of a burst are shown in Figure 1 using a 2.5-hr interval of the B 2 time series. The duration of the burst is denoted by d , and the burst size, S, is calculated by numerically integrating the difference between the signal value and the threshold over the burst duration. The final parameter is the burst waiting time, which we define as the time between subsequent upward threshold crossings. We exclude any burst with a duration of less than 5 min; as in the 1-min cadence time series they are not sufficiently well resolved. As the threshold is raised, both the number of bursts and their parameters may vary: bursts reduce in size or split up into smaller bursts. Scatter plots of burst size, S, against burst duration, d , using logarithmic axes, for bursts in B 2 during each of the minima and maxima of solar cycles 23 and 24. Bursts are thresholded at the 0.75 (red), 0.85 (green), and 0.95 (purple) quantiles of each data set. The gradient of the best linear regression of log 10 (S) onto log 10 ( d ) for each set is denoted by , and the solid black line shows this fit for bursts over the 0.85 quantile (green) for each solar cycle phase. The uncertainty on is the 95% confidence interval.
As the burst analysis requires a continuous time series, it is necessary to consider the impact of missing data in the Wind observation set. This accounts for around 10% of data in each yearlong time series that we consider here. There are two distinct types of missing data: short duration data gaps which occur randomly and are due to instrumental or sampling constraints and longer duration intervals where Wind crossed into Earth's magnetosphere; the latter are excluded from our analysis. Three methods were trialed for interpolating through the random gaps: using the most recent valid observation throughout the gap, linearly interpolating between the previous and next valid observations, and using the data immediately surrounding the gap to produce a first-order autoregressive (AR(1)) model (Box et al., 2008), seeded at the most recent valid observation. The results were insensitive to the filling method used; throughout this article the random gaps are filled using the linear interpolation method. The magnetosphere crossings occurred only during the minimum and maximum of solar cycle 23. Within these two yearlong time series, we identified the continuous intervals between the magnetospheric crossing gaps and extracted the bursts from each interval, then collected these bursts together to give the total burst sample for each year.

Results
We will discuss in detail results for bursts in the time series of B 2 ; the results for bursts in Poynting flux (estimated as vB 2 ∕ 0 ), proton density, and temperature are qualitatively similar. For each of the yearlong time series spanning the minima and maxima of solar cycles 23 and 24, we extract bursts over a range of thresholds. These thresholds are defined by the quantiles of the time series at each solar cycle phase, where the qth quantile is the value of B 2 which exceeds a proportion q of the data set. We threshold each time series at quantiles in the range 0.75 to 0.95, however, these can correspond to quite different B 2 values; for example, the 0.85 quantile in the B 2 time series during the minimum of cycle 23 is 47.7 (nT) 2 , while the same quantile in the cycle 23 maximum time series is 90.6 (nT) 2 . By using a threshold which is specific to each data set, the overall activity level of each cycle phase is accounted for in the analysis, so that we can compare the likelihood of bursts that are equally extreme relative to the activity of their cycle phase. For the q th = 0.75 threshold there are roughly 1,500 bursts in each year; raising the threshold to the q th = 0.95 quantile reduces the sample size to 400 bursts.
A power law scaling relation S ∼ d has previously been found by Uritsky et al. (2001), Wanliss and Uritsky (2010), and Moloney and Davidsen (2011) studying bursts in AE, SYM-H, and solar wind , respectively, during single phases of the solar cycle. In addition, Moloney and Davidsen (2014) found that for bursts in over extremely high thresholds, the scaling exponent did not change between periods of higher and lower activities. In Figure 2 we plot the burst size against the burst duration for bursts in the B 2 time series during each maxima and minima of solar cycles 23 and 24, using logarithmic axes. We overlay the burst samples gathered using three thresholds, the 0.75 (red), 0.85 (green), and 0.95 (purple) quantiles of each yearlong time series. By regressing log 10 (S) onto log 10 ( d ) using least squares, we obtain an estimate for the power law scaling exponent during each solar cycle phase and for each threshold. Each exponent is shown in the figure alongside its 95% confidence interval. Similar plots for bursts in the Poynting flux, proton density, and proton temperature time series are shown in Figures 3, 4, and 5 respectively. Figure 2 shows that we find the S ∼ d relationship for B 2 and that the exponent is remarkably robust; within uncertainty it does not vary either between solar maximum and minimum or between cycles 23 and 24, which are notably at very different levels of activity. Exponent is also independent of the threshold used to define the bursts; as the threshold is raised, the full scaling relation is shifted upward on the plot, increasing the intercept but not affecting the gradient of the least squares regression. The value of the exponent was found to be the same for time series of varying lengths; however, it slightly increases if we include bursts with durations under 5 min. Solar cycle invariant scaling exponents were also found for bursts in the Poynting flux ( Figure 3), proton density ( Figure 4), and temperature ( Figure 5); however, the exponent itself was different for each variable, ranging from ∼ 1.45 for the temperature to ∼ 1.65 for Poynting flux. For all four variables during each of the four solar cycle phases considered here, the regression coefficient R 2 resulting from the regression of log 10 (S) onto log 10 ( d ) exceeded 0.9; thus, more than 90% of the variability in burst size is described by its relation with burst duration.
The PDFs for bursts in B 2 over the 0.85 quantile threshold are shown on log-log axes in Figure 6. Each PDF is estimated using the kernel density method (Bowman & Azzalini, 1997), with 100 logarithmically spaced bins. For all three burst parameters, the distribution resembles a truncated power law, as found by other studies (Freeman et al., 2000a;Wanliss & Weygand, 2007), and within the power law region the PDF appears the same for all solar cycle phases. However, discerning the change in distribution between solar cycle phases by visual inspection of the PDFs is not optimal; on logarithmic axes the tail of the distribution dominates the plot, so that changes to the lower regions (where the majority of the data lie) are suppressed, while changes in the tail (where there are reduced statistics) are highlighted. The data-data QQ plots (Braun & Murdoch, 2016;Tindale & Chapman, 2016) shown in Figure 7 provide more sensitive tests of how the PDFs and their moments vary over the solar cycle.
The data-data QQ plot is a model-independent graphical diagnostic which tests whether two data sets are drawn from statistical distributions of the same functional form and if so whether the moments of the distribution differ for the two data sets. The QQ plot is constructed using the inverted cumulative distribution functions (CDFs) of the two data sets, where the CDF C(X ≤ x) is defined as the probability that variable X is less than or equal to value x. The CDFs constructed from the two data sets are inverted at the same set of probabilities between 0 and 1, and the QQ plot is realized by plotting the inverted distributions against each other with the probability as a parametric coordinate. Differences between the two CDFs are summed across their ranges, so the QQ plot can provide a more sensitive measure of the variation across the full distribution. A straight line on the plot (or a piecewise linear form) signifies that the distribution functional form is the same  for both samples, and for any standardizable two-parameter distribution model the gradient and intercept of this line capture the difference in the variance and mean of the distribution, respectively, between the two samples (Gilchrist, 2000). The empirical CDFs of the two samples, on which the QQ plot is based, can either be estimated directly using a kernel density estimation or calculated indirectly by constructing the empirical CDF and inverting it at the desired quantiles. Here we use the kernel density method to construct the QQ plot. We take the magnitude difference of these two inverse CDFs as an estimate of the uncertainty. Figure 7 shows comparisons of the distributions of (a) burst duration d , (b) burst size S, and (c) waiting time w , between (i) the minima of cycles 23 and 24, (ii) the maxima of cycles 23 and 24, (iii) the maximum and minimum of cycle 23, and (iv) the maximum and minimum of cycle 24. Results are shown for bursts constructed using five different thresholds between the 0.75 and 0.95 quantiles. The inset in each panel is a zoom of the same QQ plot, showing values that lie within the ranges 0 < d , w < 180 min and 0 < S < 5,400 (nT) 2 min. This 3-hr timescale covers turbulence-scale fluctuations and small-scale transients (Viall et al., 2010) but excludes the largest 10% of events. This timescale combined with the median exceedence of B 2 over the 0.85 quantile (around 30 (nT) 2 ) corresponds to a burst size of 5,400 (nT) 2 min. Because the tails of the burst parameter distributions for solar wind variables decay as a truncated power law, the largest 10% of bursts dominate the main panel QQ plot when displayed with linear axes.
Focusing first on the variation in the distribution of the largest burst durations (Figure 7, row a), we can see that the QQ traces are the same for all five thresholds within uncertainty. For the large events ( d ≳ 100 min), the traces are roughly linear and deviate from the y = x line, most significantly for the two comparisons of cycle maximum to minimum (Figures 7aiii and 7aiv). This suggests that the functional form of the distribution of large burst durations is the same throughout both solar cycles, though the moments of the distribution vary; in particular, the steep gradient indicates a change in the variance of the distribution between solar maximum and minimum. On examining the burst duration PDF (Figure 6i), this region corresponds to the power law roll-off, which is different for each cycle phase. The inset plots show that the smaller bursts undergo a separate, more subtle transformation, which appears either as a straight line with a different gradient (e.g., Figure 7aii inset) or as a slightly curved line (e.g., Figure 7aiii). The QQ plots for burst size in Figure 7 row b are  also linear for the large bursts, again signifying a common distribution form scaled with varying moments. Notably in Figure 7bii, unlike for burst duration, the distribution of burst sizes undergoes a significant change in variance between the maximum of cycle 23 and that of cycle 24. The likelihood of a burst of a given duration was therefore the same at both maxima; however, the reduced activity during the maximum of solar cycle 24 (McComas et al., 2013) made large excursions above the threshold less likely, resulting in smaller integrated burst sizes.
The power law distributions found previously for bursts in the time series of solar wind variables, as well as the power law scaling relation S ∼ d shown here, indicate that the solar wind time series may be considered scale free within some range. For a truly scale-free time series, the statistical character of the bursts should be independent of the threshold, as raising the threshold simultaneously causes the smallest bursts to drop out of the sample and the largest bursts to be split into smaller bursts, in a manner that leaves the overall distribution unchanged. For bursts in B 2 , within the range of thresholds q th = 0.75 to q th = 0.95, we observe a region with little dependence on the threshold for the burst duration d and waiting time w (Figure 7, rows a and c). The traces overlap up to a timescale of around 60 min before slightly diverging, more clearly in the waiting time than burst duration. Conversely, the QQ traces for burst size are threshold dependent across the full range, so the change in the moments of the burst size distribution between solar cycle phases is also dependent on the threshold used to define the bursts.
Lastly, the shape of the waiting time QQ plot (Figure 7, row c) is nonlinear, particularly at high thresholds. The waiting time PDFs (Figure 6iii) show a second peak in the waiting time distribution, which is most prominent in the distributions at solar minimum. This could be due to the large-scale structure of the solar wind being dominated by corotating streams during solar minimum, as opposed to the transients of solar origin which Journal of Geophysical Research: Space Physics 10.1029/2018JA025740 are more prevalent at maximum (Tsurutani et al., 2006). Moloney and Davidsen (2014) also found the waiting time distribution for bursts in solar wind (using high thresholds) and exhibited a shoulder whose position was dependent on the solar cycle phase.
When considered together, Figures 2 and 7 show that while the moments of the statistical distributions of burst duration and burst size independently vary with solar activity, they do so while maintaining a power law scaling relationship between them with a constant scaling exponent. These findings may appear at first to be contradictory; however, they can be reconciled when we consider how changes in variance of the burst duration and size distributions impact upon the power law scaling relation. A change in variance can be expressed as a linear change in variables from S → S S and d → d , where S and are the gradients of the relevant QQ plots. These changes may be due to solar cycle effects, the threshold used to define the bursts, or a combination of both. If the original scaling relation is expressed as then after the transformation of S and d , it becomes log 10 ( S S) = log 10 ( d ) + c.
This simply rearranges to log 10 (S) = log 10 ( d ) + [c + log 10 ( ) − log 10 ( S )] so that the scaling exponent, , remains the same; the change in distribution of the burst parameters will instead be seen in a translation of the intercept. Therefore, there exists a constant characteristic duration for a burst of a given size; however, the likelihood of occurrence of such a burst will change throughout the solar cycle as the moments of the burst size and burst duration distributions vary.
This result may then provide a relevant space weather constraint, as it implies that a characteristic duration exists for a burst of a given size in the solar wind. If we consider that a burst has a size which is proportional to a mean amplitude, A, and has duration, d , then therefore, Thus, if the observed value of B 2 increases to amplitude A above the threshold, it may be expected that the duration of the excursion will be of the order −1 d , where is independent of solar cycle phase, and takes a different value for each solar wind parameter. The scatter around the power law shown in Figure 2 also provides some estimate of the uncertainty in this characteristic duration. Energetic conditions related to various large-scale solar wind structures increase the probability of occurrence of an extreme geomagnetic event (Richardson et al., 2000), while smaller-scale fluctuations caused by turbulence also impact space weather through their effect on the transport of energetic particles (Alouani-Bibi & Le Roux, 2014). The existence of a characteristic duration for a burst of given size, for the full range of events and at any solar cycle phase, could have a potential application in the design of preventative measures to minimize risk from space weather events, in particular since we find that it holds for a broad range of event sizes and at all solar cycle phases.

Implications for Time Series Power Spectral Density and PDF
The bursty nature of time series of solar wind observations has been attributed to the existence of intermittent turbulence within the solar wind plasma (see Bruno & Carbone, 2013, and references therein). The search for other signatures of intermittent turbulence in the solar wind led to a large body of work on the distribution of observed values of the IMF and other solar wind variables (Padhye et al., 2001;Veselovsky et al., 2010), as well as their power spectra (Chapman & Nicol, 2009;Podesta et al., 2007), structure functions (Horbury et al., 1995;Marsch & Liu, 1993), and the distribution of the differences time series (Hnat et al., 2002;Sorriso-Valvo et al., 1999). In particular, the power spectrum of coordinate components of the IMF shows multiple power law regions, with a 1/f region at low frequencies, attributed to events of solar origin (Matthaeus & Goldstein, 1986) and an approximately five-thirds slope at higher frequencies that is attributed to Alfvénic turbulence In addition, the middle panels (aii, bii, cii, and dii) show the PSD of an unadjusted random phase surrogate (green), while the lower panels (aiii, biii, ciii, and diii) show the PSD of an amplitude-adjusted random phase surrogate (pink). Each original yearlong time series is split into overlapping 8-week intervals; the PSD is estimated for each interval, and the final PSD is the average of these estimations. (Kiyani et al., 2015, and references therein). However, the IMF magnitude does not exhibit a similar break in its spectrum (Hnat et al., 2011). The distribution of the IMF magnitude, as well as other solar wind variables, is often approximated by the lognormal model, though the goodness of fit of this model has been found to vary over the solar cycle, and in particular the tail of the lognormal consistently underestimates the likelihood of large values (Tindale & Chapman, 2017). Importantly, however, the time series is not completely characterized by the power spectral density (PSD) (i.e., the decomposition of the second-order moment by frequency) or the PDF of observed values. In addition to these characteristics, the phase information is also reflected in the properties of the burst series that we study here. We can use surrogate time series to study the contribution of each of these traits to the burst behavior. In the presence of multifractality we would expect higher-order moments to carry information (Sornette, 2000), although these may in practice be very difficult to measure accurately (Dudok de Wit & Krasnosel'skikh, 1996) and will not be considered further here.
First, by simply shuffling the order of occurrence of observations in the time series, the amplitude power spectrum can be altered (whitened) while exactly conserving the distribution of the original values. In Figure 8, the gray traces show the flat spectrum of the shuffled time series, relative to the power law spectra (black traces) of the original time series from the minimum of solar cycle 24. Following this shuffling, there are no bursts with extended lengths: the maximum observed burst duration is 7 min. The existence of longer bursts may therefore be attributable to the scaling property of the magnitude power spectrum, that is, its approximate f − form.
Next, we can use random phase surrogates of the time series (Dolan & Spano, 2001) to vary the PDF and the correlation between Fourier phases without affecting the amplitude power spectrum. We will consider two The original PDF is shown by the black line, the distribution of an amplitude adjusted random phase surrogate is shown by the pink dashed line, and the distribution of an unadjusted random phase surrogate is green. Each PDF is calculated using kernel density estimation, with 100 linearly spaced bins. types of random phase surrogate (RPS): the unadjusted RPS (Theiler et al., 1992) and the amplitude-adjusted Fourier transform (AAFT) surrogate (Schreiber & Schmitz, 1996). To create an unadjusted random phase surrogate, we take the Fourier transform of the time series, randomly generate a phase for each frequency, then pair this with the original Fourier amplitude and take the inverse transform to return to the time domain (Kaplan & Glass, 2000). An additional subtlety arises in the calculation of the unadjusted RPS for solar wind time series: in order for the surrogate time series to be real-valued, the Fourier spectrum must be antisymmetric; however, meeting this condition results in a surrogate with both positive and negative values. As the four variables we consider here are positive-defined, we use the magnitude of the unadjusted random phase surrogate as the proxy for each solar wind time series. To create an amplitude-adjusted random phase surrogate, we iterate over the procedure for the unadjusted RPS; at each iteration, we slightly alter the Fourier amplitudes so that the PDF of the surrogate time series is a closer approximation of the original time series PDF. After many iterations, the AAFT surrogate has a PSD and PDF that closely resemble those of the original time series, while the correlation between Fourier phases remains randomized.
The power spectra of unadjusted random phase surrogates for each of the four solar wind variables are shown by the green traces in Figure 8; these are almost identical to the original power spectra (black), with slight discrepancies at low frequencies due to the use of the magnitude of the unadjusted random phase surrogate here. Comparing the PDF of the unadjusted random phase surrogate (green) to the original time series (black) in Figure 9, we see that the distribution has changed its shape from roughly a lognormal to closer to a Gaussian. Again, as we are using the magnitude of the unadjusted RPS, we see only the positive-valued side of the Gaussian. The impact of this change in distribution is visible in the histograms in Figure 10; each histogram shows the distribution of values of the joint area-duration scaling exponent extracted from 500 unadjusted random phase surrogates for each time series. Because Wind made numerous magnetosphere crossings during both phases of solar cycle 23, we take the Fourier transform of the longest continuous interval between crossings in these years, which in both cases spanned roughly 105,000 observations (compared to 525,000 observations in a full year). The exponent is calculated using bursts over the 0.85 quantile of each surrogate time series. The burst size versus duration scaling exponents from the original time series (as displayed in Figures 2-5) are shown by the relevantly colored cross, along with their 95% confidence interval.
For B 2 and the Poynting flux, during all four solar cycle phases, the scaling exponent is reduced by roughly 7% following the random phase shuffling. The exponent for the proton density and temperature during the cycle 23 maximum (yellow) and the two phases of solar cycle 24 (green and blue) is also reduced but by a lesser amount (∼ 4%). However, during the minimum of cycle 23 (red) the confidence interval of the original scaling exponent encompasses the histogram of surrogate values. This may be related to the reduced interval of observations from which the surrogates were created for the cycle 23 minimum, which for the density and temperature may not be a true reflection of the full yearlong interval used to compute the original burst scaling exponent. The R 2 regression coefficient for the linear fit of log 10 (S) against log 10 ( d ) exceeds 0.9 for all surrogates, indicating that the variation in burst sizes is equally well described by the relation S ∼ d before and after the random phase shuffling. The heavy-tailed nature of the underlying distribution therefore impacts upon the value of the scaling exponent, with a heavier tail resulting in a higher exponent. The robustness of the observed scaling exponent to the solar cycle may then be related to the roughly constant width of the lognormal model used to approximate the distribution of solar wind variables (Burlaga, 2001;Tindale & Chapman, 2017).
We can extract the impact of varying the phase only by using AAFT random phase shuffling (Schreiber & Schmitz, 1996). As described above, surrogate time series computed using this method approximately retain Figure 11. Histograms of the burst size versus burst duration scaling exponent, , in the same format as Figure 10. Here each surrogate is calculated using the amplitude adjusted Fourier transform method.
both the PSD and PDF of the original time series, while the Fourier phases are randomized. Here we have computed AAFT surrogates using the algorithm given in Kaplan and Glass (2000), as implemented in MAT-LAB by Silva and Moody (2014). In Figures 8 and 9, the pink traces show the power spectra and distribution, respectively, of one AAFT random phase surrogate, again based on data from the minimum of solar cycle 24. The power spectra of the original time series (black), the unadjusted random phase surrogate (green), and the AAFT surrogate (pink) are sufficiently close that to allow a clear comparison between the PSD of each surrogate and the original time series, we are required to plot the panel for each solar wind variable 3 times in Figure 8.
Five hundred AAFT surrogates are created from each of the four cycle phase time series of the four solar wind variables. For the minimum and maximum of solar cycle 23, we again use the reduced, continuous intervals between magnetosphere crossings. Histograms of values extracted from the AAFT surrogates are shown in Figure 11, in the same format as Figure 10. Now that both the PSD and PDF of each surrogate are preserved from the original time series, the joint area-duration scaling exponent is less strongly modified. In particular, for B 2 and the Poynting flux (Figures 11a and 11b), each histogram of the values of the joint area-duration scaling exponent lies within the confidence interval of its value for the original time series. However, for the proton density and temperature (Figures 11c and 11d), a more complex picture emerges. The surrogate histograms are well aligned with the original value in some cases, for example, density during the maximum of cycle 23 (yellow) and temperature at the minimum of cycle 24 (green), while in other cases the surrogate values are shifted outside the confidence intervals of the original, such as density during the minimum of cycle 24 (green). This suggests that the contribution of the phase information in a time series to the properties of burst is unique to both the variable studied and the period of observation.
While the properties of the bursts extracted from a time series are inextricably linked to characteristics such as its power spectrum and the probability distribution of its values, it is clear that additional information about the underlying process is encoded into the bursts. In particular, despite having different power spectra and PDFs, the four solar wind variables are united by the scaling properties of their bursts, where for all four variables the size and duration of a burst are related by a power law whose exponent is approximately constant over the solar cycle. The exponent is, however, unique to each solar wind variable. The robustness of the exponents against solar cycle changes suggests that they are the signature of some underlying physical processes.
This may on the large scale relate to the physics of the solar corona, where power law burst distributions are seen in flare statistics (e.g., Wheatland et al., 1998) and on the small scale are related to in situ turbulence (Kiyani et al., 2007). However, an alternative interpretation may be that of similarity in plasma structures carried by the solar wind (Nieves-Chinchilla et al., 2018). If, for example, these bursts are dominated by force free flux ropes, then such a structure seen on a given length and timescale can be rescaled to structures seen at all other length and timescales. The observing satellite then sees a burst as it flies through these structures at an arbitrary angle. To confirm this interpretation requires detailed modeling of the solar wind, but our results could offer an additional check on such modeling.
An open question is how the scaling exponents of the burst size versus duration plots relate to underlying physics. By constructing surrogate time series, we have shown that the exponent depends on the form of the PSD, on the non-Gaussian nature of the PDF, and on additional information contained in correlations between the Fourier phases. This is consistent with models of intermittent turbulence in that Kolmogorov's original scaling exponent of the energy spectrum is necessary but not sufficient to specify the time series. A lognormal probability density for the magnetic field magnitudes also does not uniquely specify that the component time series are intermittent (Hartlep et al., 2000); in addition, information contained in correlations between the Fourier phases and in higher-order scaling exponents (the intermittency parameter) are needed (Frisch, 1995).
The behavior of bursts (also known as excursions; (Majumdar & Comtet, 2005)) has been studied for only a few idealized time series models (Sornette, 2000) which are mathematically tractable. These may also be compared to our observational findings. Of particular note is Brownian motion, which results from summing temporally uncorrelated white Gaussian noise. Here the burst size S scales as T 3∕2 , which can be intuitively understood by the way in which a burst integrates the √ T scaling of the characteristic amplitude of Brownian motion with respect to time. Fractional Brownian motion (fBm; Sornette, 2000) allows the incorporation of long range memory in the increments and thus changes the self-similarity parameter H from one half to lie in the range from 0 to 1. In consequence Watkins et al. (2012) noted that the joint burst area-duration scaling should generalize to S ∼ T 1+H . Numerical simulations of fBm Watkins et al., 2009) support this view. However, these are all nonstationary models, whereas long-time stationarity might be expected in the solar wind on physical grounds. For stationary models the mean is constant and so the above reasoning would lead to S ∼ T.

Conclusion
Continuous periods of high activity in the solar wind can increase the likelihood of occurrence of severe space weather events; thus, it is important to understand the frequency and magnitude of such periods. By studying the bursts in a time series, where a burst is defined as a continuous period over which a variable exceeds a threshold, it is possible to characterize the probability of occurrence of a large event of a given size or duration. This method has also been used in the past to study the nature of the dynamical process that underlies the observed time series of solar wind variables. Here we specifically study how the statistics of bursts in the solar wind vary over the solar cycle and from one solar cycle to the next. Our results provide a constraint on the likelihood and properties of a burst in the time series of solar wind magnetic energy density, B 2 , and solar wind Poynting flux, proton density, and temperature.
We compared the samples of bursts extracted from B 2 time series centered around the maxima and minima of solar cycles 23 and 24 and found a power law scaling relationship, S ∼ d , between the integrated size and the duration of bursts that holds for both the maximum and minimum of both solar cycles. We find that the joint area-duration scaling exponent is approximately unchanged both between solar maximum and minimum and between the maxima and minima of cycles 23 and 24, for bursts defined over a wide range of thresholds. This provides an important constraint: a burst of a given duration will have a characteristic integrated size, and equally a given amplitude of the B 2 time series will persist for a typical duration, independent of the solar cycle or cycle phase during which the burst occurs.
We then considered the statistical distributions of burst size and burst duration independently. For the large values of burst size, burst duration and waiting time, we find that the functional form of their distribution is constant throughout the two solar cycles; however, there is a significant change in the moments of each distribution with solar cycle phase. The variance of both distributions is larger at solar maximum than at minimum, with a more dramatic change within cycle 23 than 24. We also find when comparing the 10.1029/2018JA025740 distributions between the two cycle maxima that the moments of the burst size distribution change considerably while the burst duration distribution remains approximately unchanged. This indicates that a burst of a given duration was equally likely to occur in both maxima, but the lower solar activity during the maximum of cycle 24 reduced its probable size.
Finally, we discussed the interpretation of these results in the light of the established characteristics of solar wind time series, including the power-law spectrum and heavy-tailed distribution of the IMF. Through the use of random phase surrogates, we saw that while correlation between values in the time series allows the existence of long bursts, in our case it was the non-Gaussian nature of the PDF that dominated the scaling of burst size and burst duration. However, the relationship between the burst properties, power spectra, and PDFs is nontrivial, and importantly additional information such as that contained in correlations between the Fourier phases (as opposed to the amplitude spectrum) can be essential to characterizing the scaling properties of solar wind bursts.
We find that while the likelihood of occurrence of a solar wind burst of a given magnitude varies both within and between the different solar cycles, it has a corresponding (magnitude dependent) characteristic timescale that is independent of solar cycle phase and peak activity. The joint area-duration scaling exponent is insensitive to solar cycle changes over a broad range of burst sizes and durations, and this may suggest they share common underlying physics. This could be related to coronal activity, active in situ solar wind turbulence, physical similarity in solar wind structures such as flux ropes of different size, or some mixture of these. Large-scale bursts can have a significant space weather impact, both through increasing the likelihood of extreme events and through modulating cosmic ray propagation. Our results thus may provide an important space weather relevant constraint on the properties of these bursts.