{'fraction_non_alphanumeric': 0.02448022802040463, 'fraction_numerical': 0.02229592040596156, 'mean_word_length': 4.653765060240964, 'pattern_counts': {'":': 0, '<': 6, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 1, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
High-frequency wave power observed in the solar chromosphere with IBIS and ALMA

We present observational constraints on the solar chromospheric heating contribution from acoustic waves with frequencies between 5 and 50 mHz. We utilize observations from the Dunn Solar Telescope in New Mexico complemented with observations from the Atacama Large Millimeter Array collected on 2017 April 23. The properties of the power spectra of the various quantities are derived from the spectral lines of Ca II 854.2 nm, H I 656.3 nm, and the millimeter continuum at 1.25 mm and 3 mm. At the observed frequencies the diagnostics almost all show a power law behavior, whose particulars (slope, peak and white noise floors) are correlated with the type of solar feature (internetwork, network, plage). In order to disentangle the vertical versus transverse plasma motions we examine two different fields of view; one near disk center and the other close to the limb. To infer the acoustic flux in the middle chromosphere, we compare our observations with synthetic observables from the time-dependent radiative hydrodynamic RADYN code. Our findings show that acoustic waves carry up to about 1 kW m$^{-2}$ of energy flux in the middle chromosphere, which is not enough to maintain the quiet chromosphere, contrary to previous publications.


INTRODUCTION
Balancing the radiative losses of the non-magnetic chromosphere requires an energy input of about 4 kW m −2 (Withbroe & Noyes 1977). The two most widely accepted theoretical frameworks for chromospheric heating are the same as for the corona: wave dissipation or ubiquitous small scale magnetic reconnection (Carlsson et al. 2019). There is evidence for both mechanisms being at work in the chromosphere, but definitive observations constraining their relative importance are still elusive. In this paper we concentrate on quantifying the contribution from compressive waves in the chromosphere which produce measurable Doppler shifts in chromospheric diagnostics.
The first to suggest that the chromosphere can be kept in its thermal state by dissipation of acoustic waves were Biermann (1946) and Schwarzschild (1948). They suggested that the convective overshoot at the boundary of the the upper convective zone and the photo-Corresponding author: Momchil Molnar momo@nso.edu * DKIST Ambassador sphere drives acoustic waves with a fairly broad range of periods. Those waves with frequencies above the acoustic cutoff frequency (roughly 5 mHz or 200 seconds) can propagate upward into the chromosphere (Bray & Loughhead 1974). As the acoustic waves move into higher layers in the solar atmosphere, they find a strongly decreasing plasma density while the temperature, and hence sound speed, remains (almost) constant, which results in amplitude growth and wave steepening. Waves that steepen into shocks can dissipate their energy and supply the heat needed to maintain the chromospheric plasma in its basal state (Carlsson & Stein 1992). Intermittent shocks would ionize the hydrogen and helium in the chromosphere out of equilibrium and maintain the chromospheric ionization state away from statistical equilibrium due to the long recombination timescales (Carlsson & Stein 2002).
The magnetic field in the chromosphere allows for the existence of magnetosonic wave modes and the similarity between the sound and Alfvèn speeds allows the easy conversion between them (Cally & Goossens 2008). The multitude of magneto-acoustic wave modes and their numerous damping mechanisms expands the possible wave propagation scenarios, but does not alter the basic premise of the theory of how energy is being transported from the convection zone to the chromosphere. For a recent review on the subject of magnetosonic waves, the reader should see Jess et al. (2015).
The observational evidence constraining the energy contribution of wave heating in the chromosphere has been inconclusive. Wunnenberg et al. (2002), using Fabry-Perot imaging spectroscopy of the Fe I 543.4 nm line, inferred around 0.9 kW m −2 acoustic flux from waves with periods between 50 s and 100 s (10-20 mHz) at height of 600 km above the photosphere. They estimated that the large extent in height of the velocity response functions reduced their observed wave amplitudes by a factor of three, and arrived at an actual acoustic flux of 3 kW m −2 . On the other hand, Fossum & Carlsson (2005) used TRACE intensity-only observations of the 160 nm UV continuum (sampling the upper photosphere at height of 450 km), coupled with self-consistent simulations of the solar atmospheric oscillations, to derive an acoustic wave flux in the frequency interval 5-28 mHz of 0.4 kW m −2 . This would not be sufficient to sustain the non-magnetic chromosphere, but other authors (e.g., Cuntz et al. 2007;Wedemeyer-Böhm et al. 2007) have argued that the limited angular resolution of the TRACE observations and other model assumptions lead to significant underestimation of the wave flux in that analysis.
Measurements of velocities in two photospheric spectral lines with higher temporal and spatial resolution by Bello González et al. (2009González et al. ( , 2010 showed the presence, after taking into account the width of the velocity response functions, of significant acoustic flux (up to 3.8 kW m −2 ) in the middle photosphere between 5 and 15 mHz. More recent work by Abbasvand et al. (2020a,b) utilizing observations of the chromospheric Ca II 854.2 nm and H I Balmer-α and Balmer-β lines with a scanning spectrograph reached similar conclusions of about 5 kW m −2 flux in the chromosphere. However, the latter authors did not account for the width of the velocity response function (Mein & Mein 1980), which makes their flux estimation a lower bound for the actual wave energy flux.
Another still poorly understood aspect of the wave heating theory is the contribution from high-frequency waves, above 30 mHz. The sparsity of observations in this regime has been due to the difficulty of obtaining high-temporal-cadence spectral information of the chromosphere at the high spatial resolution required to resolve the small-scale chromospheric structures. In one of the few studies in this regime, Hansteen et al. (2000) has shown intermittent wavelet power up to 50 mHz in upper-chromospheric and transition-region lines taken with SUMER on the SOHO spacecraft.
Some of the aforementioned studies take into account the attenuation of the observed wave amplitudes when the contribution to the observed spectral signal becomes similar to or greater than the wavelength of the propagating waves (Mein & Mein 1980). Consideration of this effect is essential for inferring the flux, as it has been estimated to be a factor of between 2 to 10 for frequencies between 5 and 50 mHz (Wunnenberg et al. 2002;Bello González et al. 2009). The latter authors used a static semi-empirical 1D model to calculate the wave response in photospheric lines using a perturbative approach with sinusoidal waves. This approach is not similarly applicable in the chromosphere, where the waves are generally not sinusoidal and could be strongly affected by radiative losses. A more realistic approach is undertaken by Fossum & Carlsson (2005) and Wedemeyer-Böhm et al. (2007) who use time-dependent hydrodynamic simulations to infer the actual wave attenuation including the effect of radiative wave damping. We show in Section 5 that this approach is better as it naturally explains the high frequency signal in our observations. Furthermore, we argue that modeling based on semi-empirical 1D atmospheres might be overestimating the contributions from the high-frequency waves.  showed that Doppler diagnostics in the chromosphere have a power-law behavior from the acoustic cutoff frequency out to 20 mHz. It is not understood whether this trend is due to the true distribution of acoustic oscillations at these frequencies, or the result of a frequencydependent attenuation of the chromospheric signal.
This paper presents observations of high frequency wave Doppler velocity signal in the solar chromosphere. We obtained a data set of cotemporal observations with the Interferometric Bidimensional Imaging Spectrograph (IBIS) instrument at the Dunn Solar Telescope (DST) and with the Atacama Large Millimeter Array (ALMA) in a sparsely explored temporal regime up to 50 mHz. In Section 3 we present evidence for the presence of high-frequency power in the velocity measurements of the Hα and Ca II 854.2 nm lines as well in the ALMA brightness temperatures. To infer the wave energy fluxes from our observations we model the propagation of acoustic waves throughout the solar atmosphere with the RADYN code in Section 4. The results from the modeling are presented in Section 5 and the discussion of our results and the conclusions are summarized in Section 6.

OBSERVATIONS AND DATA PROCESSING
We obtained coordinated solar observations with ALMA (Wootten & Thompson 2009) and the DST (Dunn 1964;Dunn & Smartt 1991) on 2017 April 23. These observations were part of ALMA Project ID 2016.1.01129S/Cycle 4. See Molnar et al. (2019) for some of the initial results from this observational campaign.
Two separate fields of view (FOV) were observed on 2017 April 23 and they are referred chronologically as FOV 1, near the limb, and FOV 2, close to disk center. The targets in both cases were regions of magnetic network or plage, but they were observed at different inclination angles to the solar surface, which is essential for the discussion in Section 3.3. Context images of FOV 1 and 2 are provided from SDO/AIA (Lemen et al. 2012), SDO/HMI (Schou et al. 2012), andIRIS (De Pontieu et al. 2014) in Figure 1. Further analysis of the IRIS dataset is left for a forthcoming publication.
FOV 1 was observed with the DST at 13:50-15:14 UT at solar coordinates E 66.2 • , S 09.8 • , at an inclination of µ = 0.41 (the cosine of the angle between the line of sight and the solar surface normal). FOV 1 was cen-tered on the trailing edge of NOAA Active Region (AR) 12653. Based on the context imaging (top panel of Figure 1) FOV 1 contains some active region plage as well as internetwork regions with little magnetic field. The inclined viewing angle results in more confusion among features in the field of view because of projection effects of predominantly vertical features and longer integration along the line of sight.
FOV 2 was observed with the DST at 15:15-18:19 UT at solar coordinates E 4.9 • , N 10.9 • at an inclination of µ = 0.98. FOV 2 was centered on the leading edge of the active region NOAA AR 12651. There is a magnetic concentration in the center of the field, surrounded by a largely field-free internetwork area, especially in the southern portion of the field. There is a region of plage in the northwestern corner of the field and the leading edge of the penumbra/superpenumbra on the northeastern corner.

DST observations
The DST took observations on 2017 April 23 between 14:00 UT and 18:40 UT, under conditions of excellent to good seeing. The instrument setup included IBIS (Cavallini 2006;Reardon & Cavallini 2008), the Facility In-fraRed Spectrograph (FIRS, Jaeggli et al. 2010) and the Rapid Oscillations in the Solar Atmosphere instrument (ROSA, Jess et al. 2010), which provided thorough coverage of key spectral lines in the optical and the near-IR parts of the spectrum. All of the instruments were fed by the high-order adaptive optics system (Rimmele 2004). None of the instruments were run in polarimetric mode, as high temporal cadence was the priority for this study.
IBIS observed the spectral lines of H I Balmer-α 656.3 nm (Hα), Ca II 854.2 nm and Na I D 1 589.6 nm with an average plate scale of 0.096 pixel −1 . Each scan of a single spectral line scan took between 3 and 4 seconds with an overhead of about 1.5 seconds for changing the prefilters. We utilized two different scanning strategies during the observations. At the beginning of each ALMA observing block at a given pointing and frequency band, we ran "fast" scans of Hα and Ca II 854.2 nm. Each scan consisted of 25 and 21 wavelength points, respectively and the average profiles from those scans are shown in Figure 2. Scanning a single line avoids the overhead for changing IBIS prefilters, resulting in a cadence of about 3.5 seconds for a single spectral scan (the precise cadences are listed in Table 1). This observing strategy was adopted to closely match the ALMA twosecond sampling rate and allow us to study the high frequency wave regime (with Nyquist sampling of 130 mHz). These scans were intended to capture the fast dynamics of the chromosphere, at high temporal frequencies rarely explored with wide-field, bi-dimensional spectroscopy. Table 1 contains a summary of the single Figure 2. Average spectral profiles from the fast scans with IBIS in the Hα (656.3 nm) and Ca II 854.2 nm lines. The line cores were sampled more densely to allow for more accurate determination of the line properties derived from the line cores (velocity, intensity and width). The blue shaded regions are used for the line core fitting with a parabola after a resampling on a finer (0.005 nm) grid.
line ("fast") scans with IBIS and the ALMA time series used in this study.
Following these fast scans, we ran longer durations of a "standard" repeating cycle of all the three spectral lines. The standard scans consisted of 90 total spectral points among the three lines and had duty cycle of about 15 seconds. This type of scan thoroughly covers the solar atmospheric layers from the photosphere through the middle chromosphere, making it well suited for studying the propagation of wave energy through the lower solar atmosphere.
To remove the seeing distortions in our narrowband images we relied on destretching the cotemporal broadband images (from the white-light IBIS channel) to HMI white-light images using the sub-aperture crosscorrelation method introduced by November & Simon (1988). We applied the destretch maps to the cotemporaneous narrowband images and the resulting stability of the destretched images was on the order of one IBIS pixel (∼ 0.1 ). The wavelength samplings for the fast scans of Hα and the Ca II 854.2 nm lines are presented in Figure 2. The lines were sampled non-equidistantly to ensure better coverage of the line core, which is used for the velocity, intensity and width measurements. The maps of the observed line core intensity and width in Hα and Ca II 854.2 nm for both FOVs observed with IBIS are presented in Figure 3.
Following Cauzzi et al. (2009), we measured the width of the chromospheric line cores as the separation of the two wing positions at which the intensity reached an intensity level halfway between the core intensity and the wing intensity at a defined wavelength offset from the local core position. The wing offsets used are ± 0.1 nm for the Hα line and ± 0.06 nm for the Ca II 854.2 nm data. The resulting Hα line width maps are shown in the middle column in Figure 3. These line widths are thought to be related to the temperature of the emitting plasma under chromospheric conditions Leenaarts et al. 2012;Molnar et al. 2019).
We measured several different velocity signatures from our spectra. After sampling onto an evenly spaced wavelength grid with 0.005 nm sampling, we fitted the line minimum position and intensity with a parabola, in order to determine the Doppler velocity, which is related to the velocity at the τ =1 for the line core in the atmosphere. The fitting was done on an interval of ±0.02 nm and ±0.015 nm for around the minimum position for for Hα and Ca II 854.2, respectively. This corresponds to approximately four (six) of the originally sampled points for Hα (Ca II 854.2), indicated in Figure 2 as the blue region.
We also calculated the center of gravity (COG) velocity of over ± 0.12 nm wavelength region for the Hα line and over a ± 0.105 nm for the Ca II 854.2 nm line This velocity measure takes into account the whole line profile and might carry some information about the photospheric velocity field in the case of Hα (Socas-Navarro & Uitenbroek 2004). Finally, we used the same bisector calculation described above to determine a bisector shift for both lines at 50% level.

ALMA observations
ALMA became available for solar observations in Cycle 4 (2017) after extensive testing (Phillips et al. 2015). The continuum wavelength bands available for solar observations were Band 3 (3 mm/100 GHz) and Band 6 (1.25 mm/240 GHz) which are expected to sample the high and middle chromosphere (Wedemeyer et al. 2016). The continuum radiation in these wavelengths forms in the chromosphere under Local Thermodynamic Equilibrium (LTE) conditions, as the main source of opacity is due to free-free processes which makes the source function to be locally determined by the plasma temperature. The ALMA intensity can be interpreted as brightness temperature under the Rayleigh-Jeans limit (Wedemeyer et al. 2016). Hence, ALMA is a valuable tool to study the thermal structure of the chromosphere (some recent results include Shimojo et al. 2017;Loukitcheva et al. 2019;da Silva Santos et al. 2020).
ALMA is an excellent instrument for studying highfrequency waves in the solar atmosphere due to its fast sampling cadence of 2 seconds, spatial resolution better than 2" (depending on the array configuration), and direct sensitivity to electron temperature. However, it is important to remember that the opacity scale of the ALMA radiation is determined through the local electron density (and the ionization state of the plasma), which is thought to be far from thermodynamical equilibrium (Carlsson & Stein 2002). The time-varying opacity can complicate the interpretation of the timeseries of temperature brightness measurements (Molnar et al. 2019).
The ALMA observations discussed in this paper were obtained in configuration C40-3 and their reduction is described in Molnar et al. (2019). The millimeter observations are obtained in approximately 10-minute blocks interspersed with several minutes of off-target calibrations. For FOV 2, we analyzed temporal blocks 1 (17:19-17:29 UT) and 2 (17:31-17:41 UT) from Band 3 (also used in Molnar et al. 2019) and blocks 1 (15:53-16:01 UT) and 2 (16:03-16:13 UT) from Band 6. An ALMA data summary is presented in Figure 4. We chose these particular ALMA observing blocks because they are the closest in time to our high-cadence IBIS observations of FOV 2 listed in Table 1. The relative positions of the ALMA observing fields to the IBIS observation regions are shown in Figure 1 as the colored circles. The useful regions of the FOVs of the ALMA data are about the size of the respective circles in Figure 1: 60 for ALMA Band 3 and 20 for Band 6.

Solar feature classification
To study the typical wave characteristics in different regions of the solar atmosphere, we partitioned the solar surface in the field of view into five different classes of features that represent general chromospheric structures: penumbra, internetwork, fibrils, network and plage. We distinguished the regions following a methodology based on the properties of the Hα and Ca II 854.2 nm spectral lines.
For FOV 2, we first categorized the penumbral region by its proximity to the leading sunspot in NOAA AR 12651 and its low Hα line width (smaller than 0.1 nm) and low Doppler velocity fluctuations relative to the rest of the field. Secondly, we distinguished between the plage and the network regions, which are brightest in FOV 2. We can clearly distinguish between the two by the amount of magnetic flux and the intensity in the Ca II 854.2 nm line core intensity. We assign a pixel to be a plage region if the photospheric magnetic field strength (from HMI) is above 1200 G and the Ca II 854.2 nm line core intensity is more than 60% above the mean intensity of the whole FOV. Pixels were categorized as network if their magnetic field was above 250 G and their Ca II line core intensity is below the 60% intensity threshold. We labeled the regions with lowest intensity in the core of Ca II 854.2 nm and with the lowest Hα line core width (less than 0.105 nm) as internetwork or fibrils. We distinguished between the fibrils and the internetwork regions by the ratio of the relative power in the 3 min to 5 min power in the Doppler velocity power spectra as suggested by Vecchio et al. (2007), where the internetwork has a ratio greater than one and the fibrils have a ratio less than one between those two frequency windows.
Each of these criteria resulted in binary masks, which we blurred with a Gaussian filter (with a standard deviation of 40 pixels) to smooth the boundaries of the regions and avoid holes. In cases where the smoothing caused masks for separate classes to overlap, we chose the darker of the two classes to define that pixel. The blurring insured, for example, that isolated magnetic elements in the internetwork were smoothed out and the classified regions were largely contiguous. The resulting mask for FOV 2 which we employed throughout the rest of the paper to distinguish the different regions of the solar surface is shown in Figure 5.
We use a similar masking algorithm for the FOV 1 observations, but using only three types of solar features because there is no penumbra in this target and we can-  not easily distinguish between the internetwork and the fibrils at such highly inclined projections. We take into account in our analysis that FOV 1 is mostly a plage region and hence label all of the "quieter" regions as fibrils. We increased the Hα core width cutoff for fibrils to 0.12 nm, as the average Hα line width increases closer to the solar limb. We distinguished between plage and network by the photospheric magnetic field strength -if the magnetic field was above 1200 G (as for FOV 2) then we classified the pixel as plage, otherwise it was classified as network. We also used a Ca II 854.2 nm line core intensity threshold that was twice the average intensity of the whole FOV for distinguishing between plage and network The resulting mask for FOV 1 is presented in Figure 6.

OBSERVED HIGH-FREQUENCY POWER IN THE SPECTRAL DIAGNOSTICS
In this Section we present the high-frequency observations derived from the Power Spectrum Density (PSD) of the observed chromospheric diagnostics (intensities, velocities, etc.) described in Section 2. The PSDs are calculated as the absolute value of the Fourier amplitudes of the observed signals, where we subtracted the mean of the observed time series and then apodized them with a Hamming window with 3% of the total length of the observations on each side.

Observed high-frequency power in chromospheric spectral diagnostics
The top four rows of Figure 7 present a summary of the power spectra of the chromospheric spectral diagnostics derived from Hα (data series starting at 15:18:57 UT) and Ca II 854.2 nm line (data series starting at 15:13:21 UT) IBIS observations of FOV 2. The pixels analyzed for this plot were taken from the central 12 × 12 regions of FOV 2 (the green rectangle in Figure 3) to ensure that the same solar region is observed with IBIS and ALMA. Furthermore, the selected area (green square) spans only a network region, which simplifies Figure 6. Mask of the different solar regions in FOV 1 using the approach described in Section 2.3. No distinction between internetwork and fibrils was made for this FOV due to the large viewing angle and preponderance of magnetic flux.
the analysis to a single type of source region. The red dotted line is the mean of the distribution at each frequency and the lighter blue lines are individual power spectra of each pixel. The frequency resolution of the power spectra is on the order of δν ≈ 2 -3 mHz as the length of our data series is about 10 minutes and the cadence is on the order of 3.3 seconds (see Table 1 for the details of the observations). We do not clearly observe the low frequency power around the five-minute oscillation window with high resolution due to the short temporal extent of these fast scans.
The power spectra show a strikingly similar power law above 7 mHz in all observed spectral diagnostics, limited at the higher frequencies by white-noise floor. The level of the white noise varies for the different spectral diagnostics. Our results extend previous work by  to higher frequencies where they agree in the low frequency part (7-20 mHz) of the power spectrum. The white noise which dominates the high frequency signal is likely due primarily to the photon noise from the measured signal. In Appendix A we describe a detailed estimation of the effect from photon noise on the white-noise floor of the line core velocity measurements of the Ca II 854.2 nm data. This additional contribution to the PSD is important to characterize in order to properly quantify only the solar contribution to the integrated power.
The observed and simulated white-noise level distributions are shown in the Figure 8 and summarized in Table 2. In both cases, we see that the profiles with higher Ca II 854.2 nm line core intensities (plage, network) tend to have higher noise levels compared to deeper profiles (with more pronounced core minima). The distributions of the simulated noise levels have somewhat broader tails, indicating a possible overestimation of the noise in our model in extreme cases. However, the median values (dashed vertical lines) of the simulated distributions are all in good agreement with the observed ones. Hence, we believe that we can well characterize the noise floor separately for different classes of solar structures in our dataset.
We employ the same data analysis approach for the brightness temperatures from the ALMA Bands 3 and 6 data. The PSD of the brightness temperature from the same central 12 × 12 of the ALMA FOVs are presented in the last row of Figure 7. The red lines again indicate the average of the distribution and the light blue lines are the power spectra of individual pixels.
The power law behavior of the ALMA brightness temperature PSD (seen before in Nindos et al. (2020)) extends to the white-noise floor at about 100 mHz in the Band 3 data. In the Band 6 observations we do not see clearly the white-noise floor. We do not observe the 3 minute (5 mHz) oscillations clearly, like in Patsourakos et al. (2020). This might be due to presence of magnetic elements in the observed region as suggested by Jafarzadeh et al. (2021), or due to the limited frequency resolution of our data of about 2 mHz. However, the observed region is mostly covered under the magnetic canopy where we notice the bright fibrils (in ALMA wavelengths) dominating our field of view (see Figures 1 and 3).

Properties of the observed power spectra
We concentrate on the Ca II 854.2 nm velocity power spectra as they are the most reliable chromospheric velocity diagnostic we have, as the Hα line synthesis results depends strongly on the full 3D radiative transfer solution (Leenaarts et al. 2012). Using the feature classifications derived in Section 2.3 and shown in Figure 5, we calculated the average Ca II 854.2 nm line core Doppler velocity PSD for every type of region of the solar surface described above. The average PSDs are presented in Figure 9. The plage and internetwork have higher amounts of velocity oscillation power compared to the network and the quieter regions (fibrils and penumbra). However, the plage has a significantly higher white-noise Figure 7. Power spectra of different spectral chromospheric diagnostics. In the first four rows, the left columns correspond to spectral properties derived from the Ca II 854.2 nm line (Ca II IR) and the right column from the Hα line observations. The chromospheric diagnostics which each row is derived from are: First row: COG velocity; Second row: line center velocity; Third row: intensity of the line core (where a.u. stands for arbitrary units); Fourth row: equivalent width. The last row corresponds to PSDs derived from ALMA Band 3 (left row ) and Band 6 (right row ) intensity observations. The blue lines represent individual power spectra and the red dotted lines show the mean in each frequency bin. floor (comparable to the network one) than the one seen in the internetwork. The internetwork has high amount of wave flux in the 5-20 mHz interval but a lower whitenoise floor -similar to the one in the quieter regions (fibrils, penumbra). Figure 10 shows the spatial maps for FOV 2 of the amount of oscillatory power at each pixel integrated between 5 and 50 mHz in the line core intensity and Doppler velocity measured from the Ca II 854.2 nm and Hα lines as well the oscillatory power derived from the ALMA brightness temperatures. Throughout the paper when we refer to oscillatory power we mean the integrated power spectral density (PSD) in the specified frequency range. When estimating PSDs, we always subtract a white noise estimate based on the mean of the last 25 (high) frequency bins. We see that properties of similar diagnostics (velocities, intensities and widths) derived from different spectral lines have similar distributions. The ALMA Band 3 temperature fluctuations map (top right) correlates with the spatial distribution of the PSDs in the optical diagnostics, but the Band 6 one (bottom right) does not appear to resemble its optical counterparts. This might be due to time-varying changes in the height of formation of the ALMA intensities as well as the lower angular resolution of the ALMA observations or the very limited FOV of the ALMA Band 6 observations.
We calculated the power law slope describing the observed power spectra with least-χ 2 fitting between 10 and 40 mHz. Using the same solar feature mask, we calculated the distributions of the slopes in the Ca II 854.2 nm line core velocity in the different regions of the solar surface and the results are presented in Figure 11. The distributions of the slopes from the observed (raw) data are presented in Panel a) of Figure 11. We see that plage and network regions have shallower slopes than the fibrils, internetwork and penumbral regions. The slope distributions are overlapping and form a continuous transition between all of the different types of solar features, but their progression does follow the trend of mean Hα line-widths for those solar features (excluding the penumbra). The quantitative comparison of the power slope distributions is listed in Table 2, along with some other quantities derived from the observed power laws. The quoted regions of uncertainty are the 10th/90th percentiles of the cumulative distributions. We note that the power-law slopes we find for the network, fibril, and internetwork regions are very comparable to those found by  for the 5-20 mHz interval. To infer the true slope of the vertical velocity field in the solar atmosphere, we applied the velocity attenuation coefficient described in Section 4.3. The resulting corrected velocity power law slope distributions are presented in Panel b) of Figure 11 and their medians are summarized in Table 2. The relative order of the solar features is preserved, but the resulting distributions of slopes are shallower as the compensation for the wave attenuation makes the power laws less steep. The resulting values of the power law slopes of the corrected data are roughly between −2 and 0. These values are shallower than the ones expected from the Lighthill-Stein turbulence theory (Ulmschneider et al. 1996), which predicts power law slopes between −3.5 and −3. Eulerianbased treatments of turbulence in the solar atmosphere give slopes between −2.4 and −1.3, which are closer to our observed values (Rubinstein & Zhou 2002). Hence, our observations favor the Eulerian approach in treating the turbulent time correlations in the solar atmosphere. However, the penumbra, the network, and the plage regions which exhibit higher power law slopes than either theory predicts and require further investigation, which we leave for a following study. Figure 12 shows the correlation between the uncorrected power-law slopes and the total amount of power in each spectral diagnostic, with pixels with higher integrated power having steeper slopes. Since the PSDs are steeper, the higher power is a consequence of those regions tending to have higher overall oscillatory power in the lower portion (5-10 mHz) of the frequency interval.
We partitioned the Band 3 field of view into three regions based on the mean observed 3 mm brightness temperature: T b < 7500 K; 7500 K < T b < 10000 K; and T b > 10000 K. The resulting segmentation mask is presented in Figure 13. We used this segmentation approach because the classification described in Section 2.3 did not necessarily correspond well with the ALMA brightness temperature features. In particular, the fibril regions exhibit brightness temperatures comparable to those seen in the network regions. The average power spectrum of the ALMA Band 3 brightness temperature for each of these different solar regions are presented in Figure 14. The power spectra share similar slopes, in-dependent of the region, with a value of -1.63 ± .07 in the region of 10 to 50 mHz. The brighter regions have slightly higher oscillatory power at frequencies below 10 mHz. We did a similar analysis of the Band 6 data (using the same temperature-based segmentation based on the Band 3 data), but did not find any significant differences in the slopes or total power among the different (temperature discriminated) regions of that FOV. This might be due to the rather small FOV of Band 6 (20 ) compared to Band 3 (60 ). Figure 15 shows the correlation between brightness temperature and the relative (normalized to the mean) brightness temperature fluctuations in the ALMA data between 5 and 50 mHz. ALMA Blocks 2 were used for both bands in the figure. In order to avoid the attenuation of the sensitivity farther away from the center of the beam, we used only the central 60 of the Band 3 FOV and the central 17 for the Band 6 FOV. The red lines show the median trend of the histograms. We can see a clear correlation between the brightness temperature and oscillatory power in the Band 6 data and an almost non existent one in the Band 3 data. The relative temperature fluctuation power is related to the amount of compressive wave flux as shown in Section 4.2.

Velocity oscillations at different viewing angles
Observing the solar atmosphere at different inclinations provides a way to disentangle the longitudinal from the transverse velocity oscillations, as these two components are differently projected into line-of-sight Doppler shifts. This allows us to statistically disambiguate between the transverse and the longitudinal oscillations if we observe similar solar regions at (a minimum) two different viewing angles. FOV 1 was observed at an incidence angle of 66 degrees or µ = 0.41, close to the east solar limb.
Assuming mostly vertical magnetic field orientation to the solar surface, the observed velocity oscillatory power is to be composed of not only the longitudinal velocity oscillations v (angle brackets stands for the average root mean square value over time of the quantity), and a transverse (Alfvénic-like) v ⊥ component. The perpendicular component v ⊥ is of special interest for constraining coronal heating models, as those waves are expected to propagate readily throughout the chromosphere and the into the corona (Cranmer & van Ballegooijen 2005). In the case of perpendicular to the solar surface observations, close to disc center as in the case of FOV 2, we will detect only v as the Doppler velocity. Knowing v from the observations of FOV 2 (close to disc center), we can calculate v ⊥ from the observed velocity oscillations v 2 obs (θ) (at an incidence angle θ) in FOV 1. The observed velocity oscillations v 2 obs (θ) can be written into their components as (assuming that v and v ⊥ are not correlated): where µ = cos θ, the cosine of the incidence angle of the observation.
The average power spectrum of the Ca II 854.2 nm FOV 1 velocity (shown in Figure 16) exhibits lower oscillatory power compared to FOV 2 and has higher whitenoise floor compared to FOV 2, where the averaged FOV 1 data are shown as the black dots in Figure 9. This is further illustrated in Figure 17 where the velocity fluctuation power in FOV 1 and 2 are compared. The red curves in Figure 17 show the FOV 1 velocity fluctuation power between 5 and 20 mHz for the different solar regions as segmented in Figure 6; the yellow curves show the distributions for the solar regions for FOV 2. Under the assumption that the velocities detected in FOV 2 are all due to longitudinal displacements (since the observed region was close to disc center at µ = 0.98), we projected the observed distributions to what would be observed at an inclination angle of µ = 0.41 by proportionally reducing the velocity magnitudes (green curves). For all three types of solar structures (fibrils, network and plage), the amount of oscillatory power in the projected FOV 2 distributions is smaller than what is actually observed at that inclination in the FOV 1 data. Therefore, we believe that in FOV 1 we are observing not just projected vertical velocities (v ), but also an additional transverse (v ⊥ ) component. Based on Equation 1 and taking the median of the velocity oscillation distributions in Figure 17 as representative of the averaged oscillation power, we calculate Figure 13. The mask distinguishing the different regions of the solar surface in for the ALMA Band 3 FOV. The mask is based on the brightness temperature T b split into three categories -T b < 7 500K, 7 500 < T b < 10 000K, and T b > 10 000 K.    the magnitude of the transverse oscillations v 2 ⊥ for the different solar regions. The resulting velocity components for the different solar regions are presented in Table 3. We can see that the plage region has the highest longitudinal oscillations and lowest transverse components. The network and fibrils have similar values for the transverse oscillation power, while the network has slightly higher longitudinal oscillation power. On average, the value of transverse velocity component v ⊥ is a few times smaller than the longitudinal one v . This has been suggested in previous modeling work (for example Cranmer et al. (2007)).
Further observational studies could make use of samples at more values of µ to better confirm and constrain the contributions of the two components. In addition, a more detailed treatment of radiative transfer effects at inclined viewing angles should be made. Observations of more homogeneous solar regions (such as quiet Sun or network) would provide more assurance that similar structures were being compared and provide more accurate estimates for v 2 ⊥ .

MODELING THE OBSERVATIONAL SIGNATURES OF ACOUSTIC WAVES IN THE CHROMOSPHERE
Optically thick diagnostics, such as the Hα and Ca II 854.2 nm lines or the millimeter continuum used in this study, sample a wide range of heights in the chromosphere. Further, the interval sampled by these diagnostics are changing dynamically during the continuous evolution of the chromospheric properties which are significantly out of thermodynamical equilibrium (Carlsson et al. 2019).
The observed Doppler velocity response of a longitudinal wave-like perturbation in the solar atmosphere, with a wavelength comparable to the extent of formation region of the spectral line, will be attenuated due to the inherent mixing of signals from different phases of the oscillatory fluctuation (Mein & Mein 1980). Modeling is needed to estimate the extent of this wave signature attenuation and infer the true wave flux corresponding to the measured oscillatory amplitude.

Propagation of acoustic waves in the solar atmosphere
The dispersion relation for acoustic waves in the solar atmosphere is (following the derivation in Bray & Loughhead 1974): where ω = 2πν is the angular frequency; k h and k z are the horizontal and vertical wave numbers; ω ac = γg/(2c 2 s ) is the cutoff frequency below which no acoustic waves can propagate in the atmosphere. Under typical chromospheric conditions, the speed of sound is c s = 7 km/s and ν ac = 5.2 mHz for an adiabatic index of γ = 5/3; ω BV = √ γ − 1 g/c 2 s is the Brunt-Väisälä frequency which determines the lower cutoff frequency below which gravity waves can propagate.
The energy flux F ac carried by an acoustic wave with a velocity amplitude squared per unit frequency P (v 2 w ) in the frequency interval of ν 0 to ν 1 is: where ρ(z) is the density at height z, v gr (ν) is the group velocity. The group (propagation) velocity is the derivative ∂ω/∂k z from Equation (2) and in our case it equals the sound speed times a factor of order unity above the acoustic cutoff frequency and zero below it: Compressive waves will also show a temperature fluctuation, which is proportional to the energy flux of the wave. We can substitute in Equation (3) the wave velocity fluctuation v w with the magnitude of the wave temperature fluctuation (Bray & Loughhead 1974): where Γ 1 is the adiabatic exponent describing how pressure responds to compression; δT is the wave temperature fluctuation and T 0 is the mean temperature of the atmosphere. Based on Equation (5), we can substitute in Equation (3) to get the final expression for the wave energy flux expressed via the relative temperature fluctuations squared per unit frequency P (δT /T 0 ): where we integrate again over the frequency interval above the acoustic cutoff and below the frequency at which the white noise floor dominates our signal. ALMA is an ideal tool for the detection of the wave temperature fluctuations as it measures directly and linearly the plasma temperature in the chromosphere, with the observed brightness temperature equal to the contribution-function-weighted mean temperature over the formation interval. This is a supplementary measurement to the direct velocity measurement with IBIS in the Ca II 854.2 nm line.

Synthesizing wave observables with RADYN
To model the line formation in the presence of waves we need time-dependent models of the solar atmosphere which incorporate all of the necessary physical processes of wave propagation (including the optically thick radiative transfer effects). We used the RADYN radiative hydrodynamic code (RHD) (Carlsson & Stein 1992, 1995, 1997Abbett & Hawley 1999;Allred et al. 2015) to self-consistently model the propagation of high frequency acoustic waves in the chromosphere. This 1D code solves the hydrodynamic equations coupled with the radiative transfer equations in non-local thermodynamical equilibrium (NLTE). RADYN supports subphotospheric velocity drivers defined by the user and treat time-dependent NLTE ionization of the primary atomic species (H, He and Ca). RADYN performs the radiative transfer calculations to generate time-dependent synthetic line profiles consistent with the wave dynamics. Based on these synthetic observables, we are able to interpret our data in terms of the realistic solar plasma parameters and estimate the amount of flux carried by the acoustic waves in the chromosphere.
In our numerical setup, we use a 191 point atmosphere with 6-level atom models for Hydrogen and Calcium and 9-level atom for He -including singly and doubly ionized states. We use an open (transmitting) upper boundary condition where the corona is maintained at constant temperature 1 MK at height of 12 Mm and a bottom piston boundary condition at fixed temperature of 5944 K. The piston bottom boundary allows for driving vertically propagating acoustic waves with arbitrary properties which are defined by the user.
We ran a grid of models with an increasing amount of wave energy being injected through the bottom boundary condition. These models are used consistently throughout the paper with the name format of model XXXX, where the increasing numerical factor stands for a stronger bottom boundary driver. We used the bottom boundary vertical velocity drivers presented in Fossum & Carlsson (2006), scaling their amplitude by a multiplicative factor to achieve the desired amount of input wave energy. These drivers specify a power spectrum of sub-photospheric velocities, at a range of oscillatory frequencies from 1 to 50 mHz, that have different relative amplitudes. To synthesize the velocity time series for the bottom boundary condition, we used the power spectrum of the driver initialized with random phases. The vertical velocity PSDs of the used bottom boundary drivers are presented in Figure 18. The resulting acoustic flux in the chromosphere from those runs are presented in Figure 19. We have calculated the acoustic fluxes from Equation (3), where we have interpolated the models on a 4 000 point equidistant height grid to remove the movement of the grid points during shock passages. To compute the average acoustic fluxes, for each grid point we calculate the mean of the sound speed and the plasma density and filter the vertical velocity between 5 and 50 mHz. We have chosen time series starting about 10 minutes after the Figure 18. The vertical velocity power spectrum of the bottom boundary condition of the different models we ran. The overall shape of the PSD is taken from Fossum & Carlsson (2006) and scaled by a multiplicative factor and then supplied as a bottom boundary (vertical velocity) condition in the RADYN runs. The PSDs presented here are the Fourier transforms of the actual 3 second sampled RADYN run bottom boundary velocity. initialization of the simulation, when a steady state is reached. We have further excluded models with higher acoustic fluxes above model 3000 as the Ca II 854.2 nm line core goes into emission and do not reproduce the observations.
To compare the high-frequency signal in the synthetic observables and in the IBIS observations we took the simulated Ca II 854.2 nm line profiles from the RADYN runs at 0.5 second intervals and averaged them to the temporal and spectral sampling of the IBIS instrument. We did apply the spectral PSF of IBIS, but this did not change the final results significantly as its spectral resolution is very high (R > 200,000) for the Ca II 854.2 nm line. To measure the synthetic velocities, we used the same methods (i.e. line center fitting) used for the real data processing. 4.3. Estimating the acoustic flux in the Ca II 854.2 nm data based on RADYN Given the attenuation of the wave amplitude discussed in Section 4.1 we need to find a relationship between the observed Doppler velocities, v obs , and the true wave amplitudes v w . One approach would be to discretize the acoustic flux in the observed bins and then correct it for the wave amplitude attenuation for each frequency bin separately (as in Bello González et al. (2009)): where the coefficient T (ν i ) quantifies the attenuation of the wave velocity amplitude at frequency ν i due to the extended formation height of the spectral line. In this case the quantity v 2 obs (ν i ) is the observed amount of velocity oscillation power per frequency bin and has units of velocity squared. T (ν i ) is defined as the ratio of the observed wave amplitude to the real (physical) wave amplitude in the middle of the formation region of the spectral line. We measure v 2 obs (ν i ) from the Doppler velocities in our observations, and estimate ρ(z), v gr and T from the self-consistent RADYN simulations.
We derive T by computing the square root of the ratio between the Doppler velocity power in the synthesized Ca II 854.2 nm line (Panel c) in Figure 20 and the power of the actual velocity field in the RADYN atmosphere at the formation height of the spectral diagnostic for each frequency bin (Panel b) of Figure 20. For the effective formation height of Ca II 854.2 nm, we chose the peak of the velocity response function. For our models, presented in Figure 20, the peak of the Ca II 854.2 Doppler velocity response function was at 1150 km, which is at about the same height as the average τ = 1 surface for the line core. The extent of the velocity response function is shaded as the red region in Figure 19. We use 5 mHz frequency averaging windows as the coefficient T and the power does not change significantly over this region and the averaging removes the inherent uncertainty in the power spectra. The inferred transmission coefficient T for Ca II 854.2 nm is presented in Panel d) of Figure 20. The Doppler signal of the low frequency waves (below 10 mHz) is less attenuated as their wavelength is significantly larger than the formation layer of the spectral diagnostic. However, the high frequency wave signal is significantly attenuated, as shown before (Bello González et al. 2009).
We note that our RADYN simulations can generate a high-frequency velocity signal due to the steepening of the acoustic waves propagating from the photosphere. The waves steepening into shocks in the chromosphere have saw-tooth-like velocity profile, which when Fourier decomposed creates power law tails that extend to significantly higher frequencies than that of the driving wave (Vecchio et al. 2009). The modeling approaches in the previous studies (for example Bello González et al. 2009) using monochromatic high-frequency waves and static 1D atmospheric models overestimate the transmission coefficient T (and the inferred wave flux respectively) at higher frequencies, as suggested by our estimation of T from the RADYN simulations. This is due to the fact in our dynamic solar atmosphere models the propagating wave packets are not monochromatic and their steepening creates high frequency signal.

Estimating the acoustic flux in the ALMA data based on RADYN
There are two differences between the temperature fluctuations observed with ALMA and the ones in the right hand side of Equation (6). First and foremost, the temperature response of ALMA will be a result of the wave fluctuations convolved with the atmospheric response function. Another complication, that makes Equation (6) applicable with limited validity in its current form is the fluctuating height of formation of the ALMA continuum (Molnar et al. 2019), which makes using a particular height of formation (and local plasma density at that height) nonphysical.
To take into account of all the effects described in the the previous paragraph, we rewrite Equation (6) as: where the proportionality coefficient, called from now on the ALMA transmission coefficient, C ALM A (ν) encapsulates the local plasma properties at the formation height of the ALMA radiation. We note that the C ALM A (ν) coefficient is the proportionality coefficient between the brightness temperature fluctuation power and the wave energy flux in the atmosphere, whereas the attenuation coefficient T in Section 4.3 is the attenuation of the observed wave velocity amplitude and goes in the denominator of Equation (7). By using the RADYN atmosphere models to synthesize synthetic ALMA observables and compare them with the observations we take into account those two effects as described below. We use the RH code (Uiten-broek 2001) to synthesize the millimeter continuum from the RADYN model output. We used the RADYN atmospheric models (including the instantaneous electron densities and hydrogen level populations) as an input for the RH code to synthesize the mm-wave radiation corresponding to Band 3 (100 GHz/3 mm) and Band 6 (240 GHz/1.25 mm). The RH code takes into account the opacity from neutral hydrogen as well as H − and H − 2 , which are the main sources of opacity in the solar atmosphere at millimeter wavelengths (Zlotnik 1968). We averaged the output intensity from the 0.5 second time steps of the RADYN runs to the 2 second cadence of our ALMA observations. The extent of heights at which the ALMA Bands sample the atmospheric plasma temperature are shaded in Figure 19 as the green region for Band 6 and the blue region for Band 3. We have determined those regions as the height of τ = 1 surface for the millimeter continuum from the RADYN run most closely reproducing the observations. For Band 6 this was model 19000 and for Band 3 is model 3000 as shown in Figure 20. The mean formation height for the ALMA Band 6 is at 700 km and for Band 3 is 1150 km, while the physical widths of the formation regions shown in Figure 19 correspond the variation of the τ = 1 height in our models. Those heights are lower than the ones previously presented in the literature based on modeling (Molnar et al. 2019;Martínez-Sykora et al. 2020) and observational (Patsourakos et al. 2020) methods.
The middle and right columns of Figure 20 present the results from the ALMA spectral synthesis with the RH code from the RADYN simulations and the resulting C ALM A (ν) coefficient. The left column corresponds to ALMA Band 3 (3 mm continuum) and the right column corresponds to ALMA Band 6 (1.25 mm continuum). Panels e) and i) of Figure 20 show the temporal variation of the synthesized brightness temperature in the RADYN runs, which agree well with previous studies of the millimeter continuum based on RADYN simulations (Loukitcheva et al. 2004;Eklund et al. 2020). However, we do note that the average synthetic temperatures of the RADYN models (4250/5250 K for Band 3 and 6 respectively) are significantly lower than the observed ones (7000/8500 K for Band 6 and 3) in Figure 4. Panels f) and j) show the acoustic flux PSD at the average height of formation of the millimeter radiation (700/1150 km for Band 6/3 respectively), with a clear correlation between the amount of acoustic flux and the amplitude of the brightness temperature fluctuations. This correlation is further demonstrated in panels g) and k) of Figure 20, where we present the PSD of the modeled brightness temperature fluctuations. This correlation is not surprising, as compressive waves have temperature perturbations and ALMA measures the plasma temperature in the chromosphere.
C ALM A is presented panels h) and l) of Figure 20 where we have calculated it as the ratio of the acoustic flux density at the formation height of the millimeter radiation (panels f) and j)) and the power density of the relative temperature fluctuations (panels g) and k)). We have averaged the transmission coefficient C ALM A over 5 mHz frequency windows to smooth out the inherent noise in the power spectra. The C ALM A coefficient converges to the same values for the RADYN runs with parameters closest to our observations (models 3000-19000). The converging values for C ALM A at different acoustic fluxes confirm that our modelling approach is not strongly dependent on the wave amplitudes (Equations (8)) and evades the complications of using Equation (6) directly for estimating the wave flux.
Due to the spatial smearing from the finite PSF of the ALMA beam the observed oscillatory power is underestimated by a factor of 2 based on previous studies for Band 3 (Loukitcheva et al. 2006(Loukitcheva et al. , 2015Wedemeyer et al. 2020). We do include this factor in our analysis only for Band 3, as we don't have an estimate for it for Band 6. From the RADYN synthetic observables presented in the previous section, we had all the constituents of Equation (7) to compute the acoustic flux in our observations. The mean density at the mean formation height of the Ca II 854.2 nm line is 5×10 −9 kg m −3 (at 1150 km height in the model atmosphere). This is comparable to estimates from previous work (Abbasvand et al. 2020b). The values for the attenuation coefficient T are presented in the bottom left panel in Figure 20 for the RA-DYN models with different driver strengths. We used the values for T taken from the model producing the closest synthetic velocity/brightness temperature variations to the real data for the particular frequency bin and diagnostic. This was the model 3000 run for our IBIS observations. This was also the RADYN run with the highest piston amplitude that didn't cause the core of the Ca II 854.2 nm line to flatten out and go into emission (a condition not widely seen in our data). We see that the attenuation coefficients T from different model runs converge to similar values at frequencies under 30 mHz. Above that, the attenuation coefficient T for different models do not agree that well, but this frequency region contributes relatively little to the total acoustic flux. Therefore, the exact driver strength does not appear to be important in the determination of the attenuation coefficient and the derivation of the acoustic flux.
We computed the wave energy flux PSD in our observations based on our estimates for density and T from the RADYN simulations and the results are presented as the magenta data in the second row in Figure 20. We can see that the energy flux PSD derived from the observations is significantly higher in Ca II observations compared to the simulations. However, the shape of the wave flux frequency distribution is almost flat up to 20 mHz and resembles the one found in the RADYN models.
At each pixel, we divide the observed Ca II 854.2 nm line core velocity power in each frequency bin by the corresponding value of T for model 3000 and sum over the interval from 5 to 50 mHz to derive the total acoustic flux, with the results presented in Figure 21 for the FOV 2 observations from 17:04 -17:11 UT. We verified this result by calculating the same flux estimate using another FOV 2 dataset obtained from 15:39 to 15:46 UT, with very similar results. This suggests that the derived acoustic fluxes are not significantly affected by the seeing variations or evolution of the granular or supergranular conditions in the photosphere.
According to Withbroe & Noyes (1977), the average radiative losses in the middle (high) chromosphere are ranging from about 2 kW m −2 (0.3 kW m −2 ) for the quiet internetwork to 20 kW m −2 (2 kW m −2 ) for the plage regions. The distribution of the acoustic flux in Figure 21 suggests that over most of the FOV the acoustic flux can not be the dominant source of heating in the chromosphere. There are some disjoint regions with fluxes above 1 kW m −2 (shown in dark blue) where acoustic waves could be a significant source of chromospheric heating. These regions of enhanced wave flux are located primarily in the more fibril-free internetwork areas and also in the network and plage. The locations dominated by chromospheric fibrils (see Figure 2) almost uniformly have acoustic fluxes less than 0.1 kW m −2 .
The distributions of observed acoustic flux in the different solar regions (using the mask in Figure 5) are presented in Figure 22. The averages and the 10th/90th percentiles of the cumulative distribution are summarized in Table 2. The regions with the highest inferred acoustic flux are the plage and internetwork regions. However, after further investigation of the shape of the spectral lines in the plage regions with highest acoustic fluxes we found that the spectral line core fills in (flattens), and the Doppler velocity derived from the line core minimum becomes less reliable. The line core T coefficient for the different RADYN models. Middle column: results for Band 3 (3 mm); Right column results for Band 6 (1.25 mm). The panels in these columns share the same information for the two ALMA Bands: Panels e) and i): synthesized ALMA brightness temperature from the RADYN models; Panels f ) and j): Similar to Panel b), but for the formation height of ALMA Bands 3 and 6; Panels g) and k): Similar to Panel c, but for the observed brightness temperature fluctuations in ALMA Bands 3 and 6; Panels h) and l): The CALMA coefficient for the different RADYN models. The color coding throughout this figure is consistent and corresponds to the same models presented in the legend of the bottom row panels.  of Ca II 854.2 nm is determined by the local temperature at the formation height of the line ). Strong impulsive heating in the plage regions can be responsible for this dynamic peculiarity of the Ca II 854.2 nm line profile. Therefore, we believe that the high velocity oscillation power observed at hotter (plage) locations might not be due to genuine wave motions but because of the the flattening of the line core. Furthermore, the strong magnetic field in the plage regions might lead to MHD wave effects beyond the scope of this paper. We observe a similar phenomenon in our RADYN simulations, where if the wave driver is injecting too high velocity perturbations (acoustic flux) we see the line shape of Ca II 854.2 nm going into emission, when the shocks pass through the chromosphere. However, on average the 1D simulations produce too narrow line core profiles (as previously shown in Leenaarts et al. 2009) which could be due to insufficient microturbulence in our simulation (RADYN has a default microturbulence of 2 km/s). Hence, we leave the detailed investigation of the dynamics of plage spectral line behavior for a future study.
Further analysis of the frequency dependence of the wave flux show that about 60% of the wave energy flux is in the 5-20 mHz frequency range, 30% in the 20-40 mHz range and 10% in the 40-60 mHz range. The pixels with significant relative contribution from the 20-40 mHz frequency range have lower total acoustic flux and are mostly found in the plage and the network.
In summary, our observations show that acoustic wave dissipation is likely not the dominant heating mechanism for the middle chromosphere. However, acoustic waves could contribute significantly to the quiescent state of the upper chromosphere in the internetwork as the observed flux is of the order of magnitude of the radiative losses in that layer. Future work will constrain the wave flux and dissipation in this layer using IRIS observations.

Inferred wave flux from the ALMA observations
Using the RADYN results presented in Section 4.4 we are able to estimate the acoustic fluxes in the ALMA observations. We match the observed brightness temperature fluctuations (the magenta points in panels g) and k) in Figure 20) to the closest RADYN model by the amount of brightness temperature RMS and then use the C ALM A coefficient from that RADYN run (bottom row in Figure 20) to calculate the amount of acoustic flux. To calculate the acoustic flux for Band 3 we used the observing block obtained between 17:31 and 17:41 UT and the C ALM A values for model 3000. For Band 6 we used the observing block between 16:03 and 16:11 UT and the C ALM A values for model 19000.
The inferred acoustic flux for Bands 3 and 6 are presented in Figure 23. Circular masks were applied around the edges field of view as the noise outside of those regions are significant due to the decreasing sensitivity. We observe higher acoustic flux in Band 6 compared to Band 3, which is expected if the wave flux is being dissipated as the waves propagate upward. The region of the FOV of Band 6 is shown on the Band 3 FOV as the green circle. Since the two different ALMA Bands were observed an hour apart, they do agree to a certain extent but not fully as the solar atmosphere is changing on shorter time scales. Furthermore, the very limited FOV of Band 6 makes comparisons difficult. Examining the frequency distribution of the wave flux shows that the dominant source of signal in both ALMA Bands are the frequency range between 5 and 20 mHz, where more than 70% of the signal is found. This agrees with our results for the Ca II 854.2 nm line, that the frequency range between 5 to 20 mHz contains most of the wave flux.
To compare the results from IBIS and ALMA, we compare the derived acoustic flux over a common FOV corresponding to that of the Band 6 data. The distributions of the observed fluxes at those locations are presented in Figure 24. The ALMA Band 6 data exhibits the highest flux. They are followed by the Ca II 854.2 nm and the ALMA Band 3 data. This ordering of the amount of wave flux follows the height of formation of the diagnostics shown in Figure 19. Since Band 3 and 6 have similar systematics in their formation height we can estimate the dissipation between their formation heights. The average dissipated energy (flux difference) across the FOV of Band 6 data is around 0.7 kW m −2 . However, there is a significant high-power tail in the dissipated energy which is greater than 1 kW m −2 and is energetically significant to maintain the quiet chromosphere in some confined regions. However, on average this is not enough to sustain the quiet middle chromosphere, being too small by an order of magnitude.

CONCLUSIONS AND FUTURE WORK
We obtained an extensive data set containing spectral observations covering from the upper photosphere and the middle (Ca II 854.2 nm, Hα, and ALMA Band 6) and upper chromosphere (ALMA Band 3). The rapid cadences of our observations allow us to study high frequency dynamics of the chromosphere in a rarely studied frequency regime. Our observations have extended the work of  to show that the powerlaw distribution of the Fourier PSDs is ubiquitous in all of the observed velocity and temperature diagnostics. In particular, in this paper we have focused on the velocity diagnostics derived from Ca II 854.2 nm, as tracers of upward propagating compressive waves, as well as brightness temperature fluctuations from ALMA Band 3 (3 mm) and 6 (1.2 mm) as indicators of local heating from those waves. We found that the power law properties depend on the observed solar region. We confirmed that the white-noise level in the power spectra is consistent with the photon shot noise in our data (see Appendix A), indicating that seeing-induced crosstalk or other sytematic noise is likely not responsible for the detected Doppler velocity power. The amplitude of the noise floor increases in more "active" regions of the field (e.g. plage), perhaps due to flatter line cores typical of these features. Furthermore, the slope of the spectral diagnostics' power laws changes with the observed regions. "Hotter" regions (like network and plage) exhibit less steep power law slopes as shown in Figure 11. This flattening of the power law could be due to the filling of the core of the spectral profile of the Ca II 854.2 nm line in the hotter regions (network and plage). We believe the steeper power law is an important characteristic of the dynamic nature of the different solar surface features that has to be reproduced in future modeling efforts.
We compared the RMS values of velocity oscillations observed at different incidence angles (µ = 0.98 vs µ = 0.41). The RMS of the velocity fluctuations is significantly smaller at the higher incidence angle compared to the disc center observations, but still larger than what would be expected if the amplitudes of those fluctuations were only due to the inclined viewing of vertically propagating waves of the same magnitude as seen in our disk center observations. Modeling the observed signal as a superposition of a longitudinal (vertical) component and a transverse (possibly Alfvénic) component, we found that the transverse component had a RMS amplitude of about 0.15 km/s, compared to the vertical component amplitude of around 0.5 km/s.
To characterize the acoustic waves that could explain the power laws in our observations we used the RA-DYN code to model propagation of waves from the upper convection zone into the chromosphere. We used wave drivers (as bottom boundary conditions) similar to the ones in Fossum & Carlsson (2006) adjusted with a scaling factor. We ran the RHD models and then produced synthetic observables for both IBIS and ALMA, using the RADYN built-in radiative transfer module in the former case and the RH code in the latter case. The dynamic RADYN models are able to reproduce the features of our observations in terms of total oscillatory power and slopes. Hence, we were able to correlate our observed oscillatory power in different diagnostics to the actual acoustic flux present at different heights in the simulations.
The acoustic flux derived from the Ca II 854.2 nm line Doppler velocity data is estimated to be between 0.1 to 1 kW m −2 . The lowest amount of flux is found in the penumbra and fibril regions and the highest in the internetwork and plage regions. We believe that the inclined nature of the magnetic field in the penumbra and the fibril region plays role in the observed lower fluxes.
Furthermore, as discussed in Section 5.1, the high values in the plage region above 1 kW/m 2 require further  examination due to the changes in the Ca II 854.2 nm line profile in plage regions, which leads to spuriously high measured Doppler shifts.
Furthermore, we found that most of the contribution to the acoustic flux comes from the 5-20 mHz frequency interval (around 60% from the total acoustic flux) and that high frequency waves above 40 mHz do not contribute significantly (less than 10% of the total observed flux).
We also compared the brightness temperature fluctuations in our ALMA data with the synthetic observables from the RADYN runs and inferred the acoustic flux by using the correlation between brightness temperature fluctuations and acoustic flux in the different RADYN model runs. Based on this comparison we can infer that Band 6 has acoustic flux on average of 0.7 kW m −2 , compared to the formed higher in the atmosphere Band 3 which contains about 0.03 kW m −2 . From these two observations at different heights, we can infer that the average wave flux dissipated between two layers probed by ALMA is about 0.7 kW m −2 , which is not sufficient to heat the middle chromosphere, but is a significant contribution to its energy budget. This result agrees quantitatively with previous work by Nindos et al. (2020), who used ALMA to compute the heating in small scale chromospheric brightnenings. However, in certain regions the dissipated wave flux exceed the threshold of 2 kW m −2 which is sufficient to maintain the quiet chromosphere locally. We believe that the limited spatial resolution of the ALMA observations could lead to an underestimation of the wave flux, compared to our Ca II data (Loukitcheva et al. 2015;Wedemeyer et al. 2020) and further observations with higher angular resolution (more sparse ALMA configuration) will provide better constraints on the wave flux.
Another peculiarity between our observations and ALMA data is the fact that the Band 3 corresponds most closely to model 3000 while Band 6 to model 19000. Model 19000 has twenty times higher wave energy flux than model 3000 in the chromosphere. This discrepancy could be due to different reasons -either the heights of formation of the ALMA continuum are inaccurate in the RADYN models or there is a wave dissipation mechanisms not included in the RADYN models. If the height in the atmosphere where the ALMA continuum originates from is determined inaccurately, this will lead to incorrect acoustic flux determination and result in differing RADYN models corresponding to the observations. The latter possibility of missing physics is also very probable due the 1D hydrodynamic nature of the RA-DYN code which might be omitting the required physics to treat fully all the relevant wave damping mechanisms.
Our work raises further questions such as what is the role of the magnetic field on the wave propagation characteristics in the solar atmosphere? Future observations with multiwavelength spectropolarimetric capabilities throughout the photosphere and the chromosphere from Daniel K. Inouye Solar Telescope (DKIST) (Rimmele et al. 2020;Rast et al. 2021) will be able to address that. Also, the higher throughput of the future generation solar telescopes will help with driving down the white noise floor and provide simultaneous spectral observations at different heights in the solar atmosphere. Another interesting aspect we will pursue in a following publication is the amount of velocity oscillations in the upper chromosphere and the transition region observed cotemporally with IRIS during our April 2017 campaign.
The numerical side of our work requires further investigation as well. One important question that will be addressed in future work is the sensitivity of our results on the numerical setup that was utilized. For example, how does the transmission coefficient T depend on the model atmosphere and the number of grid points in the atmosphere? Furthermore, studying wave propagation in 3D is essential for understanding the observed signals as the nature of the observed chromospheric structures is strongly non-vertical and non-local (Carlsson et al. 2019;Eklund et al. 2021). Even though current 3D RMHD solar models have significantly differing wave propagation properties , further modeling with realistic solar atmospheres in three dimensions is essential.

A. ESTIMATING THE NOISE FLOOR IN THE POWER SPECTRA OF THE OBSERVED DOPPLER VELOCITIES
The power spectra presented throughout this paper (for example in Figures 7 and 9) exhibit white noise behavior at the high frequency limit. To interpret the derived properties from the observed power spectra correctly we need to understand the noise sources in our data that could contribute to the wave signal. Spectral line profiles with different shape would be affected to a different degree from photon noise due to their varying shapes and intensity levels. For example, deeper and narrower profiles chromospheric profiles of the Ca II IR lines would be less susceptible to Doppler velocity measurement shot noise, compared to plage or AR shallower spectral profiles whose cores fill up and flatten. Hence, we modeled the effect of photon shot noise on each different solar surface feature described in Figure 5.
To estimate the white noise floor properties due to photon shot noise to the Ca II 854.2 nm line velocity power spectra we adopted a Monte-Carlo approach. We chose 150 random spectral profiles from the each region in Figure 5. We did not choose the average spectral line profile for each region of the Sun to be representative as averaging over space and time does not represent the instantaneous realization of the spectral profiles. We computed 1500 noise realizations for each chosen spectral line profile. The noise for each wavelength point was calculated by using the ADU (2.5 e-/DN) of the camera (Andor iXon 885) used that day to calculate the number of photons. Since the signal to noise ratio was on the order of a few hundred (even for the line cores), we applied a Gaussian noise statistics to the estimated photon shot noise levels. The Doppler velocity from the simulated time series of 1500 noise realizations was measured with the same techniques used for reducing our IBIS data (described in Section 3.2). The measured Doppler velocity power spectrum density was white noise, as expected for uncorrelated noise. The bottom panel of Figure 8 presents the distribution of the median noise level in our estimation. The synthetic noise distributions match well with the observed ones. Hence, we can assume that seeing induced crosstalk is not the dominant source for the high frequency white noise floor.