© The Authors 2023

1. Introduction

The solar atmosphere is replete with distinct types of waves caused by perturbations of the plasma (Stein & Leibacher 1974). In particular, the magnetic field, ubiquitously present as diverse structures in the solar atmosphere, acts as an effective restoring force for perturbations, giving rise to different types of magnetohydrodynamic (MHD) waves (see reviews in Nakariakov et al. 2005; Van Doorsselaere et al. 2009; de Moortel 2009; Jess et al. 2015, 2023; Jess & Verth 2016). The magnetic field also provides an optimal wave conduit, as compressible and incompressible MHD-wave modes are readily generated in magnetic flux structures by recurring perturbations around the solar surface (Verth & Jess 2016). These MHD waves have long been suggested to act as a key player in the transport of energy and the resulting heating of the upper solar atmosphere (e.g., Osterbrock 1961; Porter et al. 1994; Klimchuk 2006; Srivastava et al. 2021). The focus of the current study lies in assessing the energy budget of waves and expanding on tentative existing candidates for dissipation of wave energy into local plasma heating (i.e., Ionson 1978; van Ballegooijen et al. 2011; Grant et al. 2018).

Compressible MHD wave modes are determined by their ability to perturb the local plasma density. An example of these are the compressible MHD sausage mode, whose observational signatures are periodic fluctuations in the cross-sectional area of their cylindrical waveguide, accompanied by periodic fluctuations in intensity. The mode of the oscillation can be discerned from the phase relationship between the cross-sectional area and total intensity oscillations, with slow sausage modes displaying an in-phase behaviour, as well as an anti-phase relationship for fast sausage modes (Edwin & Roberts 1983; Moreels et al. 2013). With the advent of high-resolution instruments such as the Rapid Oscillations in the Solar Atmosphere imager (ROSA; Jess et al. 2010), in-depth studies of sausage modes in the lower solar atmosphere could be conducted (e.g., Morton et al. 2011). Sausage modes have since been positively identified in a variety of magnetic structures, such as sunspots (Dorotovič et al. 2014), magnetic pores (Grant et al. 2015; Keys et al. 2018; Gilchrist-Millar et al. 2021), and chromospheric fibrils (Morton et al. 2012; Gafeira et al. 2017). More recently, Guevara Gómez et al. (2021) utilised the unique observations of the Atacama Large Millimeter/submillimeter Array (ALMA; Wootten & Thompson 2009) to reveal the signatures of fast sausage modes in small-scale features previously undetectable by intensity and velocity observations.

Incompressible MHD (transverse) kink modes manifest as motions perpendicular to the vertical waveguide axis. If the cross-sectional area of the waveguide is perpendicular to the line of sight, these modes are observed as horizontal motions with an oscillatory pattern. Various structures in the photosphere and chromosphere, such as small-scale bright points, fibrils, mottles, and spicules, have been identified as waveguides whereby kink modes are excited, having periods on the order of 30–350 s and velocity amplitudes of about 1–29 km s⁻¹ (e.g., De Pontieu et al. 2007b; Lin et al. 2007; Kuridze et al. 2012; Stangalini et al. 2013, 2014, 2015, 2017; Jess et al. 2015; Jafarzadeh et al. 2017a,b). Guevara Gómez et al. (2021) also observed transverse oscillations in three small-scale bright features observed in solar ALMA observations with periods on the order of 30–136 s and velocity amplitudes between 0.6 km s⁻¹ and 5.9 km s⁻¹, although the structures associated with the waveguides were not identified.

Evaluating atmospheric plasma parameters, such as temperature, is essential in calculating the energy flux of MHD waves. The inherent complexities in spectral line formation in the chromosphere result in simplified estimations of these parameters, such as interpolated models (Grant et al. 2015) or photospheric inversions (Gilchrist-Millar et al. 2021). However, the introduction of solar ALMA observations provides a novel insight into the chromosphere, through observations of millimetre continuum radiation formed mostly by thermal bremsstrahlung, also known as free-free emission (see, e.g., Valle Silva et al. 2021; Wedemeyer et al. 2016; White et al. 2006; Dulk 1985, and references therein). ALMA records the brightness temperature, which is analogous to the actual gas temperature of the continuum forming layer, providing unprecedented accuracy in chromospheric plasma temperature measurements. ALMA observations can also be tailored for observing oscillations and waves in the solar atmosphere, with high temporal cadences of up to 1 s (see, e.g., Jafarzadeh et al. 2021b; Patsourakos et al. 2020; Narang et al. 2022; Nindos et al. 2021; Eklund et al. 2020).

For this study, we utilise high cadence observations of brightness temperatures with ALMA for the detection and analysis of wave modes in the solar atmosphere, supported by artificial millimetre observations based on numerical simulations. The ALMA data and the numerical simulations are described in Sect. 2, followed by a description of detection method for oscillation events in Sect. 3. The analysis of the detected events and exhibited wave modes are presented in Sect. 4. A discussion and conclusions are provided in Sect. 5.

2. Data

In this study, solar ALMA observations are employed alongside synthetic continuum maps from a Bifrost simulation, both in the original spatial resolution and in a degraded resolution to match ALMA’s spatial and temporal resolution.

2.1. Observations

The ALMA Band 3 (2.8–3.3 mm) and Band 6 (1.2–1.3 mm) observations were carried out on 22nd April 2017, as detailed in Jafarzadeh et al. (2021b). Band 6 observations took place from 15:59 to 16:44 UTC followed by Band 3 observation between 17:20 and 17:55 UTC, both as part of the program 2016.1.00050.S. The reported solar P-angles, which determine the rotation of the ALMA field of view (FoV) with respect to the solar north, were −25.32° and −24.65° for bands 3 and 6, respectively. Each of the observed bands consist of four sub-bands which are centred at different wavelengths within the band range (see, e.g., Wedemeyer et al. 2020; Eklund et al. 2020). For the purpose of this study, we used the full-band image reconstruction mode in which all the sub-bands of a particular band are combined into one, resulting in time series for frequencies around 3 mm for Band 3 (100 GHz) and 1.25 mm (240 GHz) for Band 6. The full-band ALMA images have a higher signal-to-noise ratio (S/N) compared to the sub-bands S/N because of the better sampling of the spatial Fourier space during the image reconstruction. As an interferometer, ALMA measures the temperature differences on the observed source. It is therefore necessary to use the additional total power (TP) maps to obtain the absolute temperature values for the observed source. The TP maps are full-disk images of the Sun that are derived from observations with the (single-dish) total power antennas of ALMA using the fast scanning technique described by White et al. (2017). The absolute temperature ranges of the ALMA observations are [5106,15 126] K with a mean of 9317 K and a standard deviation of 1229 K for Band 3 and [2090,12 196] K with a mean of 7496 K and a standard deviation of 1014 K for Band 6.

The Band 3 observation consists of a time series split into 3 scans (blocks of observation), with a duration of about 10 min each and a cadence of 2 s. There are gaps of about 2 min between scans that ALMA uses for observing calibration sources. The pixel size chosen during the reconstruction process is 0.34 arcsec. The ALMA synthetic elliptical beam size defines the interferometric spatial resolution of the observations and depends on the observed frequency, time of day, and atmospheric conditions. For one specific frequency, the elliptical beam shape smoothly changes as time progress. In this observation lasting for about 35 min including calibration breaks, the median area of the synthetic elliptical beam was 2.94 arcsec² with a standard deviation of 0.04 arcsec² throughout the whole observation, and its dimensions were about 2.2 arcsec × 1.7 arcsec (1600 km × 1235 km assuming that 1 arcsec ≈ 727 km), with a clockwise median inclination of 28.7° with respect to the solar north.

The Band 6 observation consists of a time series split into five scans, each with a cadence of 2 s. The first four scans are ∼8 min long, with the final scan lasting ∼4 min. There are calibration gaps of about ∼2 min between scans. The pixel size chosen during the image reconstruction process was 0.14 arcsec. In this observation lasting about 45 min (including breaks), the median area of the ALMA synthetic elliptical beam was 0.45 arcsec² with a standard deviation of 0.03 arcsec². The median shape of the beam was about 0.84 arcsec × 0.67 arcsec (610 km × 488 km), with a clockwise median inclination of 86.1° with respect to the solar north.

These ALMA observations have been spatially co-aligned with observations from the Solar Dynamic Observatory (SDO; Pesnell et al. 2012), consisting of photospheric magnetograms from the Helioseismic and Magnetic Imager (HMI; Scherrer et al. 2012). The solar coordinates of the centres of the field of views at the start of observations are (x, y) = (−246, 267), (arcsec) with a radius of 33 arcsec and (x, y) = (−255, 263), (arcsec) with a radius of 16 arcsec for bands 3 and 6, respectively. Both bands sample a plage region with strong magnetic fields spanning a range of [−1084,59] Gauss in the photosphere. These data products are presented in Jafarzadeh et al. (2021b) as data sets D2 and D8, where the authors showed that the magnetic topology at chromospheric heights represents vertical or near-vertical fields in a large portion of the FoVs with more inclined field towards the edges. This magnetic field topology verifies that these data sets are good candidates to search for MHD waves. In particular, the Band 3 observations have already proved to exhibit small-scale blob-like bright features with oscillations in temperature and size which reflect the characteristics of fast-sausage MHD-modes, with a transverse component likely being kink MHD-modes (Guevara Gómez et al. 2021).

The left panels of Fig. 1 shows the contextual TP or full-disk maps for the Band 3 (top) and 6 (bottom) observations, characterised by the bright signatures of an active region. The area of interest for this study is the plage region denoted by black dashed squares, with small black circles showing the FoV of the ALMA interferometric observation. The middle-left panels show the SDO/HMI images corresponding to the black dashed square in the TP maps at the start of the observation, with red circles denoting the Band 3 FoV. This SDO/HMI image reveals the magnetic nature of the underlying photosphere. The middle-right panels display a composite of the Atmospheric Imaging Assembly (AIA; Lemen et al. 2012) at 170 nm (pink) and 17.1 nm (yellow), and shows the structuring of the coronal magnetic field. The right panels show Band 3 and 6 snapshots with the beam shape visualised in the bottom-left corner by an ellipse. The white circle in the top panel illustrates the ALMA Band 6 FoV, with the green dots are the median location of all the small-scale bright features analysed in this study (see Sect. 4). Finally, the white square with the blue cross in the Band 6 panel marks the position of a single feature which is used for illustrative purposes in this paper (see Fig. 3). Through inspection of the images, it is clear that Band 6 can resolve magnetic structures with greater fidelity in comparison with Band 3.

Fig. 1. Total power or full-disk map in ALMA Band 3 (top row from left to right). Small black dashed square delimits the region showed in SDO/HMI and SDO/AIA 1700-171 and the small black circle depicts the ALMA Band 3 interferometric FoV, SDO/HMI image showing the presence of magnetic field in the ALMA FoV indicated by the red circle, SDO/AIA 170-17.1 nm composite image of the same region, ALMA Band 3 image in which the green dots indicate the median location of all analysed feature and the white circle shows the ALMA Band 6 FoV. Bottom row: same as in the top row but for Band 6, in the ALMA Band 6 image, the small white square with the blue cross in the centre signalise the location of the example feature used in this paper. The small black ellipse in the bottom-left corner of ALMA images show their respective beam shapes.

The specific time series used in this study can be downloaded from the Solar ALMA Science Archive (SALSA)¹ as data sets D02 (Band 3) and D08 (Band 6; for further details on SALSA see Henriques et al. 2022). The Solar ALMA Library of Auxiliary Tools (SALAT)² can be used for convenient data cube reading and loading of basic properties in header (for more details on SALAT, see Jafarzadeh et al. 2021a).

2.2. Simulation

The simulation used in the current study was created with the 3D radiation magnetohydrodynamic (RMHD) code Bifrost (Gudiksen et al. 2011) and represents a magnetic enhanced network region surrounded by the quiet Sun (Carlsson et al. 2016). The simulation includes non-local thermodynamic equilibrium (non-LTE) radiative transfer and an equation of state that handles non-equilibrium hydrogen ionisation (Leenaarts et al. 2015). The computational box extends from the upper convection zone to the corona with an extent of 24  Mm × 24  Mm × 16.8  Mm in the (x, y, z) directions respectively. In the vertical direction (z), the simulation uses an adaptive grid of 496 cells whose sizes vary between 19 km and 98 km, being of about 20 km at chromospheric heights, 14.4 Mm out of the 16.8 extent above the photosphere. For the horizontal directions (x, y), the simulation has a grid of 504 cells with a constant size of 48 km, so that the pixel size of the simulation is equivalent to 0.066 arcsec.

The simulation sequence used here has a cadence of 1 s and lasts for about 1 h of solar time, reaching an state of quasi-equilibrium after 200 s (Eklund et al. 2021a). The physical set-up of the simulation evolves in time with periodic boundary conditions in (x, y) directions and open boundary conditions in the vertical direction. The magnetic field strength has an average unsigned value of 50 G at photospheric heights. The initial condition is comprised of two photospheric regions of opposed polarity with a separation of 8 Mm between them, which connect by magnetic loops in the atmosphere above (for a detailed description of the simulation see, Carlsson et al. 2016). This simulation is of interest to this study as MHD oscillations have been detected and analysed in fibrillar structures at coronal heights (Kohutova & Popovas 2021). This implies that MHD-modes may be present throughout the atmosphere, including the height range sampled by ALMA.

We used the advanced radiative transfer code (ART; de la Cruz Rodríguez et al. 2021), which solves the radiative transfer equation under the assumption of local thermodynamics equilibrium (LTE). In this work, ART uses the non-equilibrium electron densities provided by the Bifrost simulations, as this is important for the resulting continuum radiation at sub-mm and mm wavelengths (for an extended discussion of the use of simulations in the context of ALMA, see, Wedemeyer et al. 2022). The final data product consists of two time series of continuum brightness temperatures at the frequencies of ALMA Bands 3 and 6 (see Sect. 2.1). In order to ensure a fair comparison between the observations and the simulation, we had to match the time resolution to those of ALMA. For this, we have arbitrary chosen the snapshot corresponding to 700 s as starting point³ and so on removed every other frame, constructing a new time series wit the frames 700, 702, 704 and etc., until eventually we matched both ALMA observations time stamps. In the case of Band 3, the 3 scans are composed of the snapshots [700–1288, 1428–2022, 2172–2764]. For the case of Band 6, the 5 scans correspond to the snapshots [700–1182, 1294–1770, 1913–2389, 2499–2981, 3120–3354].

After matching the time resolution of ALMA, we have also matched the spatial resolution so we can compare results between the original and the degraded (ALMA equivalent) resolution for each band. The spatial resolution was amended through convolving each snapshot with the corresponding ALMA synthetic elliptical median beam to simulate ALMA’s angular resolution. Please note that this approach assumes a perfect sampling of the spatial Fourier space while, in reality, ALMA’s Fourier space coverage is sparse. Simulating these effects in detail, however, would go beyond the scope of the study presented here. Subsequently, the pixel size of the simulation was changed from 0.066 arcsec to the respective 0.34 arcsec in Band 3 and 0.14 arcsec in Band 6. The resulting time series of degraded mm brightness temperature maps altogether with the equivalent original-resolution time series are the ones used for the analysis on this paper (see Sect. 4). Figure 2 shows the same snapshot of the Bifrost simulation for the synthetic ALMA Bands 3 and 6 in the top and bottom row, respectively. The images in the left column correspond to the original simulation, which roughly spans a range of temperatures between 2700 K and 20 000 K for Band 3 and 2750 K and 16 000 K for Band 6. The absolute temperature ranges of the degraded simulations are [3954,14 280] K for Band 3 and [3005,12 726] K for Band 6. The mean temperature vales are the same in both cases, with 7337 K for Band 3 and 5779 K for Band 6. In all the simulations, it is possible to visually identify small structures and the typical fibrillar nature of the quiescent chromosphere. The images in the right column correspond to the simulation after convolving with the ALMA beams and resizing to the pixel size of the ALMA data, namely, the so-called degraded simulation. The small white ellipses in the bottom-left corners represent the interferometric beams. The most evident effect of the convolution is the degradation of the spatial resolution, which blurs the structures in the simulation, giving them a more blob-like shape. The convolution also has an impact on the temperatures, reducing the range of observed values, which has to be taken into account when searching for temperature enhancements. The green dots in all the panels mark the median location of all the small-scale bright features detected in the simulations and the blue cross in the Band 6 degraded snapshot mark the location of a feature which has been used to illustrate the analysis (see Fig. 8). In addition, it is worthwhile noting that the area of the FoV in the Band 3 observation is roughly four times the area of the FoV in the Band 3 simulation. A complete description of the impact of the different ALMA observing modes in this Bifrost simulation can be read in Eklund et al. (2021b).

Fig. 2. Snapshots at the same time from Bifrost simulation for the synthetic ALMA Band 3 and 6 are shown in the top and bottom rows, respectively. Left column corresponds to original snapshots with full resolution. Right column shows the result of convolving and matching with ALMA Beam and pixel size. The green dots mark the median position of all analysed features. The blue cross by (20,20) arcsec in the Band 6 degraded simulation shows the position of the example feature.

3. Feature identification and tracking

A visual inspection of the ALMA observations reveals the presence of transient events in the shape of small-scale bright features that appear to move while their brightness and size vary (see movies in SALSA). Moreover, Guevara Gómez et al. (2021) showed by means of wavelet transform analysis that three of these small-scale features found in the Band 3 observation were likely of magnetic nature, indicating that their oscillatory behaviour are signatures of MHD-modes. However, these isolated cases were not sufficient to make a statistically significant conclusions. Accordingly, in this study we aimed to analyse the oscillatory behaviour of as many features as possible, in both the two observed ALMA bands, (Sect. 2.1) and the two synthetic continuum map time series from the Bifrost simulation (Sect. 2.2). Therefore, it was necessary to implement an algorithm to automatically identify and track the features across the time series such that the output from this algorithm could be used to perform a statistical analysis.

Identifying and tracking varying features in a series of ALMA images is not a minor task, as the chromosphere contains a variety of features, each with different sizes, shapes, and intensities. These structures can commonly interact and overlap within the FoV, resulting in extraordinary shifts and morphological changes from one frame to the next one. To constrain the problem, it was imperative to establish some threshold characteristics that any of the bright features had to meet in order to be selected by the algorithm. In particular, the feature size was the main property taken into account. By manual inspection, we found that the sizes of the features were around the sizes of the corresponding ALMA elliptical beams (i.e., the spatial-resolution proxy). This fact was in agreement with the previous results of Guevara Gómez et al. (2021), with sizes between 0.95 and 1.27 times the median beam’s major axis.

Based on a size threshold, and following the essence of the method of Crocker & Grier (1996) to identify, track, and extract quantitative information from colloidal suspensions in noisy images (implemented to small-scale solar magnetic elements by Jafarzadeh et al. 2013), we applied to each frame a difference of Gaussians (DoG) filter with Gaussian widths of half and double the median major axis of the ALMA beam. This filter enhances the blob-like features making their detection easier (see Appendix A). Then, for each frame, we saved the position of all local maxima separated by at least 1.5 times the median major beam. Subsequently, we used the positions of the local maxima as seed points in the Flood fill method implemented in the Python Scikit-image library (van der Walt et al. 2014). The flooding method is applied to a normalised, inverted colour scale of the ALMA images, namely, those where bright pixels were made dark and vice-versa. The flooding algorithm begins at a starting point or seed position (each local maximum) and then checks if the pixel values of surrounding positions fall within the so-called tolerance range with respect to the value of the seed point. This range was chosen such that for each seed point the method would flood all surrounding pixels within ∼2.4 standard deviations, corresponding to the full width at half maximum (FWHM) understood as all the pixels where the temperature is halfway between the pixel seed position (local maximum) temperature and the median temperature of the complete FoV. This yields a mask of flooded pixels for each frame, in which it is possible to identify the isolated features as small flooded regions completely surrounded by non-flooded pixels. As we are interested in blob-like features of sizes similar to the beam size, all the features with distances between any two points on their borders greater than four times the median beam major axis were discarded. This step serves as a filter for large features, but also for features with elongated shapes. This process yields a cube of masked images in which every individual mask represent a possible feature of interest. To track the features in time, the overlapping of the masks is compared frame by frame, meaning that if the mask seen at the time t_i − 1 shares at least 2 pixels with the mask seen at the time t_i then they are considered as belonging to the same feature. In the case that two or more masks merge in the following frame, each individual feature is considered alive until the previous frame. When all the features are individualised, each of those obtains and obtains an unique identifier. This process results in a cube that contains the masks of all individual features and their unique labels. Such a cube is produced for both the observational and simulated ALMA data for both bands and has the same dimensions as the respective original data cube. Finally, the maximum brightness temperature of each unique feature in each frame is extracted and used to construct light-curves for each feature. These light curves are classified by means of the Density-Based Spatial Clustering of Applications with Noise (DBSCAN) algorithm which is part of the Python Machine Learning Scikit-learn package (Pedregosa et al. 2011). This classification results in different clusters of light-curves, namely, oscillating light-curves (see e.g., Fig. 8), light-curves with a sole transient brightening (see e.g. Fig. 1 at Eklund et al. 2020) and light-curves with any other different behaviour. We have then chosen all the features belonging to the cluster of light-curves with oscillatory behaviour to be analysed in this study.

4. Analysis and results

4.1. Occurrence and general properties

The method described in Sect. 3 was applied to the time series of ALMA observations and simulations in both bands. In total, 193 features were identified in Band 3 and 293 in the Band 6 observations. In the simulation at the degraded resolution, 24 and 204 features were identified in Bands 3 and 6, respectively. While in the original resolution 36 and 392 features in Bands 3 and 6 were identified, respectively. The raw number of features does not tell anything about how similar the results in the observation are compared to the simulation. Thus, to be able to compare the results, it is necessary to take into account the duration of each time series and the area of the FoVs. The resulting number density values, which are defined as the total number of features per square mega-metre per second, are listed in Table 1. These density values fall within the same order of magnitude for the observation and simulation in each band. However, the densities in the observational data are approximately double for Band 6 and triple for Band 3 than the densities found in the degraded simulation for those bands. While the density values found at the simulation with original resolution are closer to those in the observation, specially for Band 6. This is a direct consequence of the different spatial resolutions at each band. An additional possible explanation for at least part of this systematic difference may be due to the differences in magnetic topology of the observed and simulated field. The FoV of the observations is mostly covered by the magnetic enhanced region with long fibrillar structures (see Fig. 1), whereas in the simulation the magnetic enhanced region is comparatively smaller and confined to the central FoV region only (see Fig. 2). Corresponding differences in the occurrence of transient bright features as analysed here are then expected. We also note that the events are a factor of 5–6 more frequent in Band 6 than in Band 3, which is seen in both the observations and simulations. The higher angular resolution of Band 6 compared to Band 3 (a factor of 2.3) may play a part, allied with the height differentials of the atmospheric layers mapped by Band 3 and Band 6 (see Sect. 5).

For each individual feature, the algorithm returns its shape, area, value of the brightest pixel, average pixel value within the feature, initial and final time, and the position of the centre of gravity (of intensity) as computed from the pixel values within the feature’s border as weights. For the purpose of this paper, all the areas were converted to the diameter of a circle of equal area. This equivalent diameter is hereafter referred to as the size of the feature.The representative brightness temperature of the features has been chosen to be the value of the brightest pixel with median and standard deviation values in Table 1.

4.2. Location, lifetimes, size, and horizontal velocity

Figure 3 illustrates the outcome of the method previously described. From left to right, it shows the time evolution of one of the many features identified and traced in Band 6 which has been highlighted in the bottom-right panel of Fig. 1 as a blue cross in the middle of a white square. In Fig. 3, the range of the colour map is uniform throughout the 4 panels which correspond to the times indicated as vertical lines in the top left panel of Fig. 8. The black contours correspond to the borders of the feature. Although the shape of the feature changes with time, it remains close to a blob-like shape, making the diameter of a circle of equal area an excellent proxy for its size and its corresponding evolution. The rightmost panel shows the displacement of the centre of gravity during the lifespan of the feature, while the arrow shows the direction of time. Each displacement of the centre of gravity in x and y directions between two consecutive frames was used to calculate the instantaneous horizontal velocities, v_x and v_y, and the corresponding instantaneous horizontal velocity $v_{\mathrm{t}}=\sqrt{v_x^2+v_{\mathit{y}}^2}$ for each feature.

Fig. 3. Time evolution of a feature in Band 6 observation. Four panels correspond to times marked as vertical lines in the top-left panel of Fig. 8. In each panel, the black contours depict the feature border found by thee algorithm. The right-most panel shows the trajectory of the centre of gravity with the arrow indication the direction of time.

Figure 4 shows the distribution of lifetimes for all features found in the observation in blue, features found in the degraded simulation in grey and features found in the original simulation as a solid red line. Band 3 is shown in the top panel and Band 6 in the bottom panel, in both bands, the peaks of the distributions are between 80 and 120 s. The lifetimes of the observations are a bit shifted to the right compared with the lifetimes of the simulation.

Fig. 4. Distribution of lifetimes for all the features for Band 3 (top) and Band 6 (bottom). In both panels, the blue histograms correspond to results from observations and the grey to results from degraded simulations. The solid lines are their corresponding log-norm fits. The red solid line corresponds to the log-norm fit of the results from the simulation at the original resolution.

Table 1 summarises the range, median and standard deviation of the sizes and Fig. 5 shows the distributions of sizes of all the features at all time frames in observations and simulations. Band 3 is shown in the top panel and Band 6 in the bottom panel. In both panels, the vertical maroon dashed lines are the median minor and major axes of the ALMA beam in each corresponding band, showing that the features found by the method were completely resolved by ALMA. The grey distributions correspond to results from the degraded simulations and the blue distributions to results from observations. The blue and black thick lines shows a log-norm fit in each case. The red solid lines show the log-norm fit of the results from the original resolution simulation. The peaks of the distributions are located at larger sizes in degraded simulations compared to observations although the shapes of the distributions are alike. The distributions of sizes in the case of the original simulation are very similar to the distributions found in the observations. This may be due to the fact that the degrading of the resolution results in a wider FWHM as the median temperature remains the same but the higher temperatures are diminished. However, as discussed in Sect. 5, the results between observations and simulations do not differ significantly, implying that the difference in the peaks may be due just to a discrepancy in the spatial power spectrum between observations and simulation given by the apparently more grainy pattern present in the observations. This is likely caused to some extent by ALMA’s Fourier space coverage of the observations being sparse whereas perfect Fourier space coverage is assumed for the simulations (see also Sect. 2.2).

Fig. 5. Distribution of sizes for all the features at all frames for Band 3 (top) and Band 6 (bottom). In both panels the blue histograms correspond to results from observations and the grey to results from degraded simulations. The solid lines are their corresponding log-norm fits. The red solid line corresponds to the log-norm fit of the results from the simulation at the original resolution. The vertical maroon lines shows the minor and major beam axes.

Figure 6 shows the distribution of the horizontal velocities for Band 3 in the top row and for Band 6 in the bottom row. Again, the grey distributions correspond to results from simulations, the blue distributions to results from observations and the blue and black thick lines correspond to log-norm fits. The distributions of horizontal velocities are shown in the left column, while the distributions of the velocity amplitudes are shown in the right column. The velocity amplitudes are computed by taking the fast Fourier transform of each horizontal velocity oscillation and subsequently choosing those amplitude values greater than two standard deviations over the average. Table 1 summarises the ranges, means and median velocity amplitudes in the four cases.

Fig. 6. Distributions of total horizontal velocities (left column) and amplitude velocities (right column) for degraded simulations (grey) and observations (blue) in the features found and traced in Bands 3 (top row) and 6 (bottom row). The blue and black solid lines correspond to a log-norm fit in each corresponding case. The red solid lines correspond to the log-norm fits of the results from the simulation at the original resolution.

Figure 7 shows the distributions of formation heights at τ = 1 in each band from the degraded simulation at the top panel and the original simulation at the bottom panel. The shadow regions represent the heights inside the features while the solid lines show the heights outside the features. Interestingly, the radiation is formed higher in the solar atmosphere within the features both in the degraded and original simulations. This higher formation is more noticeable in Band 6 than in Band 3. This fact is in agreement with the behaviour found by Loukitcheva et al. (2015) and Eklund et al. (2021b) in the same simulation. In particular, the radiation formed above regions of stronger magnetic field at ALMA wavelengths comes from higher heights when compared to the radiation formed above quieter regions.

Fig. 7. Distribution of geometric heights from the original (top panel) and degraded (bottom panel) simulations inside and outside the feature borders. The grey-black colours correspond to values in Band 3 and the red colours to values in Band 6.

4.3. Wavelet analysis and periods

The time evolution of the features temperature, size, and horizontal velocity exhibits an oscillatory behaviour. In order to determine the oscillation periods, we performed a continuous wavelet analysis (CWT) for the temperature, size, and velocity for each individual feature. This allows for the analysis of a time-varying signal in both frequency and time domains by determining the dominant modes of oscillation in the power spectrum. A detailed mathematical description of the method can be found in Torrence & Compo (1998). In addition, it is also possible to perform a cross-wavelet analysis to identify the phase difference between two signals. The cross-wavelet is useful when we know that the variation of two signals is linked, for example, due to physical reasons. Wavelet analysis has been widely used to characterise oscillations in the solar atmosphere (see, e.g., Bloomfield et al. 2004; Jess et al. 2007; Jafarzadeh et al. 2017a; Guevara Gómez et al. 2021, and reference there in).

To proceed with the analysis of these oscillations, all time series corresponding to temperature, size, and horizontal velocity were detrended (using a linear fit) and apodised (using a Tukey window with a length of 0.1). Figure 8 shows the temperature (black lines) and size (blue lines) oscillations of two features in Band 6 in the top row. The Band 6 detection (left column) corresponds to the same feature shown in Fig. 3, with the degraded simulation detection (right column) exhibiting a similar morphology, highlighting the complementary nature of features detected in observations and simulation. The continuous wavelet transforms for temperature, size, and cross-wavelet are shown, where black contours highlight the power regions above a 95% confidence level (for further analysis). For these two particular features, the dominant periods for both temperature and size are around 32 s. The cross-wavelet plots also show small arrows within the confidence level. These arrows represent the phase lag of the two signals, being in phase when they point downwards and anti-phase when they point upwards, visualising the phase difference between size and temperature oscillations at the dominant periods within the confidence level. The horizontal velocity oscillations and their continuous wavelet transforms look similar to those already reported in Figs. 5 and 8 of Guevara Gómez et al. (2021) where the velocities for three features in the Band 3 observation reach values up to 20 km s⁻¹ and average oscillating periods between 43 and 72 s. Consequently, the plots regarding the velocities are not shown explicitly here.

Fig. 8. Temperature and size oscillations of an example feature in Band 6 observation (left column) and degradedsimulation (right column) are shown in the top row. From the second row to bottom the continuous wavelets for temperature, size and the cross-wavelet for the same two quantities are shown. For these panels the period is plotted using a log2 scale in the y-axis and time is shown in x-axis. The hashed black region highlights the cone of influence for each spectrum and the black solid lines mark the 95% confidence levels. The arrows in the bottom panels depict the phase difference between temperature and size oscillations, being in-phase when arrows point straight up, anti-phase when arrows point straight down, pointing to the left indicate that perturbation of temperature follows those of size and vice versa when pointing to the right.

The wavelet analysis described above is applied to each individual feature, resulting in the dominant oscillation periods (i.e., all periods corresponding to wavelet areas with a confidence level of 95% or larger and outside the cone of influence) as well as the phase differences between temperature and sizes. Table 2 summarises the mean, median and standard deviation of the dominant periods from the wavelet analysis. Figure 9 shows from left to right the period distributions of temperature, size, and velocity for simulations in grey and solid red line, and observations in blue. The distributions for Band 6 are narrower than for Band 3, with smaller periods, indicating that, generally, features oscillate faster at the chromospheric heights mapped by Band 6 as compared to the heights for Band 3.

Fig. 9. Distribution of dominant periods for Band 3 in the top row and Band 6 in the bottom row, for temperature (left column), size (middle column), and horizontal velocity (right column). The degraded simulation’s distributions are shown in grey, the observation’s distributions in blue and the original simulation’s distribution is show as a red line.

4.4. Phase relations

Figure 10 shows the distribution of phase differences between temperature and size for all the features in Band 3 (top) and Band 6 (bottom). In both bands, the distributions for the observations (blue-shaded histograms) clearly exhibit a predominant anti-phase relation (−180° or 180°) between size and temperature for both bands. For the degraded simulation (grey-shaded histograms), there is also a noticeable number of events with anti-phase relations. However, in Band 3, the distribution also shows substantial peaks at phases around −130°, while in Band 6 this behaviour is not present. The red solid lines shows the distributions for the features found in the original simulation. In the case of Band 3, the anti-phase relation is more clear than in the degraded simulation. Please note that the original simulation does show a peak around −90° which is likely enhanced and extended when the convolution with the ALMA beam is applied. In addition, the distributions for Band 3 and 6 of original simulations show a small bump around 0°. This ‘in-phase’ behaviour is also seen in the observations, but their relevance compared to the anti-phase behaviour is minimal. The phase relationships between structural properties have been used to derive the mode of the MHD wave (Moreels et al. 2013) and would indicate that a higher combination of slow and fast sausage modes were present in the simulation compared to the observation. However, these models depend on density and structural uniformity seen only in the photosphere, and have been shown to no longer be valid in the chromosphere (Grant et al. 2022). Instead, they suggest that that the physics present in the simulations, although similarly to the physics of the observed plage region, cannot completely account for this particular portion of the real Sun.

Fig. 10. Distribution of phase angles between temperature and size for degraded simulation (grey), observation (blue) and original simulation (red) in ALMA Bands 3 (top) and 6 (bottom).

4.5. Energy flux estimation

The time averaged flux of energy ⟨F⟩ carried by a wave through a surface perpendicular to a magnetic flux tube can be estimated using the following equation originally provided by Moreels et al. (2015):

$$\begin{array}{c}\hfill ⟨F⟩=\frac{1}{2}f(1+ln(1/f)){\rho }_{\mathrm{e}}{⟨v_{\mathrm{r}}⟩}^2{v}_{\mathrm{ph}},\end{array}$$(1)

where f is the filling factor, ρ_e is the mass density outside the flux tube, v_r is the radial velocity of the perturbation, and v_ph is the phase speed of the wave. The filling factor is computed as the averaged ratio of the area occupied by the features and the total area of the corresponding FoV. It is worth noting that a non-unity filling factor can result in over- or under-estimation of the energy flux (Van Doorsselaere et al. 2008; Goossens et al. 2013). The densities are taken directly from the simulation as the average of the mean values within circles of three times the radii of each feature excluding the densities inside the feature, being ≈2.5 × 10⁻⁸ kg m⁻³ for Band 3 and ≈4.0 × 10⁻⁷ kg m⁻³ for Band 6. The factor of 1/2 was introduced to account for the time-averaged energy over one period. For this study, the estimated energy fluxes correspond to those carried by the sausage modes.

The radial velocity, v_r, has previously been derived to estimate the properties of MHD sausage waves observed in pores (Freij et al. 2016). It is derived from the radial change of the feature under inspection (δr) and the oscillation period of the feature size (P). For the energy flux estimate presented here, we set δr to the averaged amplitude of the radius oscillation across all features and P to the averaged mean periods of size oscillations across all features. This procedure is applied to both observations and simulations separately. The radial velocity can then be written as:

$$\begin{array}{c}\hfill v_{\mathrm{r}}=\frac{2\pi \delta r}{P},\end{array}$$(2)

the value of v_r corresponds to the maximum amplitude and the time averaged value to be used in the Eq. (1) is given by $⟨{\mathit{v}}_{\mathrm{r}}⟩={\mathit{v}}_{\mathrm{r}}/\sqrt{2}$ (Morton et al. 2012). The phase speed of the wave v_ph can be expressed as:

$$\begin{array}{c}\hfill v_{\mathrm{ph}}=\frac{\omega }{\kappa }=c_{\mathrm{s}}\sqrt{\frac{\pm A-1}{\pm A-1+(\gamma -1)(h\nu /k_{\mathrm{B}}T)}},\end{array}$$(3)

where γ is the ratio of specific heats, h is the Planck constant, ν is the frequency of observation (100 GHz for Band 3 and 240 GHz for Band 6), k_B is the Boltzmann constant, T is the temperature within the feature, and c_s is the sound speed estimated using the average temperature of the observation or simulation in each case. The dimensionless amplitude ratio A is defined as:

$$\begin{array}{c}\hfill A=\frac{\delta T/T_0}{\delta S/S_0},\end{array}$$(4)

where δT is the amplitude of the temperature perturbation, T₀ is the mean temperature, δS is the amplitude of the cross sectional area perturbation, and S₀ is the mean area. Grant et al. (2015) defined and used Eq. (3) to estimate the phase speed of upwardly propagating oscillations in a pore at specific periods observed in distinct wavelengths. In their formalism, fractional changes in emergent intensity from ROSA images are used. Given that brightness temperature measured from ALMA lies within the Rayleigh-Jeans limit, it can be applied as a proxy for emission intensity.

Table 3 summarises the variables involved in the estimation of flux energy. Each value is shown with its respective uncertainty from the error propagation analysis. The large area perturbations in comparison to previous works (Grant et al. 2015) tentatively suggest again that the dominant observed waves correspond to fast-sausage modes (Edwin & Roberts 1983). The flux of energy in the Band 3 observation is 1839 ± 464 W m⁻² while the energy flux in the Band 6 observation is 5341 ± 1348 W m⁻². For the degraded simulation, the energy fluxes are 621 ± 407 W m⁻² and 3640 ± 1557 W m⁻² for Bands 3 and 6, respectively. For the original simulation, the energy fluxes are 453 ± 289 W m⁻² and 5485 ± 1587 W m⁻² for Bands 3 and 6, respectively. In all the cases, the energy flux decreases from Band 6 to Band 3. Using the median formation heights of the degraded simulation in Table 1 of 1563 km for Band 3 and 1260 km for Band 6, the slope of the energy fluxes between the two heights in the observation is 12 ± 5 W m⁻² km⁻¹.

5. Discussion and conclusions

We have analysed 193 and 293 small bright features in ALMA observations, 36 and 392 features in corresponding synthetic continuum maps from Bifrost in ALMA Bands 3 and 6 at original resolution, and 24 for Band 3 and 204 for Band 6 at the degraded resolution. It should be noted that the spatial resolution is critical for this type of study, and that any further loss of resolution would have considerable impact on the quality of data and the resulting detection of such features. The method used for the study was designed to ensure reliable detection of small-scale events down to sizes similar to the size of the ALMA beam. However, the method systematically underestimates the true size of the detected features by ∼20% for most common sizes in Band 3 and 6. We note that it does not affect the derived oscillation periods or the measured temperatures and horizontal velocities. However, it might affect the radial velocity in Eq. (2) and therefore the estimated energy flux by a factor of $\sqrt{2}$. The reliability of the method is discussed in more detail in Appendix A. The fact that the method performed equally well for the observations and simulations tested in the study clearly demonstrates the overall reliability of this approach.

The direct observation of MHD modes in the detected features implies that they are magnetic in nature. This interpretation is supported by the magnetic-field topology of the datasets analysed here, as demonstrated in Jafarzadeh et al. (2021b), indicating that the bright features are the result of magnetic fields branching from the solar surface into the chromosphere. Additional support is provided by the fact that the simulations exhibit very similar features which are clearly connected to magnetic field structures (i.e., they are cross sections of magnetic flux tubes). The properties of the features, namely, temperature, intensity, and horizontal velocity change with time (see Fig. 3) and exhibit an oscillatory behaviour (see Fig. 8). The high temporal resolution and unprecedented spatial resolution (at these wavelengths) of the ALMA observations allow reliable characterisation of the observed oscillations. It is found that the number of features per square mega-metre per second in the observations and simulations falls in the same order of magnitude for each corresponding receiver band, respectively (see Table 1). Interestingly, there are between five to six times more features in Band 6 than in Band 3 in both observation and simulation. This may be due to several factors: (i) features in Band 3 with sizes comparable to the majority of features detected in Band 6 would not be observable because of the lower spatial resolution of Band 3; (ii) on average, the continuum radiation in Band 6 is formed below the continuum in Band 3, while magnetic features tend to expand with height, directly affecting detection at a given resolution; and (iii) as a result of different formation heights, ALMA Band 3 and 6 are mapping different layers and thus different parts of the magnetic structures that are connected to the detected features. Correspondingly, different properties of the mapped magnetic environment would affect the occurrence of small-scale bright features.

The periods of oscillations in temperature, size, and velocity are identified by means of wavelet analysis (see Sect. 4 and Fig. 8) and summarised in Table 2. The horizontal velocity amplitudes are summarised in Table 1. The determined values of velocity amplitudes and oscillation periods are comparable between those obtained from the simulations and observations in each respective band. On top of that, the calculated parameters fall within typical values of MHD modes observed in distinct features, that is, respectively, 5–29 km s⁻¹ and 37–350 s for spicules (De Pontieu et al. 2007b; He et al. 2009; Jess et al. 2012; Pereira et al. 2012; Bose et al. 2019), 1–10 km s⁻¹ and 100–250 s (Jafarzadeh et al. 2013; Stangalini et al. 2013) for bright points, 8–11 km s⁻¹ and 120–180 s for mottles (Kuridze et al. 2012, 2013), and for fibrils 1–7 km s⁻¹ and 94–315 s (Pietarila et al. 2011; Morton et al. 2012). Although it is not possible to unambiguously determine what type of the aforementioned phenomena we are observing in these ALMA datasets, the agreement among the derived oscillation periods and velocity amplitudes supports the conclusion that the detected oscillations are indeed associated to MHD modes. Furthermore, we speculate that the oscillations in horizontal velocities may be associated to MHD transverse (kink) waves.

Moreover, oscillation periods larger than 100 s (corresponding to frequencies lower than 10 mHz) seem to be absent in Band 6 but present in Band 3 (see Fig. 9) with the mean and median periods being accordingly shorter in Band 6 (see Table 2). This finding is consistent with other studies using other observing facilities at different wavelengths and different spatial resolutions. For instance, high-frequency waves (> 10 mHz) in the lower-to-mid chromosphere are observed in slender fibrillar structures (see e.g., and reference there in Gafeira et al. 2017; Jafarzadeh et al. 2017a), whereas lower frequencies (< 10 mHz) are often reported to be present in the high chromosphere, for instance, for fibrillar structures observed in the Hα and Ca II 854.2 nm lines (see e.g., De Pontieu et al. 2007a; Kuridze et al. 2012; Morton 2012). However, we should note that no direct observational evidence for the dissipation of energy associated with high-frequency waves has been found to date. While the estimated decrease of the energy flux with height of the identified sausage waves from Band 6 to Band 3 could suggest dissipation of wave energy between the respective formation height ranges of the two bands, the fact that the ALMA data for Band 3 and Band 6 used here were not obtained simultaneous prevents drawing any strong conclusion regarding the differences in the found oscillations period distributions and its implications for any related decrease in the energy flux.

The analysis of phase differences between temperature and size oscillations shows a predominant anti-phase behaviour (see Fig. 10). The phase differences together with the variation of the cross-sectional area of the features could indicate the presence of compressible MHD fast-sausage modes (Moreels et al. 2013; Morton et al. 2012; Guevara Gómez et al. 2021), while the in-phase relation may be an indication of slow-sausage modes. In the real Sun, there is a mix of superimposed modes that are challenging to disentangle when observed, but whose particular signatures have nevertheless been found in the ALMA observations. This is reflected as phases out of the anti-phase and in-phase behaviours. From the simulations, we have found that these phases are mostly coming from the regions with low magnetic-field strength, that is, towards the fainter parts in Fig. 2. Thus, future works should consider analyses of other physical quantities that the simulations provide and the observations lack (such as gas density, pressure, and magnetic fields), to further explore the properties of MHD modes. Additionally, simulations can provide further information on how spatial resolution can influence the detection of such wave modes co-existing across the same small-scale magnetic structures.

In addition, we present the first estimation of the energy flux of sausage modes in numerous structures at chromospheric heights using ALMA observations. For this, we have made use of Eq. (1) in which all quantities are derived directly from our measurements except for the mass density. Hence, we have used the densities directly extracted from the simulations as explained in Sect. 4.5. As expected, the energy fluxes in Band 6, i.e., lower in the atmosphere, are larger than the energy fluxes in Band 3, that is, higher in the atmosphere. This implies that, on average, the energy transported upwards through the chromosphere (along the small magnetic structures) is damped at these heights, thus suggesting the propagation of these waves in the solar chromosphere. The propagation of MHD waves between different layers corresponding to two frequencies within the Band 3 observations has been proven by Guevara Gómez et al. (2022) who showed that transverse waves propagate with speeds of the order of 96 km s⁻¹ and carry a flux of energy of the order of 4000 W m⁻². Although the propagation of waves cannot be measured directly, the observations presented here imply an estimated energy flux of the fast-sausage mode of 1838±429 W m⁻² and 5341±1448 W m⁻² in Bands 3 and 6, respectively. The energy fluxes from the degraded simulations were approximated as 621±407 W m⁻² in Band 3 and 3640±1557 W m⁻² in Band 6. For the simulation at the original resolutions the energy fluxes were estimated as 453±289 W m⁻² in Band 3 and 5485±1587 W m⁻² in Band 6. The energy fluxes computed from both observations and simulations in Band 6 are favourably comparable within their uncertainties. The relatively higher difference in the energies between the observations and degraded simulations in Band 3 may be due to the lower spatial resolution which consequently has an impact on the measured quantities, but this may only result in relatively larger uncertainties. In fact, the energy flux in Band 6 at the original resolution is only 3% greater than the energy flux in the Band 6 observations.

The structures presented in this paper, exhibit lower energies to wave modes detected in similar small chromospheric flux tubes. Alfvén waves exhibiting energy flux of 15 kW m⁻² (Jess et al. 2009) and kink modes with energies in excess of 30 000 W m⁻² (Bate et al. 2022) have been detected in fibrils and spicules across the chromosphere. In particular, Morton et al. (2012) estimated the wave energy flux of compressible MHD sausage waves, presented in 103 chromospheric structures observed in the Hα line-core, as 11 700±3800 W m⁻². With the Hα spectral line forming in the mid-to-upper chromosphere, close to the heights sampled by ALMA (Rutten 2017; Molnar et al. 2019, 2021; Brajša et al. 2021), this infers that the detected waves are less energetic than their Hα counterparts. However, a robust comparison of ALMA band formation heights compared to other chromospheric diagnostics is lacking (Molnar et al. 2019), so the energy differential may be influenced by formation height offsets. Moreover, it is worth noting that kink modes with relatively low energies on the order of 310 W m⁻² detected in small-scale magnetic elements (Jafarzadeh et al. 2013) as well as 15 kW m⁻² in slender fibrillar structures (Jafarzadeh et al. 2017a) have also been reported for observations in the low-to-mid chromosphere, sampled by the SUNRISE Ca II H filtergrams (Solanki et al. 2010, 2017).

In light of estimated chromospheric radiative losses of between 1000 and 10 000 W m⁻² (Withbroe & Noyes 1977), it is unlikely that these waves alone can provide sufficient energy to balance the chromosphere. However, as part of a plethora of energetic chromospheric waves modes, they may play a part. Indeed, a difference in energy between the waves when observed in ALMA Bands 3 and 6 could be a signature of damping (see Sect. 4.5). Although, with the observations not being co-temporal, and with no propagation statistics, it cannot be concluded whether this is damping in identical wave trains due to upward propagation. It is still a reasonable possibility that this wave energy is damped in the chromosphere, as significant damping has already been observed in sausage modes bridging into the chromosphere (Grant et al. 2015; Gilchrist-Millar et al. 2021). In terms of chromospheric heating, any damped energy would need to be dissipated into localised plasma heating for these waves to have an influence. The obvious candidate would be leaky modes which naturally deposit energy (Cally 1986), as fast modes can readily take this form, though the fast mode classification of waves in this study is only tentative. Apart from this, compressible waves have been proposed to dissipate their energies by multiple mechanisms, notably mode conversion (Ulmschneider et al. 1991), resonant absorption (Goossens & De Groof 2001), thermal conduction (Ofman & Terradas 2003), and turbulent mixing (van Ballegooijen et al. 2011). These mechanisms act on scales below current resolution capabilities, so remain a viable, but hypothetical, method for direct dissipation. There is incentive therefore for future detailed studies of wave propagation in ALMA observations. In particular, wave mode characterisation, alongside propagation between sub-bands would be invaluable to accurately derive wave energetics and their influence on plasma heating, given the large population of these small flux tubes across the solar atmosphere.

Further research on how to exploit ALMA’s capability for measuring temperatures, as well as the use of diagnostics at other wavelengths, is necessary to unveil the potential role of the energy carried by these waves in the processes of heating. Furthermore, the possibility of co-simultaneous observations of the same target at two different ALMA bands is also a direct desired conclusion of this study, as it would facilitate the study of the same feature at different atmospheric heights, resulting in more direct conclusions regarding the energy transport and heating in such small-scale structures.

Acknowledgments

J.C.G.G wish to thank Carlos José Díaz Baso for his valuable input on the ML clustering method. This work is supported by the SolarALMA project, which received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No. 682462), and by the Research Council of Norway through its Centres of Excellence scheme, project number 262622. HE was supported through the CHROMATIC project (2016.0019) funded by the Knut and Alice Wallenberg foundation. This paper makes use of the following ALMA data: ADS/JAO.ALMA#2016.1.00050.S. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC(Canada), MOST and ASIAA (Taiwan), and KASI (Republic of Korea), in co-operation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO and NAOJ. We are grateful to the many colleagues who contributed to developing the solar observing modes for ALMA and for support from the ALMA Regional Centres. S.D.T.G. is grateful to the UK STFC for the consolidated grant ST/T00021X/1. Collaboration between authors was facilitated by support from the International Space Science Institute, Bern within the International Teams program (Team 502 – WaLSA: Waves in the Lower Solar Atmosphere). Scientific discussions with the WaLSA team (https://www.WaLSA.team) are acknowledged.

References

Appendix A: Feature tracking method

One of the main concerns of studying small-scale structures with ALMA observations relays on the level of reliability that should be placed on the found results. In this particular case, to the sizes of the found features because these are comparable to the ALMA beam during observations. The first question that arises considers whether the features are artefacts produced by either the observation process itself or by the image reconstruction of the interferometric data. To address the former possibility it is safe to assume that the quality assessment of the data by ALMA ensures that the artefacts due to the observation process are minimal. For the latter possibility, we rely on the extensive used of the interferometric techniques over the course of its history, as well as on the latest studies carried out with solar ALMA observations which have proven ALMA capabilities to distinguish the big picture, but also small dynamics structures with excellent S/N values. The analysis done in this research has proven the physical nature of the small-scale features found with the algorithm described in Section 3. Therefore, the intention of this appendix is to asses the limitation of the algorithm with regard to the definition of the features borders. To this end, we designed a control experiment consisting of a toy model image filled with Gaussian-like fake features.

The procedure, equivalent for Bands 3 and 6, is as follows. Fake 2D Gaussian-like features with known FWHM (assumed as the size) were randomly placed in a box with similar size to the respective ALMA band under study. Each 2D Gaussian has a FWHM of 0.5. The true position of each feature within the box and its respective size (FWHM) were stored to be used later on the process. The sizes of the fake features vary randomly between 80% of the average minor beam size and five times the average major beam size. Representing features of sizes in the ranges 1220-8000 km for band 3 and 400-3100 km for Band 6. Each feature is randomly assigned to a position within the box using a Poisson disk sampling method, in which the size of the box was defined as 81.6x81.6 arcsec for band 3 and 33.6x33.6 arcsec for band 6. It was also set that features should maintain a minimum separation distance between positions of at least the average major beam size. They way of placing all the features together consisted just in summing them up, that is when two features lie too close their peaks get mixed and they might become indistinguishable. In addition, a random noise with a standard deviation of 5% the intensity at the FWHM was added to each box. The resulting image is subsequently convolved with the corresponding ALMA beam for each band. The algorithm is then applied to this toy model expecting that it will identify and delimit each feature.

Figure A.1 shows in the top panel the toy model images for Bands 3 and 6. The white circles represent the FWHM of those features that the algorithm did not find, due to the impossibility of distinguishing between two different features whose signals were spatially mixed because of their proximity. The green circles represent those features that were successfully found by the algorithm. The success rate is roughly 60% for both bands, meaning that the algorithm is capable of picking six out of ten features present in the image. The second panel from the top shows the difference of Gaussians which is a crucial step in the algorithm as it is used to define the local maxima which would be use to feed the flooding method. The third panel shows the resulting masks, being the green circles same as in top panel. It is possible to perceive that the masks, despite having irregular shapes, are still blob-like shapes. The bottom panel shows the relation between the real size of the features, that is, the FWHM that were kept when the fake features were created and the sizes found by the algorithm. As a conclusion, it is possible to say that the algorithm is underestimating the size of the features. However, this does not change the results found in this paper, conversely, reinforce the results because staying within the borders of the features ensures that the amount of noise from outside is reduced. Hence, the algorithm has been tested to be an efficient way to delimit small scale structures which in addition to bright features can also be used to find dark features.

Fig. A.1. Toy model images for bands 3 and 6 (top panel). Blue crosses represent the centre location of fake features, green circles show those features picked up by the algorithm and white circles those features that were not picked up. The small white ellipses in the bottom left corner show the ALMA beam shape. Second panel shows the difference of Gaussians in each case. Third panel shows the resulting masks and green circles as same as in top panel. Bottom panel shows the relation between real size of the feature and the size found by the algorithm.

All Tables

All Figures

Fig. 1. Total power or full-disk map in ALMA Band 3 (top row from left to right). Small black dashed square delimits the region showed in SDO/HMI and SDO/AIA 1700-171 and the small black circle depicts the ALMA Band 3 interferometric FoV, SDO/HMI image showing the presence of magnetic field in the ALMA FoV indicated by the red circle, SDO/AIA 170-17.1 nm composite image of the same region, ALMA Band 3 image in which the green dots indicate the median location of all analysed feature and the white circle shows the ALMA Band 6 FoV. Bottom row: same as in the top row but for Band 6, in the ALMA Band 6 image, the small white square with the blue cross in the centre signalise the location of the example feature used in this paper. The small black ellipse in the bottom-left corner of ALMA images show their respective beam shapes.

In the text

Fig. 2. Snapshots at the same time from Bifrost simulation for the synthetic ALMA Band 3 and 6 are shown in the top and bottom rows, respectively. Left column corresponds to original snapshots with full resolution. Right column shows the result of convolving and matching with ALMA Beam and pixel size. The green dots mark the median position of all analysed features. The blue cross by (20,20) arcsec in the Band 6 degraded simulation shows the position of the example feature.

In the text

	Fig. 3. Time evolution of a feature in Band 6 observation. Four panels correspond to times marked as vertical lines in the top-left panel of Fig. 8. In each panel, the black contours depict the feature border found by thee algorithm. The right-most panel shows the trajectory of the centre of gravity with the arrow indication the direction of time.
In the text

	Fig. 4. Distribution of lifetimes for all the features for Band 3 (top) and Band 6 (bottom). In both panels, the blue histograms correspond to results from observations and the grey to results from degraded simulations. The solid lines are their corresponding log-norm fits. The red solid line corresponds to the log-norm fit of the results from the simulation at the original resolution.
In the text

Fig. 5. Distribution of sizes for all the features at all frames for Band 3 (top) and Band 6 (bottom). In both panels the blue histograms correspond to results from observations and the grey to results from degraded simulations. The solid lines are their corresponding log-norm fits. The red solid line corresponds to the log-norm fit of the results from the simulation at the original resolution. The vertical maroon lines shows the minor and major beam axes.

In the text

Fig. 6. Distributions of total horizontal velocities (left column) and amplitude velocities (right column) for degraded simulations (grey) and observations (blue) in the features found and traced in Bands 3 (top row) and 6 (bottom row). The blue and black solid lines correspond to a log-norm fit in each corresponding case. The red solid lines correspond to the log-norm fits of the results from the simulation at the original resolution.

In the text

	Fig. 7. Distribution of geometric heights from the original (top panel) and degraded (bottom panel) simulations inside and outside the feature borders. The grey-black colours correspond to values in Band 3 and the red colours to values in Band 6.
In the text

Fig. 8. Temperature and size oscillations of an example feature in Band 6 observation (left column) and degradedsimulation (right column) are shown in the top row. From the second row to bottom the continuous wavelets for temperature, size and the cross-wavelet for the same two quantities are shown. For these panels the period is plotted using a log2 scale in the y-axis and time is shown in x-axis. The hashed black region highlights the cone of influence for each spectrum and the black solid lines mark the 95% confidence levels. The arrows in the bottom panels depict the phase difference between temperature and size oscillations, being in-phase when arrows point straight up, anti-phase when arrows point straight down, pointing to the left indicate that perturbation of temperature follows those of size and vice versa when pointing to the right.

In the text

	Fig. 9. Distribution of dominant periods for Band 3 in the top row and Band 6 in the bottom row, for temperature (left column), size (middle column), and horizontal velocity (right column). The degraded simulation’s distributions are shown in grey, the observation’s distributions in blue and the original simulation’s distribution is show as a red line.
In the text

	Fig. 10. Distribution of phase angles between temperature and size for degraded simulation (grey), observation (blue) and original simulation (red) in ALMA Bands 3 (top) and 6 (bottom).
In the text

Fig. A.1. Toy model images for bands 3 and 6 (top panel). Blue crosses represent the centre location of fake features, green circles show those features picked up by the algorithm and white circles those features that were not picked up. The small white ellipses in the bottom left corner show the ALMA beam shape. Second panel shows the difference of Gaussians in each case. Third panel shows the resulting masks and green circles as same as in top panel. Bottom panel shows the relation between real size of the feature and the size found by the algorithm.

In the text