Plumelets: Dynamic Filamentary Structures in Solar Coronal Plumes

Solar coronal plumes long seemed to possess a simple geometry supporting spatially coherent, stable outflow without significant fine structure. Recent high-resolution observations have challenged this picture by revealing numerous transient, small-scale, collimated outflows ("jetlets") at the base of plumes. The dynamic filamentary structure of solar plumes above these outflows, and its relationship with the overall plume structure, have remained largely unexplored. We analyzed the statistics of continuously observed fine structure inside a single representative bright plume within a mid-latitude coronal hole during 2016 July 2-3. By applying advanced edge-enhancement and spatiotemporal analysis techniques to extended series of high-resolution images from the Solar Dynamics Observatory's Atmospheric Imaging Assembly, we determined that the plume was composed of numerous time-evolving filamentary substructures, referred to as"plumelets"in this paper, that accounted for most of the plume emission. The number of simultaneously identifiable plumelets was positively correlated with plume brightness, peaked in the fully formed plume, and remained saturated thereafter. The plumelets had transverse widths of 10 Mm and intermittently supported upwardly propagating periodic disturbances with phase speeds of 190-260 km/s and longitudinal wavelengths of 55-65 Mm. The characteristic frequency (3.5 mHz) is commensurate with that of solar p-modes. Oscillations in neighboring plumelets are uncorrelated, indicating that the waves could be driven by p-mode flows at spatial scales smaller than the plumelet separation. Multiple independent sources of outflow within a single coronal plume should impart significant fine structure to the solar wind that may be detectable by Parker Solar Probe and Solar Orbiter.


INTRODUCTION
Solar plumes have been observed since humans first gazed in wonder at solar eclipses, yet their origin and dynamics are not yet thoroughly understood (see reviews by Wilhelm et al. 2011;Poletto 2015). These common features are bright collimated structures within darker, larger coronal holes; after formation, they persist from hours to days. Plumes are both denser and cooler than their surroundings, and emanate from nearly unipolar magneticflux concentrations at the photosphere. At sufficiently high spatial resolution and temporal cadence, plumes are not monolithic, stationary structures. Instead, they exhibit both transverse fine structure and dynamic features that travel along bright threads. The nature of the propagating disturbances remains controversial, as it is difficult to determine whether they represent bulk motions, waves, or some combination of both.
Early evidence for filamentary structures within plumes at 10 ′′ in size was reported from analyses of data from Skylab (Karovska et al. 1994) and the Solar and Heliospheric Observatory (SOHO; DeForest et al. 1997;DeForest & Gurman 1998). Ofman et al. (1997) and DeForest & Gurman (1998) measured oscillations in plumes with periods of several minutes and propagation speeds above 100 km s −1 , which they interpreted as evidence for compressive sound or slow magnetoacoustic waves. Ofman et al. (1999) replicated key features of the observations with analytical and numerical calculations of magnetohydrodynamic waves in a model plume. Based on these and other observations, Wang (1998) proposed that plumes result from interchange magnetic reconnection between the ambient unipolar magnetic field of coronal holes and bipolar flux emerging from below the photosphere.
Subsequent observations from Hinode and the Solar -Terrestrial Relations Observatory (STEREO) revealed abundant small-scale structures and dynamics at the base of plumes. Raouafi et al. (2008) found strong associations between polar coronal-hole jets and/or coronal bright points and the formation/enhancement of plumes. Gabriel et al. (2009) concluded that diffuse, large-scale plumes are aggregates of numerous smaller-scale plumes. Applying wavelet analysis to coronal images of plumes from the Solar Dynamics Observatory (SDO) enabled Raouafi & Stenborg (2014) to extract even finer-scale structures and dynamics, which they referred to as "jetlets" and plume transient bright points. These features are reduced-scale versions of the jets and bright points that had been suggested previously as sources of plume plasma. Observations show that coronal jets and plumes could be connected in at least two different ways. X-ray jets precede some plumes (Raouafi et al. 2008), while small-scale jets (i.e., jetlets) may sustain long-term plume evolution (Raouafi & Stenborg 2014). This dynamic activity was associated with mixed magnetic polarity, although the large-scale magnetic flux at plume footpoints is dominantly unipolar. On the other hand, more recent studies have found that most coronal-hole jets do not result in plumes (e.g., Kumar et al. 2019), so the relationship remains unclear.
Other recent work has reinforced the association between the plumes and small-scale footpoint activity, including investigations of a bright-point/plume association over a 40-hour interval (Pucci et al. 2014), combined SDO and Interface Region Imaging Spectrograph (IRIS) observations of tiny jets and bright points (Pant et al. 2015), and strong correlations between plume emission and converging flows in the magnetic network (Wang et al. 2016). In the last study, small-scale closed loops were observed at the base of plumes even when the minority-polarity magnetic flux was weak, diffuse, or undetectable, suggesting that mixed magnetic polarity is a crucial ingredient in the activation and maintenance of plumes.
To gain further insights intothe mass source and physical mechanism generating plumes, we have applied stateof-the-art image-processing techniques (DeForest et al. 2016;DeForest 2017) to a long-duration, high-resolution, high-cadence sequence of SDO observations of a plume in a low-latitude coronal hole on 2016 July 2-3. Our detailed, quantitative analysis of this unprecedented data set has yielded key new information on the plume filamentary structures, which we refer to as "plumelets," and their quasiperiodic oscillations.
We find that our coronal plume contains well-developed transverse structure, which is masked by diffuse background emission and becomes evident after proper image processing. The number of detectable plumelet edges increases systematically with the average plume brightness, reaching more than 30 once the plume is fully formed. The plumelets exhibit a persistent pattern of upward propagating quasiperiodic disturbances, which could be a manifestation of repetitive plasma outflows, slow mode waves, or both. The frequency of the disturbances is consistent with that of solar p-mode oscillations, while the phases of the disturbances observed in different plumelets are effectively random. Our analysis indicates that plumelets could be formed by bursty reconnection in a topologically complex lower corona, driven by a combination of random and periodic photospheric motions.
Plumelets, as defined in this paper, represent a fundamental and persistent attribute of the coronal-hole plume that we observed, and may well be common to all plumes. For instance, a close look at the images of three plumes analyzed by Raouafi & Stenborg (2014) reveals ubiquitous plumelet-like structures similar those studied here. Although such fine structure has been known to exist in plumes for a long time, we are not aware of any previous studies focusing on their quantitative properties such as their characteristic spatial and temporal scales, or on their relationship with the long-term evolution of the entire plume in which they are embedded. Our present paper addresses these important aspects of coronal plumelets based on a rigorous statistical analysis of an extended set of high-resolution images of a representative plume system observed over several consecutive days. We are currently analyzing the impulsive jets and bright points at all scales at the base of this plume to determine their connections to the plumelets (Kumar et al. 2020, in preparation).
We present the data set and its processing in §2, describe our results in §3, discuss their implications in §4, and provide our conclusions in §5.
2. DATA Figure 1 reproduces two Level-1 (L1) images that show the overall corona, coronal hole, and location of the plume that we examined, which was observed from 2016 July 2 00:00 UT through 2016 July 3 17:00 UT. A large data gap followed this observing period, and by July 4 the plume was too diffuse to yield useful results. We used the full-time-resolution image sequence from SDO/AIA (Lemen et al. 2011;Lemen et al. 2012) in the 171 Å, 193 Å, and 211 Å channels, post-processed to enhance sensitivity. Plumes are bright in these wavelengths, whose collisionalexcitation temperature responses peak in the 0.8-2.0 MK range that dominates in most cool coronal structures (DeForest et al. 1991;O'Dwyer et al. 2010;Lemen et al. 2011). We also examined the 335 Å channel (Figure 2, lower-right panel) but the count rates were negligible ( 7 counts per frame at the plume base), so we excluded those images from further analysis. The data were post-processed in four distinct steps: The region of interest containing the selected plume. We analyzed the plume evolution from 2016 July 2 00:00 UT through 2016 July 3 17:00 UT; this image is near the center of our study, in Interval 5 (below). "haze" by removing predetermined/calibrated stray light from the images, and requires the full frame before cropping. Figure 2 shows the resulting images of the plume near the beginning of the data set on 2016 July 2.  Figure 1, at 2016 July 2 19:35UT in three AIA channels as marked, revealing the relative brightness and detail visible in each channel. These images show L1 data that have been PSF-corrected but otherwise unmodified.    Figures 2-4, at the same time, reveals the ionization temperature structure at a glance. Cooler (bluish) colors are dominated by the 171 Å passband at or below ∼ 1MK, while warmer (reddish) colors are dominated by the 211 Å passband at or above ∼ 2 MK. Both direct transfer (left) and minsmooth unsharp masking with noise gating (right) show multiple temperatures, but only the latter reveals the fine-scale structure in the plume. See V. M. Uritsky et al. (2020) for the animated version of the processed image set.
2. We processed the data to create the images using the noise-gating technique described by DeForest (2017).
Noise gating separates additive noise from an image sequence without imposing the loss of resolution inherent in smooth kernel convolution. It treats the image sequence as a 3D data set, and works by identifying features that are coherent across space and/or time. Features are discriminated from constant or image-dependent additive noise, based on the known statistical and variational properties of the noise. Because AIA noise levels are dominated by the Poisson statistics of photon counting ("shot noise"), we used a gating threshold scaled with the image brightness B in each 12×12×12 pixel subregion, to account for the B 1/2 dependence of shot noise. We set the gating factor of the algorithm to γ = 3. This process reduced the photon shot noise by more than a factor of 10 in each AIA channel, while preserving faint structures in the image sequence.
3. We unsharp-masked the noise-gated data using the minsmooth algorithm (DeForest et al. 2016). Minsmooth is an image operator that produces a smooth background model from an image by finding the minimum (or low percentile) pixel value in the neighborhood of each pixel in the image; the neighborhood size is set by an aperture radius parameter. Unlike convolutional smoothing, minsmooth produces an estimated minimum background model based on feature scale. Subtracting this model yields an approximately positive-definite image that contains only features smaller than the aperture used for the operator. We applied minsmooth with an aperture radius of 20 AIA pixels. The minsmooth unsharp masking highlights the spatially sharp and/or variable components of the plume. Figures 3 and 4 show the importance of the noise gating step by contrasting simple unsharp masking with minsmooth (our step 3) with noise gating followed by unsharp masking (our steps 2 and 3), respectively. The noise gating reduced shot noise by a factor of 10 or more, without reducing spatial or temporal resolution. This enables analysis of the temporal variations in the narrow, transient structures in this plume cluster.
4. We combined the 211 Å, 193 Å, and 171 Å channels into a single color movie sequence ( Figure 5) for further inspection and analysis, by placing each channel into the corresponding (red, green, blue) channel of an RGB triplet and encoding the triplet for display (Celik et al. 2018). Figure 5 is a typical RGB composite image from approximately 25,000 such composites during the 3.5-day interval under study. The relative line strengths have been adjusted by their mean strengths over the entire disk, to produce "relative scaled" color. Bluish features emit mainly in the 171 Å line (0.8 MK). Green and red (or orange) features emit primarily in the hotter 193 Å and 211 Å lines, respectively. The color scheme immediately reproduces the familiar result that closed and quasi-closed structures (such as the quiet corona and streamer at right) are hotter than open coronal structures (such as the plume and coronal hole at left).
The left panel of Figure 5 shows direct, cleaned, extreme ultraviolet (EUV) radiances; the right panel uses minsmooth unsharp masking to highlight small structures. A broad plume (bluish) is visible near the center of the region of interest, extending to the left. The plume itself is readily seen to comprise both a broad, long-lived structure and narrow, more transient substructures (denoted "plumelets" in this paper).An overall "coronal haze" pervades the region, even in the surrounding coronal hole.
With any EUV images visible haze could be due to in-instrument scattering or real coronal effects (DeForest et al. 2009). These images were corrected using the validated, instrument-specific, inverse PSFs developed by Poduval et al. (2013), and hence have instrument scattering reduced by a factor of ∼30 compared to AIA Level 1 data. We conclude that the observed haze in our processed frames is a real effect, not due primarily to instrumental scatter.
We analyzed these images to determine the origin, structure, and fluctuation characteristics of the plume, as described in §3.
The digital movie associated with Figure 4 reveals many small motions that are not apparent in the raw data. We analyzed these motions using 2D time-series analysis (e.g., DeForest & Gurman 1998) as described in 3.2. The procedure relies on a time-distance representation of quickly propagating image features, in which it is essential that the sampling interval ∆t between the analyzed image frames be constant. Figure 6 shows that the complete data set was not fully continuous, with gaps as long as ∆t = 408 s. We selected several intervals of at least 500 images that were missing at most one image at 12-s cadence during the interval. For our spatiotemporal analysis, we identified nine such intervals in the analyzed image set (labeled 0 through 8 on Figure  6 and in the following discussion). Since data gaps in the selected intervals accounted for less than 0.2% of the interval duration, no time interpolation was necessary to fill the gaps.   Table 1 provides information on the nine selected time intervals: the start and end date and time, the duration, and the number of AIA 171 Å images collected. To characterize the small-scale filamentary structure in the plume, we enhanced the bright rays with the Interactive Data Language (IDL) version of the widely used Roberts cross edge-detection operator (Roberts 1963): where m and n are the discrete coordinates of image pixels and R is the Roberts transform of the original image L. The first and the second terms on the right hand side are obtained from the image convolution using the masks and approximate the absolute value of luminosity gradient in the directions parallel and perpendicular to the principal diagonal of the L matrix, respectively. The edge-enhancement performance of the IDL routine is known to be the same as that of the original algorithm based on the root-mean-square value of the two gradient terms (Roberts 1963). We chose the Roberts operator for its robustness in the presence of high noise levels (Bonny & Henno 2018), which is essential for solar image processing. The resulting Roberts transform applied to AIA images is an uncalibrated proxy for the magnitude of the image-plane gradient of the squared column plasma density: where n e is the volumetric electron number density and the integral is performed along the line of sight. Figure 7 shows Roberts transforms of AIA 171 Å images taken in the middle of eight contiguous time intervals (labeled with index n) identified in Figure 6. Interval 0 is very similar to Interval 1 and is not shown. The grayscale coding represents the R range between 0 and 30 units; for comparison, the maximum value of L is ≈ 1000 units. A highly complex, filamentary structure of the interior region of the plume is evident. The thickness of the smallest R features, presumably representing boundaries of the individual plume rays, are close to the image resolution; the largest features reach more than 10 pixels across.
Another notable feature of the Roberts-transformed images is long-term systematic evolution of the plume structure over the entire period covered by the image set. Intervals 1 to 3, covering the first ≈ 15 hours from the start on 2016 July 2, are characterized by a sparse and clustered ray pattern, with some of the plume sectors containing multiple features and others showing no such features. Interval 3 reveals the smallest number of identifiable ray edges. After that, an increasingly refined and uniform ray structure gradually develops, reaching its mature state by Interval 8 during which rays permeate nearly the entire plume. Later, the fine-scale structure becomes less pronounced, which could in part be explained by changes in the plume orientation relative to the line of sight.
We constructed 1D profiles of the R maps along a set of arc-shaped virtual slits, labeled from a to e in Figure 8, with the center of curvature chosen to ensure that most of the slits cross the plumelet edges at an approximately 90 • angle.  Figure 8 shows five such profiles for time Intervals 3 and 8 exhibiting, respectively, the least and most developed ray structure. To reduce the noise level, the profiles were averaged over ±5 pixels across the slits, normalized by their maximum values, and shifted along the vertical axis for easier comparison. Signals from the slits labeled a and b, shown in gray, are less reliable because they were contaminated by the underlying coronal moss seen through the optically thin plume material.
Compared to the earlier Interval 3, Interval 8 shows much more structure in the 1D profiles, which reveals a multitude of small-scale rays over the entire length of the slits (Figure 8). The transverse scale of the plumelets (i.e., the distance between adjacent peaks in R) during Interval 8 ranges from ∼ 1 ′′ to ∼ 10 ′′ ; depending on the detection approach, one can identify 30-35 individual features in slits c-e, thus yielding a characteristic transverse scale of 5-6 ′′ .
This scale estimate was verified by a Fourier analysis of an extended set of arc-shaped profiles averaged over multiple time steps. The region between slits c and e in Figure 8 was resampled by 11 equally spaced slits that were used to construct 11 1D profiles, smoothed over ±3 pixels across the slit. For each profile, we computed a spatial power spectrum using a Hanning-filtered Fast Fourier Transform. The procedure was repeated for every fifth image collected during the specified time interval, after which the power spectra obtained for the individual image slits and time steps were averaged:P Here, P m,n is the local power spectrum computed for the mth slit (m = 1, ..., 11) of the nth image, k is the spatial frequency, andP is the averaged power spectrum. Figure 9 shows two pairs of the slit-and time-averaged spatial power spectraP for the original (L(x, y)) and the Roberts-transformed (R(x, y)) images during Intervals 3 and 8. The spectra computed using the original images show no identifiable small-scale structure across the plume, consistent with the very small dynamic range of the rays in the L(x, y) images compared to the background signal. In contrast, the spectra of the edge-enhanced images have welldefined characteristic scales that are different for the two intervals. These scales are manifested in the form of a spectral "break" separating frequency ranges described by small and large log-log slopes. This type of spectral signature is indicative of a stochastic spatial pattern, as in a random telegraph signal for example, as opposed to isolated spectral peaks produced by periodic patterns. The spectral break characterizing Roberts-transformed plume images in Interval 8 is located near the spatial frequency ≈ 0.2 arcsec −1 , again yielding a transverse scale of about 5 ′′ . This spectral estimate is consistent with the shape of the profiles in Figure 8.
Scale measurements based on the edge-enhanced R images should be interpreted with care. Because R(x, y) is an unsigned gradient-based measure, the characteristic sizes of the structures in the original L(x, y) images should be  about two times larger that those in the transform. Thus, for instance, a periodic structure L ∝ sin(k ⊥ s) in the direction . Fourier power spectra of AIA 171 Å intensity fluctuations across the plume for Intervals 3 and 8 before and after applying the Roberts transform. The spectra were averaged over multiple time steps, and over 11 arc-shaped slits positioned between the slits labeled "c" and "e" in Figure 8 to avoid moss contamination near the plume base. The spectra of the original images (thin lines) are shifted up by two decades for easier comparison with the spectra of Roberts-transformed images (thick lines) used to estimate the characteristic spatial scale. s perpendicular to the plume's main axis would be transformed into the structure R ∝ [1 + cos(2k ⊥ s)] 1/2 , described by the wavelength π/k ⊥ that is half that of the structure in the original image.
Therefore, the actual ray scale is about twice the spatial separation between the ray edges derived from the Robertstransformed images (Figure 8) or the characteristic spatial scale derived from the spectral analysis ( Figure 9). The estimated transverse Roberts scale of 5-6 ′′ therefore corresponds to a characteristic ray scale of 10-12 ′′ .
The average intensity of the 171 Å emission from the entire plume region (as defined by the outer boundary of the adjustable curvilinear grid described in the next section) versus the number of detectable plumelets on the edgeenhanced Roberts-transformed images is plotted in Figure 10. As mentioned above, the filamentary structure of the plume evolves slowly through two phases. During Intervals 0 and 1, a moderate number of edges produce a plume of nearly peak brightness. Intervals 1 through 6 manifest a decrease followed by an increase in both the number of ray boundaries identified in the Roberts-transformed images and the average plume brightness. Finally, the plume reaches a saturation stage during Intervals 6 to 8 when an increase in the number of plumelets is no longer accompanied by an average brightness increase. Because we used stray-light-deconvolved AIA data, the observed dependence of the midscale diffuse brightness on the number of small-scale bright features is unlikely to be attributable to PSF effects in the EUV data. Regardless of the temporal sequence, Figure 10 clearly shows that the two plume properties are strongly correlated throughout the early evolution of the plume structure. While concurrent evolution of two empirical parameters does not necessarily imply a casual connection, the correlation between the amount of plume filamentation and brightness demonstrated in Figure 10 suggests that well-developed fine structure may be required to form an optically intense plume.
Several alternative edge-enhancement methods, including the Sobel-Feldman operator and a simple shift-difference filter, were also tested (the results are not shown). The plumelet structures identified with these additional methods were found to be nearly indistinguishable from those obtained using the Roberts transform, confirming that the detected Figure 10. Average brightness of the plume (relative units) during the first nine time intervals (Intervals 0 through 8) shown in Figure 6, versus the approximate number of plumelet edges identifiable in the Roberts-transformed images of the plume. Vertical error bars represent standard errors of average brightness measured in different longitudinal slices (see Figure 11); the horizontal error bar shows a typical uncertainty in counting the edges.
fine-scale structure of the plume is independent of the particular edge-enhancement algorithm used and represents an objective physical phenomenon.

Plume evolution
To characterize the longitudinal and transverse dynamics of the plume rays in each interval from Figure 6 systematically, we constructed by hand a fan-like feature-matched coordinate system, following Uritsky et al. (2013). The coordinate system is overlaid on the plume during Intervals 1, 5, and 7 in Figure 11. The feature-matched coordinates x ′ and y ′ run across and along the visually identified "grain" of the plume in the image plane. We divided the x ′ and y ′ axes into 10 windows numbered 1 to 10, as shown in Figure 11, each covering an identical range of polar angles. Figure 11. Examples of adjustable curvilinear grids used to construct time-distance plots in longitudinal and transverse directions relative to the plume. The background snapshots are AIA 171 Å images in the middle of each time interval, preprocessed as described in §2. The dark blue (bright red) colors of the chosen rainbow palette correspond to 0 (1000) AIA DN units.
To visualize the propagating disturbances and structures present in each of the slices defined by the feature-matched grid, we averaged the pixel brightness over y ′ for each transverse slice, which resulted in 10 transverse time-distance plots. We also averaged pixel brightness values over x ′ for the 10 longitudinal slices, to obtain time-distance plots reflecting the dynamics along the plume. Figure 12 shows an example of time-distance plots in the longitudinal and transverse directions relative to the main plume axis during Interval 7. To extract significant features above the background emission, the plots were detrended in the spatial and temporal directions by subtracting respectively 5th and 3rd order polynomials fitted to the data, following a previously tested methodology (Uritsky et al. 2013). The residual fluctuations comprise between 9 and 25% of the total emission in the slits.  Figure 11. The plots were detrended in order to enhance transient features.

Transverse structure
In the transverse time-distance plots (left panels in Figure 12), the narrow elongated features aligned with the horizontal axis are the evolving signatures of the rays discussed in §3.1. The number of features per plume crosssection varies significantly with slit location and the observation time. Slits with indices i = 3 − 7 show the most pronounced filamentary structure, with a typical number of identifiable rays ranging between 9 and 13. Considering the average size of these slits (≈ 120 ′′ ), the characteristic scale of these features should lie between 9 ′′ and 13 ′′ . This agrees well with the estimated widths based on our analysis of 1D profiles from the Roberts-transformed images ( §3.1). The slits with i = 1, 2, located near the plume base, contain fewer detectable rays (only 7 to 8), possibly due to a stronger projection effect resulting in multiple overlapping rays. The slit with i = 10, farthest from the plume base, has the weakest detrended signal, making it difficult to identify any features.
The positions of the rays on the transverse time-distance plots are relatively stable; however, their brightness seems to vary slowly over a time scale of about 5000-6000 s (80-100 min). This period cannot be reliably measured because it is comparable to the duration of the entire observation. In addition, it may be affected by the polynomial detrending applied to the plots. A close inspection of the time-distance plots also reveals a quasi-periodic, checkboard-like pattern in the maxima and minima of the detrended brightness, indicating that the emission crests and troughs at neighboring locations across the plume are out of phase.

Longitudinal structure
The longitudinal time-distance plots for Interval 7 (right panels in Figure 12) reveal an intense, consistent pattern of quasi-periodic, upwardly propagating disturbances. In §4 we discuss whether these disturbances are bulk motions, waves, or a combination. The characteristic time lapse T between the features ranges from T ≈ 200 s in the longitudinal slits with indices j = 2, 3 to T ≈ 750 s in slits with j = 6, 7. The slits indexed j = 8, 9 exhibit the most stable and intense observed over the entire slit length, with period T ≈ 200-400 s as estimated by a visual inspection of the time-distance plots. We derive more precise estimates of T from a quantitative analysis of these plots in §3.2.3.
The slits showing consistent dynamic activity across a wide range of projected altitudes manifest a clear tendency for the phase speed of the propagating front to increase with distance from the plume base. An example of such dispersion is shown in the top panel of Figure 13, which provides a zoomed-in view of the longitudinal time-distance plot in slice j = 9 during Interval 7. The front highlighted with two straight line segments, indicating its approximate local slopes, is significantly steeper at high altitudes than at low altitudes; its speed apparently increases from about 140 km s −1 near the plume base to about 330 km s −1 in the upper half of the slit. A similar acceleration versus height can be seen in other fronts. The observed speed increase likely represents a partial rotation of the velocity vector caused by the curved shape of the magnetic flux tubes near the base of the plume, as illustrated in Figure 13 (bottom). By approximating the curved portion of the plume with a circular arc, the apparent speed change can be related to the angular velocity and thus to the geometry of the plume: where v 1 and v 2 are the image plane projections of the phase speed observed respectively before and after the front entered the curved path, t 1 and t 2 are the corresponding times, r c is the radius of curvature of the lower portion of the plume, and v is the true speed of the front. We solved Equation 5 numerically to estimate v. By substituting v 1 = 140 km s −1 , v 2 = 330 km s −1 , and t 2 −t 1 ≈ 250 s for the front featured in Figure 13 (top) and r c = 160 Mm for the average curvature radius as suggested by shape of the plumelet edges (see Figure 7), one obtains v ≈ 400 km s −1 . This speed is greater than v 2 because of projection effects. The corresponding plume angle characterizing its approximate orientation with respect to the image is about 34 • . This angle is substantially smaller than cos −1 ( x 2 +ȳ 2 /R S ) ≈ 62 • representing the local normal to the Sun's surface, wherex andȳ are the average plume coordinates in the image plane and R S is the solar radius. Therefore the plume is not radially oriented but is tilted away from the local normal by about 30 • .

Longitudinal fluctuations
To extract longitudinal dynamics from the data, we produced detrended charts of the evolution along specific locations in the plume. To ensure that we captured spatial variations, we detrended only in time: each pixel location in the flat-fielded, tracked movie was treated as a single time series. We modeled each time series with a cubic polynomial and subtracted that polynomial from the data. This yielded absolute-brightness fluctuations in intensity units of corrected AIA DN. The resulting plots, with a common color table and dynamic range centered on 0, are displayed in Figure 14 for four time intervals. A simplified vertical cut through the plume base, where the magnetic field expands rapidly before becoming more collimated, illustrating the propagation geometry assumed in Eq. 5. The solid black curve represents a section of a flux tube defining the low plumelet geometry. The red arrows illustrate the expected change of the velocity direction due to the plume base curvature, corresponding to the change in slope shown above.
The columns in Figure 14 reflect different stages of long-term plume evolution, from the weak, sparse filamentary structure in Interval 3 to the well-developed structure evident in Interval 8. Multiple time scales and multiple speeds co-exist in single panels. For example, during Interval 5 at j=5, disturbances lasting of order 1-2 ks co-exist with others lasting only about 100 s. The disturbances at the two time scales propagate at speeds that differ by about a factor of two, as illustrated in Figure 13. Both time scales are present throughout the data set, although the relative amplitude changes across intervals and across j values within each interval, with the more intense high-frequency component typically present near the base.
To further investigate the short-term time variability of the plumelets caused by the propagating disturbances, we used an adaptive numerical technique (Uritsky et al. 2009;Keiling et al. 2012;Uritsky et al. 2013) designed to identify wave signals in time-distance plots with low signal-to-noise ratios, such as the plots in Figure 14. The technique is based on analysis of the velocity-dependent "surfing average" signal S(t, u) illustrated in Figure 15. The surfing average is defined as a time-dependent mean of the time-distance plotL(y ′ ,t) computed along a set of straight parallel world lines running through the starting position y ′ = 0 at different times, with a fixed slope u representing an assumed phase speed: Here, t is the running time, y S = u (t ′ − t) is the spatiotemporal averaging path, and H is the propagated distance in the image frame.  Figure 11. The plots highlight the longitudinally propagating disturbances and capture the differences between the rising and mature phases of the plume development. The disturbances show the most consistent periodic pattern across all plume sectors during Interval 8, when the plume was extremely filamented, whereas the pattern is barely detectable during Interval 3, when the plume contained only a few substructures (see Figure 7 for reference). Intervals 1 and 5 exhibit propagating disturbances inside the most filamented sectors of the plume. IfL(y ′ ,t) contains a periodic disturbance propagating with speed v, then the dynamic range of S(t, u) maximizes at u = v since in that case some of the surfing lines will be sitting on the troughs and the crests of the wave or flow. For a sufficiently strong and consistent signal, maximizing the value of S allows for an automatic evaluation of v (Uritsky et al. 2013). We were unable to invoke this approach here because the activity in many of the studied time intervals was irregular and weak. Instead, we detected the fronts in the detrended time-distance plots by eye. For each time interval, we selected the longitudinal slice in which the propagating disturbance was the most pronounced in terms of its amplitude, temporal stability, and spatial extent. The identified fronts were used to calculate the average characteristic value of the phase speed v based on the maximization of the dynamic range of the surfing signal (Eq. 6). The period T of the propagating disturbances was estimated from the average time delay within adjacent pairs of successfully detected wave fronts as well as from the position of the main peak on the slit -averaged Fourier spectrum.
Examples of signals obtained using the surfing average technique are shown in Figure 16. The propagating disturbances are substantially stronger in Interval 8 than in Interval 3, which has less spatial filamentation (see §3.1). The amplitude in Interval 8 varies significantly across plume sectors, and is lowest in the middle of the plume ( j = 5, 6) where the plumelet edges are the least intense (see Figure 8 for comparison).
Surprisingly, the low-amplitude oscillations in Interval 3 and the stronger oscillations in Interval 8 have similar shapes and the same well-defined characteristic frequency ≈ 3-4 mHz. The Fourier power spectra of the two signals are shown in Figure 17. These spectra were obtained by first Fourier-transforming signals from the individual slits displayed in Figure 16 and then averaging the power spectra over the slits for a better signal-to-noise ratio. For both intervals, the main spectral peaks are well above the noise floor and have comparable widths. We note that the frequency at the spectral peak is very close to the typical frequency of the p-mode oscillations that are ubiquitous in the solar photosphere (see, e.g., Uritsky & Davila 2012).
We applied the surfing average technique to all nine time intervals of plume observations. Table 2 shows the derived parameters of the dynamic activity. In addition to the apparent velocity v, we report the corrected velocity v corr = v/ cos (34 • ), taking into account the average plume angle with respect to the image plane estimated in §3.2.1. The predicted values of the apparent wavelength λ = vT and the corrected wavelength λ corr = v corr T , corresponding to the measured period T , are also provided.
The anomalously long periods in Intervals 1 and 2 reflect an unstable disturbance amplitude during these intervals. Some of the weaker fronts are missed by the method, so the average estimated T value became larger. The other intervals presented in Table 2 are characterized by shorter and mutually consistent T values ranging between ≈ 260 and 330 s, with the average period T = 288 ± 10 s. The characteristic phase speed varied between ≈ 160 and 220 km s −1 , or over the range v corr ≈ 190 to ≈ 260 km s −1 after accounting for the projection angle. The average derived wavelength in the image (plume) frame is 50 (60) Mm. The mean values shown in the bottom line of the table exclude Intervals 1 and 2; the reported uncertainties are the standard errors.
The overall average period T reported in Table 2 is in very good agreement with the peak spectral frequency for the examples shown in Figure 17, and is similarly consistent with the dominant p-mode oscillation period. . Propagating-front signals S(t, u) obtained using the surfing average technique for observed Intervals 3 (left) and 8 (right) in each of the 10 longitudinal slits. Note the significantly higher amplitude during Interval 8 compared to Interval 3. The enhancement of the propagating disturbance in the plume tends to coincide with its spatial filamentation in both space and time.

Transverse fluctuations
The concurrent dynamic activity in nearly all plume sectors provides an additional means of evaluating the transverse characteristic scale of the plume substructure, complementary to the methods based on the analysis of static plume images described in §3.1. Depending on the underlying physics, the propagating disturbances within morphologically distinct plumelets may or may not be temporally correlated. If the quasi-periodic propagating disturbance is driven by the large-scale dynamics of the photosphere, then one would expect the behavior in different plumelets to be coherent across the plume. If, on the other hand, the activity is driven by the small-scale dynamics occurring above intrusions of minority-polarity magnetic flux sprinkled across the plume base, then one would expect the behavior in different plumelets to be largely independent and not coherent across the plume.
To determine whether collective behavior is at work in the plume, we calculated linear cross-correlation coefficients C jk ( j = 1, ..., 10; k = 1, ..., 10) characterizing the degree of temporal coherency of propagating signals between all Figure 17. Slit-averaged Fourier power spectra of propagating-front signals during the two intervals shown in Figure 16.  pairs of longitudinal plume slits j and k defined by our adjustable curvilinear coordinate systems ( Figure 11) with zero time lag: Here, S j and S k are the "surfing" signals ( §3.2) in slits j and k, obtained using the speed estimates in Table 2; σ j and σ k are the standard deviations of the two signals; and the angular bracket denotes averaging over the time interval of interest. C is an unsigned measure of temporal coherence that varies between 0 (completely decorrelated behavior) and 1 (perfect correlation). By definition, C jk = 1 for j = k and the matrix is symmetric (C jk = C k j ). Figure 18 presents two-dimensional images of the cross-correlation matrices computed for Intervals 1 to 8. The value of C quickly decreases with distance from the main diagonal j = k. The two adjacent diagonal lines (suband super-diagonals) characterizing the correlation between nearest-neighbor slits ( j = k ± 1) typically show C values below 0.5, indicating weak temporal coherence. The slits separated by larger distances show effectively no crosscorrelation for most of the time intervals plotted in Figure 18. The seemingly longer correlation scale characterizing slits 7-10 during Intervals 7 and 8 is an artifact caused by the curved plumelets crossing the grid boundaries, so that the same wave packet or bulk flow appeared in more than one slit. With this exception, the cross-correlation results suggest that the characteristic length scale describing the interactions between plumelets is not larger than the transverse size of a single longitudinal slit. This size averages 20 ′′ and ranges from 8 ′′ to 30 ′′ depending upon the distance from the plume base and the specific time interval. The correlation coefficient of the adjacent slits tends to be below 0.5 and is effectively 0 for more remote slits, suggesting a characteristic transverse correlation length on the order of slit spacing. The increased correlation scale in slits 7-10 during Intervals 7 and 8 is an artifact caused by the curved shape of those plumelets that crossed the boundaries between the slits near the plume base during this time (see the last panel in Figure 11).
Our scale estimate obviously is affected by the grid size and, therefore, should be considered an upper limit on the actual physical scale controlling interactions between the plumelets. The true correlation length, as suggested by the analysis of the plume morphology ( §3.1), may be on the order of 10 ′′ , consistent with the spectral analysis ( Figure 9). Alternatively, it may be closer to the measured minimum size of the filamentary substructure, on the order of 1 ′′ . A definitive determination of the inherent scale of the dynamic behavior and its relationship to plume morphology is an ambitious undertaking that is left for future work.

DISCUSSION
We have analyzed a lengthy sequence of observations obtained with the Solar Dynamics Observatory's Atmospheric Imaging Assembly, in which a plume originating in a mid-latitude coronal hole evolved over a few days in July 2016. Nine intervals of virtually continuous, high-resolution images at 12-s cadence were identified; their durations ranged from nearly two to four hours each. We used noise-gating (DeForest 2017) and minsmoothing (DeForest et al. 2016) techniques to process the images, greatly reducing the noise and sharpening the structure to increase the effective resolution of the data set (DeForest et al. 2018). Then we applied the classic edge-detection Roberts (1963) transform to identify and track numerous filamentary structures in the resulting noise-gated, minsmoothed image sequence. Last, we used a surfing technique (Uritsky et al. 2009(Uritsky et al. , 2013Keiling et al. 2012) to extract durations, periods, and phase speeds from the longitudinal fluctuations of the identified plumelets. Altogether, these new measurements yield key insights into the short-and long-term behaviors of solar plumes.
The most important result of our investigation is that the plume, which appears to be monolithic and more or less uniform at low spatial and temporal resolution, is instead highly structured and dynamic. While earlier studies, including those cited in the Introduction, have reported both substructure and dynamic behavior (waves or flows) in plumes, to our knowledge this prior work did not quantitatively assess the longitudinal and transverse structure, oscillation periods and amplitudes, and flow speeds as functions of time in all sectors of the plume, throughout nearly two days of observations. The descriptions to date have largely treated plumes as stationary features subject to perturbations. By more fully analyzing cleaned images of the plume with high spatial and temporal resolution, our investigation reveals that the dynamical portion of the plume is dominant over the stationary structure. This conclusion turns the conventional view of plumes upside down.
We found that the plume brightness was correlated with the number of edges detected (as many as about 30) within the filamentary structure of the plume (see Figure 10). Over the analysis period, the plumelet count decreased for 15 hours, then increased for 15 hours, while the brightness varied in direct proportion. Thereafter, the count continued to increase for another 10 hours while the brightness remained essentially unchanged. During the very late evolution, the filamentary structure fragmented into more numerous, generally narrower structures (see Intervals 6 through 8 in Figure 7).
In addition to counting the number of contained edges as a quantitative measure of substructural complexity, we measured their characteristic spatial scale (see Figure 9). The average plumelet width is about 10 ′′ (≈ 7 Mm). This is commensurate with the average size of minority-polarity intrusions of magnetic flux, EUV and X-ray bright points, and the width of EUV and X-ray jets within coronal holes. These features range in size from 2-3 times our nominal plumelet width down to scales that are much smaller, and generally they are quite numerous. Previous studies have suggested that small-scale coronal-hole jets (e.g., McIntosh et al. 2010) and even smaller-scale, more-transient jetlets (e.g., Raouafi et al. 2008;Raouafi & Stenborg 2014) are major sources of plume mass and kinetic energy. The agreement in the characteristic spatial scales of these features and our plumelets is consistent with the hypothesized connections between plume formation and persistence and the much more dynamic jets and jetlets. To test this relationship in further detail, we are performing a separate analysis of the small-scale, transient brightenings at the base of the observed plume and their associations with the plumelets reported here (Kumar et al. 2020, in preparation).
By aligning grids along and across the overall plume orientation (Figure 11), detrending the longitudinal data in time to remove the background slow evolution (Figure 14), and performing a surfing-average analysis of the residual signals (Figure 15), we determined characteristic values of about 1500 s and 200 km s −1 for the duration and speed, respectively, of the disturbances (see Table 2). The former is typical for the duration of larger-scale coronal-hole jets; the latter is typical of the dense outflows in such jets. This agreement with jet lifetimes and speeds strengthens the case for direct links between jets/jetlets and plumes. At the same time, however, the measured speeds of the longitudinal disturbances are fully consistent with their being slow magnetosonic waves. In early analyses of observations from the Solar and Heliospheric Observatory's Extreme-Ultraviolet Imaging Telescope, Ofman et al. (1997) and DeForest & Gurman (1998) identified similar long-lived, slow-moving disturbances in coronal holes and plumes as such waves. The surfing-average analysis reported here provides additional compelling evidence that timeperiodic wave activity occurs on the plumelets (see Figure 16).
The dichotomy of jet-like versus wave-like behavior of the observed plumelets can be resolved straightforwardly if the plumelets are generated by jets and jetlets at the plume base. Each plumelet extends along the plume at a rate determined by the speed of the outflowing jet plasma, and the lifetime of the plumelet structure is determined by the duration of the underlying jet from the source. Clearly, the dense jet is a compressible flow, which would be expected to generate slow magnetosonic waves. The resulting dense plumelet also could support waves generated by external forcing mechanisms.
Our surfing-average analysis further determined the characteristic oscillation periods of the longitudinal fluctuations in our plumelets. We found an average period across the entire data set of 288 ± 10 s (see Table 2). By taking power spectra of the signals along individual plumelets to extract their frequencies, we found a prominent peak near the corresponding frequency of 3.33 mHz (see Figure 17). These results strongly suggest a connection with the pressuredriven p-modes that flex the solar photosphere, although in principle this agreement could be coincidental. For the 10-Mm characteristic width of the plumelets, the corresponding latitudinal and longitudinal standing-mode numbers are ℓ, m ≈ 200. Substantial power is present in the p-mode spectrum at these and shorter wavelengths (e.g., Duvall et al. 1997). Hence, the photospheric undulations driven by the p-modes can be much smaller than the plume base under study, and comparable in size to the plumelets and to the typical source regions of coronal-hole jets. The undulation periods (≈ 300 s) are significantly shorter than the typical total duration of the plumelets and the jets (≈ 1500 s). This suggests a scenario in which each plumelet source region rides up and down multiple times on the p-mode waves as the jet/jetlet is generated and the plumelet is formed. Consequently, the dense plasma that becomes a distinct plumelet oscillates at the frequency of the underlying photospheric oscillations, and this modulation appears in the plumelet wave signals.
The cross-correlation analysis of different longitudinal strips demonstrates that the plumelet waves were not correlated at transverse scales of 20 ′′ (15 Mm) and above (see Figure 18). This is consistent with the p-modes providing the photospheric modulation, whose spatial variations occur on scales as small as the typical separation between the plumelets. Although the underlying frequency of the p-modes is common to all plumelets, there is no fixed relationship between the phases of the observed oscillations in any one plumelet and in its neighbor. Hence, the plumelet wave motions are uncorrelated.
High-cadence, high-resolution observations of coronal-hole jets, especially those that exhibit helical outflows along the jet spire, provide strong evidence that nonlinear Alfvén waves accompany jet launch and lead the flow of dense 20 plasma into the outer corona (e.g., Kumar et al. 2019), consistent with our jet simulations (e.g., Karpen et al. 2017;Uritsky et al. 2017;Roberts et al. 2018). Our analysis did not reveal any Alfvén waves traveling along the plumelets. This is not a strong null result, however, nor is it very surprising, for three reasons. First, as mentioned already, the outflow of dense plasma trails the Alfvén waves at the jet front, and increasingly so at higher altitudes due to their speed differential. The dense plasma flow should be expected to provide the main contribution to emission from the elongating plumelet. Second, although an Alfvénic shock wave may lead the entire procession from the jet source region, the density enhancement at the shock is not necessarily large, and the Alfvén wave itself is incompressible even though it is nonlinear. The contribution of the Alfvén wave to the plume emission, therefore, may be insignificant. Third, Alfvénic undulations of the plumelets could be difficult to detect in the relatively low-cadence (12-s) SDO/AIA observations. For a typical coronal Alfvén speed of 1 Mm s −1 and a source-region size of 10 Mm or less, the expected Alfvén-wave periods are 10 s and below, so these high-frequency waves cannot be properly resolved by SDO/AIA.
Our results strongly indicate that the filamentary structure is imprinted dynamically on the plume at its base, and that this structure and dynamics persist out to at least 100 Mm above the surface. It seems unlikely that this coherence can be explained by any mechanism other than its association with the structure of the plume magnetic field, which guides the plasma outflow and the plasma waves. So guided, the plumelet structure and dynamics that we have observed may persist well out into the heliosphere, perhaps to where they can be measured by Parker Solar Probe (PSP; Fox et al. 2016) and/or Solar Orbiter (SolO; Müller et al. 2013). The plumelets' plasma density is enhanced relative to that in the surrounding open corona, and this enhancement should be sustained in the absence of a mechanism that selectively accelerates and rarefies the gas at the front of the plumelet. The characteristic transverse scale (∼ 10 Mm) of the plumelets in the corona should expand geometrically, due both to the diverging spherical geometry of space (∝ R/R S ) and by the square root of the super-radial factor (≈ 2-4) by which the area of open magnetic flux on the solar surface expands to fill all 4π sr of space sufficiently far from the Sun. The anticipated scale thus ranges from roughly 200 Mm at 10 R S up to 600 Mm at 30 R S . It is an enticing prospect that PSP and/or SolO might directly detect the 5-minute oscillations imprinted on the plumelets. If our identification of the solar p-modes as the driving mechanism for these waves is correct, this would mean a direct detection of the global dynamics of the solar interior at least 10 R S out into the heliosphere. We look forward to seeing the exciting new data that come back from these missions and comparing them to our expectations, which are founded on the detailed analysis of a coronal plume that was observed so diligently by SDO.

CONCLUSIONS
Filamentary structures and motions in plume images have been reported for many years (e.g., Raouafi & Stenborg (2014) and refs therein). Here we presented the first in-depth quantitative investigation of these structures, which we denote "plumelets". Using a large set of high-resolution, high-cadence solar coronal images, we quantified the highly dynamic nature of the plumelets, and demonstrated that their impulsive behavior may, in fact, dominate the largescale behavior of the "host" plume.
We processed almost 40 hours of nearly continuous observations of a typical solar coronal plume by SDO/AIA on 2016 July 2-3. Image-processing, edge-detection, and signal-analysis techniques enabled us to obtain the following results: 1. The plume comprised numerous (10 to 15) filamentary substructures, which we refer to as "plumelets", that accounted for most of the variable plume brightness over its lifetime.
2. The width, length, and duration of the plumelets averaged approximately 10 Mm, 100 Mm, and 1500 s, respectively.
4. The longitudinal fluctuations were uncorrelated at transverse scales of 20 Mm and above, i.e., from one plumelet to another.
5. The plumelet oscillation period agrees very closely with the peak-power period of the solar p-modes (∼3-5 min).
6. The plumelet width, duration, and speed are consistent with those of the dense outflows in typical coronal-hole jets.

21
We are currently investigating the suggestive connections between the plumelets and jetlet activity observed by SDO at the plume base.