Wavelet-Based Characterization of Small-Scale Solar Emission Features at Low Radio Frequencies
Akshay Suresh, Rohit Sharma, Divya Oberoi, Srijan B. Das, Victor Pankratius, Brian Timar, Colin J. Lonsdale, Judd D. Bowman, Frank Briggs, Roger J. Cappallo, Brian E. Corey, Avinash A. Deshpande, David Emrich, Robert Goeke, Lincoln J. Greenhill, Bryna J. Hazelton, Melanie Johnston-Hollitt, David L. Kaplan, Justin C. Kasper, Eric Kratzenberg, Mervyn J. Lynch, S. Russell McWhirter, Daniel A. Mitchell, Miguel F. Morales, Edward Morgan, Stephen M. Ord, Thiagaraj Prabu, Alan E. E. Rogers, Anish Roshi, N. Udaya Shankar, K. S. Srivani, Ravi Subrahmanyan, Steven J. Tingay, Mark Waterson, Randall B. Wayth, Rachel L. Webster, Alan R. Whitney, Andrew Williams, Christopher L. Williams
aa r X i v : . [ a s t r o - ph . S R ] J un Draft version September 19, 2018
Preprint typeset using L A TEX style AASTeX6 v. 1.0
WAVELET-BASED CHARACTERIZATION OF SMALL-SCALE SOLAR EMISSION FEATURES AT LOWRADIO FREQUENCIES
A. Suresh , R. Sharma , D. Oberoi , S. B. Das , V. Pankratius , B. Timar , C. J. Lonsdale , J. D. Bowman ,F. Briggs , R. J. Cappallo , B. E. Corey , A. A. Deshpande , D. Emrich , R. Goeke , L. J. Greenhill ,B. J. Hazelton , M. Johnston-Hollitt , D. L. Kaplan , J. C. Kasper , E. Kratzenberg , M. J. Lynch ,S. R. McWhirter , D. A. Mitchell , , M. F. Morales , E. Morgan , S. M. Ord , , T. Prabu ,A. E. E. Rogers , A. Roshi , N. Udaya Shankar , K. S. Srivani , R. Subrahmanyan , , S. J. Tingay , , ,M. Waterson , R. B. Wayth , , R. L. Webster , , A. R. Whitney , A. Williams , C. L. Williams Indian Institute of Science Education and Research, Pune-411008, India; [email protected] National Centre for Radio Astrophysics, Tata Institute for Fundamental Research, Pune 411007, India Indian Institute of Science Education and Research, Kolkata-741249, India MIT Haystack Observatory, Westford, MA 01886, USA California Institute of Technology, Pasadena, CA 91125, USA School of Earth and Space Exploration, Arizona State University, Tempe, AZ 85287, USA Research School of Astronomy and Astrophysics, Australian National University, Canberra, ACT 2611, Australia Raman Research Institute, Bangalore 560080, India International Centre for Radio Astronomy Research, Curtin University, Bentley, WA 6102, Australia Kavli Institute for Astrophysics and Space Research, Massachusetts Institute of Technology, Cambridge, MA 02139, USA Harvard-Smithsonian Center for Astrophysics, Cambridge, MA 02138, USA Department of Physics, University of Washington, Seattle, WA 98195, USA School of Chemical & Physical Sciences, Victoria University of Wellington, PO Box 600, Wellington 6140, New Zealand Department of Physics, University of Wisconsin–Milwaukee, Milwaukee, WI 53201, USA Department of Atmospheric, Oceanic and Space Sciences, University of Michigan, Ann Arbor, MI 48109, USA CSIRO Astronomy and Space Science (CASS), PO Box 76, Epping, NSW 1710, Australia ARC Centre of Excellence for All-sky Astrophysics (CAASTRO) National Radio Astronomy Observatory, Charlottesville and Greenbank, USA Instituto di Radioastronomia, Instituto Nazionale di AstroFisica, Bologna, Italy SKA Organisation, Macclesfield SK11 9DL, UK School of Physics, The University of Melbourne, Parkville, VIC 3010, Australia
ABSTRACTLow radio frequency solar observations using the Murchison Widefield Array have recently revealedthe presence of numerous weak, short-lived and narrow-band emission features, even during moder-ately quiet solar conditions. These non-thermal features occur at rates of many thousands per hourin the 30 .
72 MHz observing bandwidth, and hence, necessarily require an automated approach fortheir detection and characterization. Here, we employ continuous wavelet transform using a motherRicker wavelet for feature detection from the dynamic spectrum. We establish the efficacy of thisapproach and present the first statistically robust characterization of the properties of these features.In particular, we examine distributions of their peak flux densities, spectral spans, temporal spans andpeak frequencies. We can reliably detect features weaker than 1 SFU, making them, to the best of ourknowledge, the weakest bursts reported in literature. The distribution of their peak flux densities fol-lows a power law with an index of -2.23 in the 12 −
155 SFU range, implying that they can provide anenergetically significant contribution to coronal and chromospheric heating. These features typicallylast for 1 − − Keywords:
Sun: corona — Sun: radio radiation INTRODUCTIONThe new generation radio arrays are revealing thepresence of previously unappreciated variety and com- plexity in non-thermal solar emission features at lowradio frequencies (Oberoi et al. 2011; Morosan et al.2015; Tun Beltran et al. 2015). The observations from
Suresh et al. the Murchison Widefield Array (MWA) reveal thepresence of numerous short-lived, narrow-band emissionfeatures, even during what are conventionally regardedas moderate and quiet solar conditions. In terms ofmorphology in the MWA dynamic spectra (DS), thesenon-thermal features appear like miniature versions ofsolar type-III bursts, with spectral and temporal spansof about a few MHz and a second, respectively. Earlierradio imaging studies (Oberoi et al. 2011) of suchfeatures have found their brightness temperatures to besimilar to those expected for type-III bursts, implying acoherent emission mechanism behind their production.The seemingly ubiquitous presence of these featuresraises the possibility that they might correspond toobservational signatures of nanoflares. Characterized byenergies in the range of 10 –10 ergs, nanoflares werehypothesized by Parker (1988) as a plausible solution tothe coronal heating problem. At high frequencies (EUVand X-ray), the observable electromagnetic signaturearises from thermal emission due to local heating of theplasma to very high temperatures by nanoflares. At lowradio frequencies, the emission associated with theseenergetic electrons arises from coherent plasma emissionmechanisms, thus, allowing even a low energy eventto give rise to a much larger observational signature.This advantage makes low radio frequencies the bandof choice for investigating signatures of weak coronalenergy release events. In order to contribute effectivelyto coronal and chromospheric heating, the power law( dN/d W ∝ W α ) index, α , of flare energies ( W ) mustsatisfy the condition that α ≤ − α ≤ − α ≈ − − Jy) for type-I bursts.Type-I bursts, also referred to as radio noise storms,generally consist of short-lived ( . .
10 MHz) bursts that usually last for extended peri-ods and are accompanied by an enhanced broadbandcontinuum emission. Spectral and imaging observationsof radio noise storms, performed by Gergely & Kundu(1975) and Duncan (1981), reveal strong similaritiesbetween Type-I and decametric type-III sources. Onthe basis of a survey of 10,000 type-III bursts observedusing the Nan¸cay Radioheliograph, Saint-Hilaire et al.(2013) report a power law with α ≈ − . − SFU)of type-III bursts. However, unlike these type-III burstsand the type-I bursts investigated by Mercier & Trottet(1997), the small-scale features observed in the MWADS are weaker with typical fluxes of about 1-100 SFU.As the presence of such weak features in the MWA solar data has been established (Oberoi et al. 2011)only comparatively recently, their detailed obser-vational characteristics in terms of distributions oftheir spectral and temporal widths, energy content,and slopes in the frequency-time plane are yet to bedetermined. Such a statistical characterization of theproperties of these features would be the first steptowards understanding them and evaluating theircontribution towards solar coronal heating. However,their high occurrence rate of thousands of featuresper hour in the 30.72 MHz bandwidth MWA DSnecessitates an automated approach for their detectionand subsequent parameter extraction from the DS.Here we present a wavelet-based automated techniquefor robust detection and characterization of theseweak features under conditions of quiet to moderatesolar activity. Though the current implementationis tuned for the MWA DS, the technique itself ismore general and can be applied to DS from otherinstruments. As new state-of-the-art observationalfacilities flood the community with unprecedentedlarge volumes of high-quality data, the need for au-tomated data mining and analysis techniques of thesort presented here is only expected to grow more acute.Section 2 of this paper describes the observational ca-pabilities of the MWA and the data selected for sub-sequent analysis. Section 3 details the wavelet-basedapproach for automated feature detection. A statisticalanalysis of the properties of these features is presentedin section 4. The physical significance of the results ob-tained is discussed in section 5. Finally, a summary ofthe results obtained and the conclusions from our studyare presented in section 6 of this article. OBSERVATIONS AND PRE-PROCESSINGThe MWA is a low frequency radio interferometeroperational in the frequency range from 80-300 MHz.It is a precursor to SKA-Low and is located in theradio-quiet environment of the Murchison Radio Ob-servatory in Western Australia. The MWA consistsof 2048 dual-polarization dipoles arranged as 128tiles, wherein each tile is a 4 × ) for this day, this utomated Solar Feature Recognition and Characterization .
36 MHzto 232 .
96 MHz in 5 steps of 30 .
72 MHz, spending 4 min-utes at each frequency band. The entire 30.72 MHzbandwidth in each data set is comprised of 24 coarsespectral channels, each 1 .
28 MHz wide. Each coarsespectral channel is further composed of 32 fine spectralchannels with a resolution of 40 kHz each. The time res-olution of the data collected is 0.5 seconds. The MWAinterferometric data above 100 MHz is flux-calibratedaccording to the prescription developed by Oberoi et al.(2016) and Sharma et al. (2017). This flux calibrationtechnique provides estimates of the solar flux densitiesand brightness temperatures by accounting for knowncontributions from the sky, the receiver, and groundpickup noise to the system temperature. The receivertemperatures and ground pickup temperatures are ob-tained from a mix of laboratory and field measurements.Estimates of the sky temperature are obtained using theHaslam et al. (1982) 408 MHz all-sky map, scaled witha spectral index of 2.55 (Guzm´an et al. 2011), as a skymodel. The need to keep the Sun unresolved for ap-plication of this flux calibration technique constrainsus to using only short baselines. This non-imagingstudy uses data from one such short baseline of phys-ical length 23.7 m between tiles labeled “Tile011MWA”and “Tile021MWA”. The outputs from the flux calibra-tion technique described in Oberoi et al. (2016) form theinputs for our study. Here, we present the analysis fordata collected in the XX polarization alone, that for theYY polarization is analogous. METHODOLOGYFigure 1 depicts a sample raw MWA DS of normal-ized cross-correlations on the left and its flux-calibrationversion on the right. The features of interest in thiswork appear as short-lived, narrow-band vertical streaksagainst a broadband background continuum.3.1.
Removal of instrumental artifacts
The horizontal features are instrumental artifactsarising due to the poor instrumental response at theedges of coarse spectral channels and need to be re-moved. These artifacts are corrected for by performinglinear interpolation across the systematics-affectedchannels. As the coarse channel edges at the very startand end of the observing band cannot be corrected by interpolation, these are simply discarded. Recordingglitches sometimes affect the beginning and end of datarecording. To avoid contamination from such issues,we routinely discard the first six and the last nine timeslices of data as well.Though the MWA is located in a region with verylittle radio frequency interference (RFI), radio waves re-flected from aircraft can occasionally interfere with theradio signals picked up by a tile and thereby, corrupt thedata collected. Manual RFI-flagging followed by linearinterpolation across RFI-affected segments of the DS iscarried out to ensure a RFI-free DS for efficient featuredetection. The left panel of Fig. 2 displays an instru-mental artifact-free version of the DS shown in Fig. 1.3.2.
Background continuum subtraction
The solar radiation can be thought of as a superposi-tion of sporadic non-thermal radio features with a spec-trally varying, broadband background continuum. Spec-tral variations in the background flux density can oftendistort the spectral profiles of features in the DS. For im-proving our efficiency at picking up small-scale featuresfrom the DS, it is, therefore, necessary to disentanglespectral flux density variations arising from these fea-tures from that associated with the background.As the day of our observations was characterized bymedium levels of solar activity, it seems reasonable toexpect that the thermal quiet Sun emission forms thedominant component of the background continuumemission in our data. We find the temporal variationof the background flux density to be negligible over theduration of individual observing scans of four minuteseach. This allows us to then ignore the time dependenceof the background flux density and treat it as as afunction of frequency alone. As the flux densities ofthe weakest radio bursts detected in our data sets areonly a few percent of the background flux, an accurateand robust means of determining and subtracting outthe spectral variation of the background component isrequired.In this work, the Gaussian Mixtures Model (GMM)routine provided by Scikit-Learn (Pedregosa et al.2011) is applied for estimation of the backgroundflux density ( S ⊙ ,B ) as a function of frequency. Asthe background is expected to vary smoothly, the DSis divided into contiguous groups of 4 fine spectralchannels each. The data in each of these groups is thendecomposed as a sum of Gaussians using the GMMroutine. As there exists no unique way of representinga given function as a sum of Gaussians, the BayesianInformation Criterion (Burnham & Anderson 2002) has Suresh et al. ν ( M H z ) N o r m a li s e d c r o ss c o rr e l a t i o n s ν ( M H z ) S ⊙ ( S F U ) Figure 1 . Left: A sample MWA DS of normalized cross-correlations. Right: Flux-calibrated version of the same DS. ν ( M H z ) S ⊙ ( S F U ) ν ( M H z ) S ⊙ − S ⊙ , B ( S F U ) Figure 2 . Left: A version of the flux calibrated DS free of instrumental artifacts. Right: Background-subtracted version of thesame DS. Features of interest can be easily identified in this processed DS. We note that these features also overlap in manyinstances. One feature that appears to be relatively isolated from the others is marked by a red circle in the right panel. been employed to determine the optimum number ofGaussians required to fit the data. Since the thermalquiet Sun component forms the baseline emissionlevel on top of which non-thermal radio emission isdetected, it is reasonable to assume that the Gaussiancorresponding to this background continuum must bethe one with the lowest mean and the highest weight.For every group of fine spectral channels, the value ofthe mean of this Gaussian is noted as the backgroundflux density ( S GMM ( ν )) at the respective frequency andis shown by the red circles in top sub-panels in Fig. 3.Presence of strong, frequent radio bursts that outshinethe background component in a DS degrade the abilityof GMM to determine a value of the background flux ateach observing frequency. Our observations were takenon a day with moderate solar activity, allowing for theuse of GMM to determine the background flux density at several frequencies in most data sets.A degree-4 polynomial is then used to fit the large-scale smooth spectral trend in the background flux den-sity and is subtracted from the DS. The right panel inFig. 2 depicts the DS obtained after background re-moval from the DS depicted in the left panel. The suit-ability of a degree-4 polynomial fit to the backgroundcan be quantified by estimating the residual percent-age between S GMM ( ν ) and the flux densities ( S fit ( ν ))predicted from the best fit polynomial at the same fre-quency. The residual percentage is given by:Residual %( ν ) = S GMM ( ν ) − S fit ( ν ) S GMM ( ν ) × utomated Solar Feature Recognition and Characterization B a c k g r o un d f l u x ( S F U )
110 115 120 125 130 135 140 145 ν (MHz) −4−3−2−10123 F i t R e s i d u a l ( % ) (a) B a c k g r o un d f l u x ( S F U )
140 145 150 155 160 165 170 175 ν (MHz) −2.5−2.0−1.5−1.0−0.50.00.51.01.52.0 F i t R e s i d u a l ( % ) (b) B a c k g r o un d f l u x ( S F U )
170 175 180 185 190 195 200 205 ν (MHz) −4−3−2−10123 F i t R e s i d u a l ( % ) (c) B a c k g r o un d f l u x ( S F U )
200 205 210 215 220 225 230 235 ν (MHz) −5−4−3−2−10123 F i R e s i d u a l ( % ) (d) Figure 3 . Degree-4 polynomial fits to the spectral trend in the background continuum and the residuals to the fits. The 4 panelscorresponds to 4 different data sets with frequency ranges: (a) 110 . − .
54 MHz, (b) 141 . − .
26 MHz, (c) 171 . − .
98 MHz, and (d) 202 . − . study. A degree-4 polynomial is adequate to describethe spectral variation observed in the background fluxto within a mean absolute error of 3-4%.3.3. Wavelet-Based Feature Detection
Continuous wavelet transform (CWT) provides a nat-ural way of obtaining a time-frequency representation ofa non-stationary signal through the use of a wavepacketwith finite oscillation, i.e, a wavelet. In this work, oursignal is the 2D MWA DS containing the features of in-terest. The efficiency of CWT at reliable detection offeatures from the DS depends upon our choice of the 2Dmother wavelet and is maximized for a mother waveletwhich closely matches the shape of the spectral and tem- poral profiles of these features.3.3.1.
Choice of mother wavelet
From the right panel in Fig. 2, it can be seen thatwhile there do exist features that appear isolated in theDS, several features in the DS appear to bunch together.Figure 4 depicts the spectral and temporal profiles of oneseemingly isolated feature indicated by a red circle in thebackground-subtracted DS depicted in Fig. 2. A closelook at such isolated features in the DS reveals a charac-teristic smooth, unimodal nature to their temporal andspectral profiles. Assuming that each atomic featurein a DS possesses unimodal spectral and temporal pro-files, any multi-modal spectral or temporal distributionof flux densities observed can be interpreted as a su-
Suresh et al. perposition of contributions from constituent unimodaldistributions. This allows for a 2D Ricker wavelet to beemployed as a suitable mother wavelet for CWT. Mea-sured in pixel units, the features of interest usually haveaxial ratios of about 10 - 50. To best match featuresof this nature, we use a variable separable version of a2D Ricker (also called the Mexican Hat) wavelet withanalytical form : R ( t, ν ) = 43 √ π (cid:16) (1 − t ) e − t (cid:17) (cid:16) (1 − ν ) e − ν (cid:17) (2)as the mother wavelet. From this mother wavelet,scaled wavelets are constructed according to the defi-nition given by Antoine et al. (2004) as follows: R s t ,s ν ,τ t ,τ ν ( t, ν ) = 1 √ s t s ν R (cid:18) t − τ t s t , ν − τ ν s ν (cid:19) (3)The peak of a 2D scaled Ricker wavelet is located at( t, ν ) = ( τ t , τ ν ). The scales s ν and s t correspond tohalf the support of its positive lobe along the frequencyand time directions respectively. Using the scaled 2DRicker wavelets, wavelet coefficients of the DS are thencomputed according to the following definition: γ ( s t , s ν , τ t , τ ν ) = ¨ ν,t DS ( t, ν ) R s t ,s ν ,τ t ,τ ν ( t, ν )dt d ν (4)For ease of notation, let us denote the wavelet coef-ficients γ ( s t , s ν , τ t , τ ν ) by the symbol γ ( s t , s ν , t, ν ). Fora given feature peaked at ( t, ν ) in the DS, γ ( s t , s ν , t, ν )shall be maximized when s ν and s t match with the spec-tral and temporal extents of the feature respectively.Thus, the 2D Ricker wavelet acts as a peak and supportdetection filter. This then enables us to determine thepeak flux densities as well as the temporal and spectralextents of features in the DS.3.3.2. Construction of a composite matrix
Owing to the fact that the 2D CWT introduces twoadditional degrees of freedom through transformationfrom a 2D DS space to a 4D wavelet-coefficient space, alarge number of wavelet coefficients computed for a givenDS carry redundant information. The non-orthogonalityof a set of scaled Ricker wavelets further preserves thisredundancy. This aspect can then be exploited to re-construct the DS using a basis different from the set ofscaled wavelets. Torrence & Compo (1998) give an ex-plicit expression for 1D signal reconstruction from thewavelet coefficients using a basis of δ − functions. Ex-tending this formula to the 2D CWT used here, a com-posite matrix, A ( t, ν ), of wavelet coefficients that ex-actly reconstructs the DS, barring a constant normal-ization factor, is given by: A ( t, ν ) = X s ν > X s t > γ ( s t , s ν , t, ν ) √ s t s ν (5) As the wavelet coefficients are nothing but a convolutionof the DS with the scaled wavelets, it is expected that A ( t, ν ) should be a smooth reconstruction of the DS.Local maxima in A ( t, ν ) then correspond to peaks offeatures in the actual DS. However, there are two issueswith using A ( t, ν ) for feature identification. At smallscales, our measurements are dominated by noise. AsEq. 5 involves a sum over all values of s ν and s t , it alsotries to incorporate the measurement noise in A ( t, ν ).Further, bunching of features leading to overlapping ofspectral and temporal profiles of adjacent features inthe DS can hinder the ability of CWT to resolve twoclosely spaced features from one another at large valuesof s ν and s t . Hence, it is necessary to work with an in-termediate range of scales for constructing a compositematrix, M ( t, ν ), that captures details of the features ofinterest while avoiding being influenced by the inherentmeasurement noise at small scales and bunching of fea-tures at large scales. M ( t, ν ) is therefore, constructedusing the following expression: M ( t, ν ) = s ν,upper X s ν = s ν,lower s t,upper X s t = s t,lower γ ( s t , s ν , t, ν ) √ s t s ν (6)In the time domain, the features of interest are al-ready present at the resolution of the data, forcing usto set s t,lower to 0.5 seconds. Careful visual inspec-tion of a large number of DS revealed that there ex-ist few features with bandwidths less than 0.5 MHz,leading us to a choice of 0.5 MHz for s ν,lower . Again,guided by careful visual inspection of several DS, weset s t,upper = 3 seconds and s ν,upper = 5 MHz in or-der to provide both the ability to detect atomic featurespresent within a bunch of features and the capability toidentify relatively long-lived or broadband features. Thevalues chosen for s t,upper and s ν,upper in fact enable us toreliably reconstruct features with spectral and temporalextents as large as 26.04 MHz and 15 seconds respec-tively. M ( t, ν ) is then computed using the choices ofscales mentioned above. Local maxima picked up from M ( t, ν ) correspond to locations of the peak flux den-sities of different features contained in the DS. Figure5 illustrates the ability of CWT to distinguish betweenclosely spaced features despite overlaps in their flux den-sity profiles along the frequency and time axes.Panel (a) in Fig. 5 depicts a comparison between aspectral slice taken from both the DS and M ( t, ν ) atthe same time. The location of peaks of features in the M ( t, ν ) spectral profile closely agree with their corre-sponding peaks in the DS. For a given feature, we findthat its spectral extent is matched well by the distancebetween the two local minima in the M ( t, ν ) spectralslice that straddle its peak. We use the distance between utomated Solar Feature Recognition and Characterization
140 145 150 155 160 165 170 175 ν (MHz) −505101520 F l u x d e n s i t y o f f e a t u r e , S ⊙ , F ( S F U )
86 87 88 89 90
Time (seconds) F l u x d e n s i t y o f f e a t u r e , S ⊙ , F ( S F U ) Figure 4 . Left: Spectral profile of the feature marked by a red circle in the background-subtracted DS depicted in Fig. 2. Notethat two other weaker features are present at the upper and lower frequency ends in the time slice corresponding to the peakof this feature. Right: Temporal profile of the same feature. The red circles in the left and the right panels of this figure markthe location of the feature peak, as shown in the right panel of Fig. 2, along the frequency and time axes respectively.
200 205 210 215 220 225 230 235 ν (MH ) −15−10−5051015 F l u x d e n s i t y o f f e a t u r e , S ⊙ , F ( S F U ) Spectral profile from DSSpectral profile from M(t, ν) (a)
Time (seconds) −20020406080100 F l u d e n s i t y o f f e a t u r e , S ⊙ , F ( S F U ) Temporal profile from DSTemporal profile from M(t, ν) (b)
Figure 5 . (a) A spectral profile taken from both the DS (blue) and M ( t, ν ) (green) at the same time. The presence of threelocal maxima in the M ( t, ν ) profile indicates the ability of CWT to successfully detect the three features seen in the DS profile.The locations of the local minima of M ( t, ν ) enable us to distinguish individual features from one another despite overlaps intheir spectral profiles in the DS. (b) Panel illustrating ability of M ( t, ν ) to reproduce temporal widths of features reliably, whileproviding the resolution to distinguish between features located close together. The presence of multiple peaks in M ( t, ν ), onecorresponding to each local maximum in the DS, clearly demonstrates the ability of M ( t, ν ) to detect these features. these local extrema as the spectral extent of the feature.The lower extremum is then taken to be start frequency( ν start ) of the feature. The temporal extent and starttime ( t start ) of a feature are similarly estimated. In or-der to obtain estimates of a quantity similar to the halfpower width of a feature, we define the spectral andtemporal widths of a feature respectively as:∆ ν = 0 . × Spectral extent of feature∆ t = 0 . × Temporal extent of feature. For the purpose of quantifying any symmetry present inthe spectral profile of a feature with peak at frequency ν , we define its spectral symmetry parameter as follows: χ ν = ν − ν start ν (7)The value of this parameter lies in the range from 0 to1. A spectral symmetry parameter value of 0.5 for a fea-ture represents a perfectly symmetric frequency profilewhile departures from 0.5 indicate skewness present inthe spectral profile. The temporal symmetry parameter Suresh et al. ν ( M H z ) S ⊙ − S ⊙ , B ( S F U ) ν ( M H z ) S ⊙ − S ⊙ , B ( S F U ) Figure 6 . Peaks of features detected from the background-subtracted DS depicted in Fig. 2 are depicted as green circles. Theleft and the right panels differ only in the color bar range. While the left panel illustrates the ability of the CWT algorithm topick up bright features, the right panel shows the ability of the CWT code to pick up relatively weaker features as well. ( χ t ) of a feature is similarly defined.3.4. Correction of peaks detected
As seen from panel (a) in Fig. 5, peaks of featurespicked up from M ( t, ν ) do not always coincide with theircounterparts in the DS. However, since M ( t, ν ) peaks lieclose to their corresponding DS peaks, this discrepancyis easily corrected by first growing a region around a M ( t, ν ) peak and then, identifying the DS peak withinthis region. The admissibility criterion used to grow aregion S starting from a M ( t, ν ) peak is that the waveletcoefficient of the neighboring pixel under considerationis within a minimum threshold (T) percentage of thepeak wavelet coefficient. The region growing algorithmterminates when no more pixels on the boundary of Ssatisfy this criterion. Since M ( t, ν ) is only an approxi-mation to the actual DS, the temporal and spectral pro-files of a feature in the DS are reproduced exactly onlywithin a small neighborhood around its M ( t, ν ) peak.Hence, a value of T as high as 95% has been chosen toensure that all pixels contained in the region S aroundthe peak of a feature actually belong to this feature.For the features detected from all DS used in this work, M ( t, ν ) peaks show average offsets of 0.16 seconds and0.57 MHz from their corresponding DS peaks. Afterpeak correction, the peak flux density of a feature is ob-tained from the flux density at the location of its peakin the background-subtracted DS.3.5. Elimination of false detections
Since M ( t, ν ) only approximates the DS, it is possiblefor it to contain some spurious peaks which do notcorrespond to real features in the DS. Only a peakin M ( t, ν ) having a corresponding peak in the DS isregarded to be a real feature. In order to weed out false peaks, the root mean square flux density ( σ ) isestimated across quiet patches in the DS as a functionof frequency. A Signal-to-Noise Ratio (SNR) for everypeak is then defined as the ratio of the peak fluxdensity to the root mean square background noise atthe frequency corresponding to the location of the peak.The spectral and temporal profiles of all peaks detectedin M ( t, ν ) were visually examined using figures similarto Fig. 5 to check for a corresponding peak in the DS.We find false detections to constitute about 24% of thetotal number of peaks detected in M ( t, ν ), all of whichhave peak flux densities, S ⊙ ,F < σ . In all, about 26%of our detections lie below the 5 σ threshold. In orderto eliminate all false positives, we reject all peaks with S ⊙ ,F < σ .Figure 6 depicts the locations of the peaks of all fea-tures detected using this automated wavelet-based ap-proach. In order to estimate the efficiency of this ap-proach at picking up features reliably from the DS,8 laypersons were presented with plots of differentbackground-subtracted DS similar to Fig. 6 and re-quested to estimate the false positive and false negativerates. According to their estimates, the CWT pipelinesuccessfully picks up features from the DS with a zerofalse positive rate but with a false negative rate of about4–6%. A total of 14,177 features were detected from 67background-subtracted DS used for this work. RESULTSThe wavelet-based analysis, yielding a large numberof features, allows us to build statistically stable dis-tributions of their properties - their peak flux densitiesand morphology in the DS. The following sub-sectionspresent the distributions of various quantities of physical utomated Solar Feature Recognition and Characterization -1 S ⊙,F (SFU) -5 -4 -3 -2 -1 d N / d S ⊙ , F ( S F U − ) α = -2.23 ± 0.02 Figure 7 . Histogram of peak flux densities on a log-log scale. interest for these features.4.1.
Peak flux densities of features
Figure 7 shows the histogram of peak flux densities( S ⊙ ,F ) of features. While this histogram extends uptonearly 307 SFU at its upper end, it touches peak fluxdensities as low as 0 . / d S ⊙ ,F ∝ S ⊙ ,F α ) fit to this histogram yields apower law index α = − .
23 over the 12 −
155 SFU range.The flux range for this power law fit overlaps with thatof the power law fits to the flux density profile done inMercier & Trottet (1997). The upper end of this fluxrange approaches the lower end of the flux range for thepower law fits by Saint-Hilaire et al. (2013). The valueof α obtained here is intermediate between the corre-sponding values obtained by Saint-Hilaire et al. (2013)( α ≈ − .
66 to − .
8) and Mercier & Trottet (1997)( α ≈ − . − . − . − .
3. The value of α obtained hereis also much lower than the power law index of − . α ≈ − . α ≈ − . − .
100 120 140 160 180 200 220 ν (MHz) P e a k f l u x d e n s i t y , S ⊙ , F ( S F U ) l o g ( N ) Figure 8 . Two-dimensional histogram depicting distribu-tions of peak flux densities and peak frequencies. The coloraxis is in log units. ∼ −
160 SFU, where we expect the distributionto be complete.Figure 8 shows a two-dimensional histogram of thedistributions of the peak flux densities and the peakfrequencies( ν ) of the features. While the peak flux den-sities of a majority of features appear to be indepen-dent of ν , a sub-population of them seem to show afrequency-dependent variation in the peak flux density.For this sub-population, the peak flux density appearsto increase with ν from 100 MHz to 150 MHz, remainnearly constant with ν between 150 MHz and 200 MHz,and then decline for ν ≥
200 MHz.4.2.
Spectral and temporal widths
Figure 9 depicts histograms of the spectral and tempo-ral widths of features. Both ∆ ν and ∆ t follow smooth,unimodal distributions. The ∆ ν distribution peaks atabout 4-5 MHz, well above the 40 kHz frequency reso-0 Suresh et al. ∆ν (MHz) N o . o f f e a t u r e s (a) ∆t (seconds) N o . o f f e a t u r e s (b) ∆t (seconds) ∆ ν ( M H z ) l o g ( N ) (c) Figure 9 . (a) Histogram of spectral widths, ∆ ν . (b) Histogram of temporal widths, ∆ t . (c) Two-dimensional histogram showingthe distributions of ∆ ν and ∆ t . The color axis is in log units. lution of our data. On the other hand, the peak in the∆ t distribution lies at 1-2 s, quite close to the 0.5 s tem-poral resolution of these data. Fig. 9(c) further showsthat the distributions of ∆ ν and ∆ t arrange themselvesin a single well-formed cluster peaking at about 4-5 MHzand 1-2 seconds. While the bandwidths of these featuresare two orders of magnitude smaller than that for type-III bursts, they are comparable to the typical frequencyspan, ∆ ν .
10 MHz (Mercier & Trottet 1997), reportedfor type-I bursts.The left panel of Fig. 10 shows a two-dimensionalhistogram of the distributions of ∆ ν and the peak fre-quency ( ν ). The prominent peak and valley-like struc-tures are artifacts arising from the limited bandwidthof observations. While valleys occur at the edges of the observing bandwidth, peaks occur at its centre. Thereseems to be a hint of a trend for a small increase in ∆ ν with increase in ν ( ≈ .
02 MHz increase in ∆ ν per unitincrease in ν ). The original data set has an equal num-ber of observations at each of the observing bands. Thealgorithm used to determine the background continuumis designed for periods of medium or low levels of solaractivity and hence, worked effectively for most of thedata. However, it was not suitable for about 24.3% ofthe data which were characterized by high solar activ-ity (typically periods immediately following occurrencesof B and C class flares) and hence, discarded from thisanalysis. This effectively leads to different observing du-rations for different observing bands and is reflected inthe left panel of Fig. 10. In order to arrive at the truespectral distribution of features, the feature occurrence utomated Solar Feature Recognition and Characterization
100 120 140 160 180 200 220 ν (MHz) ∆ ν ( M H z ) l o g ( N )
100 120 140 160 180 200 220 240 ν (MHz) E v e n t r a t e p e r un i t b a n d w i d t h ( × − s − M H z − ) Figure 10 . Left: Two-dimensional histogram showing the distribution of peak frequencies, ν , and spectral spans, ∆ ν . The coloraxis is in log units. Right: Histogram of feature occurrence rate per unit bandwidth. Spectral symmetry parameter, χ ν N o . o f f e a t u r e s Temporal symmetry parameter, χ t N o . o f f e a t u r e s Figure 11 . Left: Histogram of spectral symmetry parameter. Right: Histogram of temporal symmetry parameter on a semi-logscale. rate per unit bandwidth is computed as a function of fre-quency. As shown in the right panel of Fig. 10, the spec-tral distribution of features appears to remain flat in thefrequency range from 140-210 MHz and declines at lowerfrequencies. Below 140 MHz, the galactic backgroundtemperature rises sharply while the intrinsic solar emis-sion becomes weaker (Oberoi et al. 2016). This leads toa drop in the SNR of our detections at these frequenciesand consequently, their being under-represented in thespectral distribution. The true spectral distribution offeature occurrences is expected to be flatter than thatshown in the right panel of Fig. 10.4.3.
Spectral and temporal profiles
An interesting finding about the nature of the spec-tral and temporal profiles of the features of interest isobtained through the histograms of χ ν and χ t (defined in section 3.3.2) shown in Fig. 11. While features largelyappear to possess symmetric frequency profiles, theirtemporal profiles display no inherent symmetry. Thepeaks at the extremes of the χ ν and χ t histogram rangearise due to the presence of features with peaks locatedclose to the edges of the DS.4.4. Background flux densities at peak frequencies
Figure 12 shows a two-dimensional histogram of thebackground flux density as a function of frequency. Thebackground continuum emission shows the expectedmonotonic increase with frequency due to the broadbandthermal radiation from the 10 K coronal plasma. RSTN( )solar flux measurements estimate the median fluxdensity on the day of our observations to be 20 SFU at245 MHz. We also note that over the course of these2
Suresh et al.
100 120 140 160 180 200 220 ν (MHz) B a c k g r o un d f l u x d e n s i t y ( S F U ) l o g ( N ) Figure 12 . Two-dimensional histogram showing distribu-tions of the background flux densities at the locations ofpeaks in the DS and the corresponding peak frequencies.The color axis is in log units. observations, the spectrally smooth background fluxdensity is observed to vary by a rather large amount. DISCUSSION5.1.
Feature energies
Having estimated the peak flux, bandwidth and du-ration of each feature detected in the DS, estimates oftheir energy can be obtained if the solid angle into whichemission is radiated is known. Assuming isotropic emis-sion for these features, W = 4 πD ∆ ν ∆ tS ⊙ ,F gives thetotal energy radiated for a feature when observed from adistance D = 1 AU. As this definition of W uses only thepeak flux density of a feature without accounting for anyreduction in flux density within its shape and assumesisotropic emission, it overestimates the actual energy ofa feature. We note that while most earlier works useconstant bandwidths and durations to estimate the en-ergy radiated, we use the spectral and temporal spanscorresponding to individual features for this purpose.The histogram of total feature energies is shown in Fig.14. The typical energies of these features lie in the rangeof 10 − ergs. These features are, hence, weakerthan both the type-III bursts ( W ≈ − ergs)investigated by Saint-Hilaire et al. (2013) and type-Ibursts ( W ≈ ergs) studied by Mercier & Trottet(1997).The best fit power law to the tail of the histogramin Fig. 13 yields a power law index of − .
98. Hudson(1991) has shown that for weak flare emissions to playa significant in coronal and chromospheric heating, thepower law distribution describing their occurrence musthave index α ≤ −
2. Within the uncertainty of the fit,the features studied here meet this criterion, and hencecan be expected to play an interesting role in coronalheating. Subramanian & Becker (2004) had estimatedthe ratio of the radiative power output from noise storm W (ergs) -22 -21 -20 -19 -18 -17 -16 d N / d W ( e r g s − ) α = -1.98 ± 0.03 Figure 13 . Histogram of total energies. continua to the total power input provided to the ac-celerated non-thermal electrons producing these burststo be about 10 − − − . Using this efficiency esti-mate, the typical energies of the non-thermal electronsproducing the features of interest lie in the range from10 − ergs. This agrees well with the estimateof 10 − ergs obtained by Subramanian & Becker(2004) for the energy transferred to the non-thermalelectron population that cause noise storm continua.This hints at a possible correlation between the proper-ties of these features with that of type-I bursts. On thebasis of observational studies, Ramesh et al. (2013) alsoreport non-thermal electron energies of about 10 ergsfor radio noise storm bursts.5.2. Comparison with type-I bursts
Our statistical analysis shows that the features of in-terest appear to ride on a broadband background con-tinuum. These findings closely agree with observationsof type-I bursts present against a continuum emission,giving rise to the speculation that these features mightbe weak type-I bursts. These results would then sup-port the theory proposed by Benz & Wentzel (1981)and Spicer et al. (1982) that describes type-I burstsas observational signatures of scattered small-scale en-ergy release events in the solar corona. The very elec-trons accelerated in such small magnetic reconnectionevents might give rise to the broadband backgroundcontinuum (Benz & Wentzel 1981). Investigations oftype-I bursts in the 160-320 MHz frequency band byDe Groot et al. (1976) suggest an average frequencydrift rate of -10 MHz/s for type-I bursts. The small-scale features detected in the MWA DS appear as verti-cal streaks with no perceptible frequency drift. However,they might possess small frequency drifts which cannotbe measured at the time resolution of the MWA data. utomated Solar Feature Recognition and Characterization
10 20 30 40 50
Background flux density (SFU) P e a k f l u x d e n s i t y , S ⊙ , F ( S F U ) l o g ( N ) Figure 14 . Two-dimensional histogram of the peak flux den-sities of features against the background flux density at theirpeak frequency. The color axis is in log units. The dashedand solid black lines represent respectively the first and thethird quartile in the distribution of peak flux densities withina background flux density bin. Figure 14 shows a two-dimensional histogram of theirpeak flux densities and the background flux density attheir peak frequency. The dashed and solid black linesin Fig. 14 respectively depict trends in the first quartileand the third quartile in the distribution of peak fluxdensities as a function of the background flux density.The 25th percentile of the distribution of peak fluxesincreases from 0.76 SFU at a background flux of 2 SFUto 10.18 SFU at 38 SFU background flux. Similarly the75th percentile increases from 1.59 SFU to 44.39 SFUover the same range of background flux densities. Thisdemonstrates a tendency for an increase in feature peakflux density with an increase in background flux density.Note that we are limited by statistics at backgroundflux estimates greater than 38 SFU.As shown in Fig. 12, the background flux density atany frequency varies by a factor of ∼
2. Such large vari-ations are seen over time scales as short as 30 minutesand are not likely to reflect changes in thermal emis-sion from the 10 K coronal plasma. The observed in-crease in the peak flux of the features with increase inbackground flux density suggests a possibility that thisenhanced background continuum could arise due to a su-perposition of a large number of features which remainunresolved at the time resolution of these data. Ob-servations of such small-scale features with finer timeresolution are required to understand them better.5.3.
Comparison with type-III bursts
Gergely & Kundu (1975) and Duncan (1981) findclose similarities between sources of type-I bursts andthat of decametric type-III bursts. Benz & Wentzel (1981) claim that electrons accelerated at magneticreconnection sites, if trapped along closed field lines,produce type-I bursts and their associated continuum.If untrapped, these electrons propagate along open fieldlines and produce type-III storm bursts. Assuming atype-III-like emission mechanism for the small-scalefeatures observed in the MWA data, we arrive at a one-to-one correspondence between their peak frequenciesand the electron densities at their heights of productionin the solar corona. We assume a 4 × Newkirk (1961)density profile in the solar corona in order to translatefrom electron densities (Li et al. 2009) to heights ( h )in the solar corona. Having computed ν start , ∆ ν and∆ t for every feature, we can also determine a heightband (∆ h ) and a propagation speed ( v = ∆ h/ ∆ t )forevery feature. Assuming emission at the local plasmafrequency, we find that the features of interest mostlypossess propagation speeds of about (0 . − . h ≈ (0.01 - 0.03) R ⊙ centered at h ≈ (0.20 - 0.50) R ⊙ above the photosphere. The typical electron speeds dropsteadily with frequency, decreasing from (0 . − . . − . ∼ . SUMMARY, CONCLUSION AND FUTUREWORKWe have carried out the first detailed statisticalcharacterization of the small-scale features observed inthe MWA solar DS. Owing to their large event rates,it is very hard or impractical to manually attempt toanalyze their properties. A robust, automated tech-nique is, hence, necessarily required for our purpose.We have developed a suitable wavelet-based approachto identify, extract and characterize these features.Individual features in the DS possess unimodal spectraland temporal profiles, and a 2D Ricker wavelet is veryeffective in locating and characterizing them. A total of14,177 features have been picked up from all DS usedin this work. Though our current implementation isadapted for the MWA data, it is quite general and can4
Suresh et al. be applied to DS from other telescopes as well.The CWT algorithm enables us to reliably detectand characterize features with peak flux densities aslow as 0 . −
155 SFU flux range.We estimate the total radiated energy of these featuresto be in the range of 10 - 10 ergs. Hence, theyare much weaker than the widely studied solar type-Iand type-III bursts. Their energy distribution is fittedwell by a power law with index − .
98. Within theuncertainty of this fit, this suggests that they couldcontribute in an energetically significant manner tocoronal heating.We find these features to be quite short-lived andnarrow-band with typical durations of 1-2 seconds andbandwidths of 4-5 MHz respectively. Interestingly,while their temporal profiles display no structuralsymmetry, their spectral profiles are largely symmetricabout the peak frequency. The distribution of theiroccurrence rate remains nearly flat in the 140 - 210 MHzfrequency range. Quite analogous to type-I bursts, theyare also found to reside on an enhanced backgroundcontinuum. We speculate that these features mightcorrespond to weak type-I bursts. Since type-I burstsand decametric type-III bursts show close associations(Gergely & Kundu 1975; Duncan 1981), it is possiblethat some of these features could be weak type-IIIbursts as well.Sensitive high-time resolution observations aimed atsearching for a frequency drift and a harmonic coun-terpart for these features would hopefully provide us with crucial information for understanding them bet-ter. Imaging studies to determine their distribution onthe solar surface and investigate any correlations withother solar features will further help explore their con-tributions to coronal heating. We also hope that thisdetailed and statistically robust characterization of non-thermal emission features will engender interest in thetheory and simulation community to understand thembetter.This scientific work makes use of the MurchisonRadio-astronomy Observatory, operated by CSIRO. Weacknowledge the Wajarri Yamatji people as the tradi-tional owners of the Observatory site. Support for theoperation of the MWA is provided by the AustralianGovernment (NCRIS), under a contract to Curtin Uni-versity administered by Astronomy Australia Limited.We acknowledge the Pawsey Supercomputing Centrewhich is supported by the Western Australian and Aus-tralian Governments. AS would like to thank Atul Mo-han (graduate student at the National Centre for RadioAstrophysics - Tata Institute of Fundamental Research(NCRA-TIFR), Pune, India) for thought-provoking dis-cussions and providing timely, constructive criticisms.AS would also like to thank NCRA-TIFR for the com-putation and infrastructure support, and the KishoreVaigyanik Protsahan Yojana scheme under the Depart-ment of Science and Technology, Government of Indiafor the financial support provided during the period ofthis work.
Facilities:
MWA
Software:
Python ( ),NumPy (van der Walt et al. 2011), Scipy(van der Walt et al. 2011), Matplotlib (Hunter 2007),Scikit-Learn (Pedregosa et al. 2011) and Scikit-Image(van der Walt et al. 2014).REFERENCES
Antoine, J. P., Murenzi, R., Vandergheynst, P., & Ali, S. T.2004, Two-dimensional wavelets and their relatives(Cambridge University Press)Benz, A. O., & Wentzel, D. G. 1981, A&A, 94, 100BBowman, J. D., Cairns, I. H., Kaplan, D. L., et al. 2013, PASA,30, 31BBurnham, K. P., & Anderson, D. R. 2002, Model selection andmultimodel inference: a practical information-theoreticapproach (Springer)De Groot, T., Loonen, J., & Slottje, C. 1976, SoPh, 48, 321DDuncan, R. A. 1981, PASA, 4, 230DGergely, T. E., & Kundu, M. R. 1975, SoPh, 41, 163GGuzm´an, A. E., May, J., Alvarez, H., & Maeda, K. 2011, A&A,525A, 138GHaslam, C. G. T., Salter, C. J., Stoffel, H., & Wilson, W. E.1982, A&AS, 47, 1HHudson, H. S. 1991, SoPh, 133, 357H Hunter, J. D. 2007, Computing in Science & Engineering, 9, 90Iwai, K., Miyoshi, Y., Masuda, S., et al. 2014, ApJ, 789, 4ILi, B., Cairns, I. H., & Robinson, P. A. 2009, JGRA, 114, 2104LLonsdale, C. J., Cappallo, R. J., Morales, M. F., et al. 2009,IEEE Proceedings, 97, 1497Mercier, C., & Trottet, G. 1997, ApJ, 474L, 65MMorosan, D. E., Gallagher, P. T., & Zucca, P. et al. 2015, A&A,580A, 65MMugundhan, V., Ramesh, R., Barve, I., et al. 2016, ApJ, 831, 154Newkirk, G. J. 1961, ApJ, 133, 983NOberoi, D., Matthew, L. D., Cairns, I. H., et al. 2011, ApJ, 728L,27OOberoi, D., Sharma, R., & Rogers, A. E. E. 2016, Accepted bySoPhParker, E. N. 1988, ApJ, 330, 474PPedregosa, F., Varoquaux, G., Gramfort, A., et al. 2011, Journalof Machine Learning Research, 12, 2825 utomated Solar Feature Recognition and Characterization15