Spatio-temporal analysis of chromospheric heating in a plage region

Our knowledge of the heating mechanisms that are at work in the chromosphere of plage regions remains highly unconstrained from observational studies. The purpose of our study is to estimate the chromospheric heating terms from a plage dataset, characterize their spatio-temporal distribution and set constraints on the heating processes that are at work. We make use of NLTE inversions to infer a model of the photosphere and chromosphere of a plage dataset acquired with the Swedish 1-m Solar Telescope. We use this model atmosphere to calculate the chromospheric radiative losses from H i, Ca ii and Mg ii atoms. We approximate the chromospheric heating terms by the net radiative losses predicted by the inverted model. In order to make the analysis of time-series over a large field-of-view computationally tractable, we make use of a neural network. In the lower chromosphere, the contribution from the Ca ii lines is dominant and located in the surroundings of the photospheric footpoints. In the upper chromosphere, the H i contribution is dominant. Radiative losses in the upper chromosphere form an homogeneous patch that covers the plage region. The net radiative losses can be split in a periodic component with an average amplitude of ampQ = 7.6 kW m^{-2} and a static (or very slowly evolving) component with a mean value of -26.1 kW m^{-2}. Our interpretation is that in the lower chromosphere, the radiative losses are tracing the sharp lower edge of the hot magnetic canopy, where the electric current is expected to be large. In the upper chromosphere, both the magnetic field and the distribution of net radiative losses are room-filling, whereas the amplitude of the periodic component is largest. Our results suggest that acoustic wave heating may be responsible for one third of the energy deposition in the upper chromosphere, whereas other heating mechanisms must be responsible for the rest.


Introduction
The heating of the solar chromosphere and corona remains one of the foremost questions in solar and stellar physics. The chromosphere is on average radiating 4 kW m −2 in the quiet Sun and 20 kW m −2 in active regions (Vernazza et al. 1981;Withbroe & Noyes 1977). At least, that energy must be transported and deposited into the chromosphere at any time by heating mechanisms. Although we cannot measure chromospheric heating terms directly, we can estimate them by assuming that they equal the radiative losses in the main chromospheric coolers, typically strong chromospheric lines and continua from the H i, Ca ii and Mg ii atoms.
The physics and heating of plage regions have puzzled our community since the 70's. Three recent independent studies have attempted to infer the strength and stratification of the magnetic field in plage targets (Morosin et al. 2020;Pietrow et al. 2020;Ishikawa et al. 2021). They found amplitudes of approximately |B |∼ 300 − 400 G, depending on the spectral line and target under analysis. In particular Morosin et al. (2020) reconstructed the canopy effect of the magnetic field in the chromosphere. The magnetic field is very concentrated in the intergranular lanes in the photosphere and it expands horizontally as we move up in the atmosphere, forming a hot magnetic canopy over the photosphere. Because of the sharp lower boundary of the canopy, the authors speculated that current sheets should be present in this boundary, purely from the application of j = ∇ × B/µ. Those currents could lead to Ohmic dissipation at the lower boundary of plage, causing heating at the base of the chromosphere.
Furthermore, modelling chromospheric lines from plage observations have typically required very large values of microturbulence, typically up to 10 km s −1 (e.g., Shine & Linsky 1974;Carlsson et al. 2015;De Pontieu et al. 2015;Carlsson et al. 2019). A more recent study using ALMA and IRIS observa-Article number, page 1 of 15 arXiv:2203.01688v2 [astro-ph.SR] 30 Mar 2022 A&A proofs: manuscript no. main tions have further constrained this value to an average value of 5 km s −1 (da Silva Santos et al. 2020). Whether these large values of microturbulence are related to turbulent velocity fields, sharp gradients along the line-of-sight induced by a hot magnetic canopy above the photosphere (Sanchez Almeida & Martinez Pillet 1994;Buehler et al. 2015;Morosin et al. 2020), answering this question is entangled with the enhanced values of radiative losses that have been reported in plage in comparison with the quiet Sun.
In the chromosphere magnetic forces equal those from pressure gradients, leading to a complex and very dynamic force balance. Therefore magnetoacoustic waves, turbulent Alvén wave dissipation, magnetic reconnection, Ohmic current dissipation, viscous heating and ambipolar diffusion can all contribute to the heating of the chromosphere (see, Hasan & van Ballegooijen 2008;van Ballegooijen et al. 2011;Khomenko & Collados 2012;Priest 2014;Martínez-Sykora et al. 2017;Priest et al. 2018;Brandenburg & Rempel 2019;Yadav et al. 2020;Díaz Baso et al. 2021b;da Silva Santos et al. 2022, and references therein). The exact set of processes that are at work in plage and their contribution to the energy budget remains highly unconstrained from observations. A great discussion and exhaustive review about the potential contribution of different heating mechanisms in plage is presented by Anan et al. (2021). They also analyzed the correlations of the radiative flux in the Mg ii h&k lines and with the magnetic field strength that was inferred from observations in the He i 10830 Å line. Their conclusion was that Alfvén waves or ion-neutral collisions could be heating plage regions. They could not find a clear correlation with the electric current in their results.
Inversion methods allow for the inference of a model atmosphere of the photosphere and chromosphere by iteratively modifying the physical parameters of a model atmosphere to reproduce the observed full-Stokes spectra. They can be used for a variety of spectral lines and it is also possible to include non local thermodynamical equilibrium (NLTE) effects (Asensio Ramos et al. 2008;Socas-Navarro et al. 2015;Milić & van Noort 2018;de la Cruz Rodríguez et al. 2019;Ruiz Cobo et al. 2022). NLTE inversions are computationally expensive and time consuming, so studying the evolution of the atmosphere parameters of an entire time series observation could become prohibitive. The great advantage of inversions is that we can use the inferred model atmosphere to calculate radiative losses in the chromosphere (see e.g., Abbasvand et al. 2020;Díaz Baso et al. 2021b) and thereby, obtain a lower-limit estimate of the chromospheric heating terms. To our knowledge, there is no other way to estimate radiative losses from observational data.
In this study, we have made use of a subset of a long time series to calculate plage models from NLTE inversions. We have used the resulting model atmospheres to train a neural network that can quickly predict the model atmosphere for the rest of the dataset, in a similar way to Asensio Ramos & Díaz Baso (2019) or Kianfar et al. (2020). The underlying assumption for this approach to work is that the training set is statistically representative of the entire time series. The quantities included in a model atmosphere are the gas temperature T , the line-of-sight velocity v LOS , the turbulence velocity v turb , and the parallel and perpendicular components of the magnetic field, respectively B || and B ⊥ .
We have calculated the net chromospheric radiative losses for all time steps of the series in order to better understand the distribution of radiative losses in a plage target, as well as the time evolution. In our analysis, we study some correlations with other physical parameters in order to suggest what heating mechanisms could be at work.

Observations and data reduction
The target of interest is the NOAA2591 and it was observed with the Swedish 1-m Solar Telescope (SST; Scharmer et al. 2003) on the 14th September 2016 at 08:26 UT. It is a plage region located at (X,Y) = (424 ,-16 ), that corresponds to a viewing angle of µ = 0.90. The CRisp Imaging Spectro-Polarimeter (CRISP; Scharmer et al. 2008) and the CHROMospheric Imaging Spectrometer (CHROMIS; Scharmer 2017) were used in order to obtain observations in Ca ii 8542 Å(full-Stokes), Fe i 6302 Å (full-Stokes) and Ca ii H&K (intensity only).
Ca ii 8542 Å was sampled at 21 wavelength positions. The distance between each position was ∆λ = 85 mÅ, except for two points in the far wings of the line located at ±1700 mÅ from line center. Fe i 6302 Å was recorded at 16 wavelength positions, of which, one at line center and two of them located in the red part part of the line at ∆λ = 40, 80 mÅ. The other 13 points are located at ∆λ = −40 mÅ from line center. The data detected with CRisp have a cadence of ∆t = 37s. In total the observation with the instrument lasts for t ∼ 22min.
The Ca ii H&K lines were sampled at 21 wavelength positions, of which 19 at a distance of ∆λ =79 mÅ from each other, and the two extra-wing points at a distance of ∆λ =1.25 Åfrom line center. A continuum position at 4000 Åwas also recorded. CHROMIS has a cadence of ∆t = 16s. In order to match the observations with the two instruments, for each snapshot of CRisp, we chose the closest snapshot of CHROMIS in time.
After the acquisition, the data have been reduced using the SSTRED pipeline Löfdahl et al. 2021). In order to take into account atmospheric effects, the data have been also processed using the Multi-Object Multi-Frame Blind Deconvolution method (MOMFBD) described in van Noort et al. (2005). Then the dataset has been properly aligned since the two instruments have a different pixel scale. We also have used the python package ISPy to handle the data and metadata of the data cubes (Díaz Baso et al. 2021).
An overview of the observed active region is presented in Fig. 1. The panels depict the observations in Ca ii 8542 Å, Fe i 6302 Å and Ca ii K. Stokes I and V/I of Ca ii 8542 Å are shown at ∆λ =85 mÅ from line center. The upper left panel of Fig. 1 illustrates the plage region in the chromosphere, with typical features as fibrils, extending from the center of the region towards the outside. There are two different polarity patches in the field-of-view (FOV). The Ca ii K wing image in the lowerright panel shows low-lying bright structures connecting both polarities.

Inversions
To estimate the thermodynamic and magnetic properties of the region we have performed non-LTE inversions, using the STockholm inversion Code (STiC) (de la Cruz Rodríguez et al. 2016. It is a modified version of the radiative-transfer RH code (Uitenbroek 2001) and it includes a fast approximation to calculate the effects of partial redistribution (PRD, see Leenaarts et al. 2012). The inversion engine of STiC includes an equation of state extracted from the Spectroscopy Made Easy (SME, Piskunov & Valenti 2017). The radiative transport equation is The spectra in the Ca ii spectral lines were calculated in NLTE by assuming statistical equilibrium and plane-parallel geometry. Furthermore, PRD effects were explicitly included in the calculations of the H&K lines. The Fe i lines were synthesized assuming local thermodynamical equilibrium. We chose to perform the inversions in a column-mass scale. In comparison to optical-depth, column-mass allows computing the gas pressure scale directly, without involving the equation of state or background opacities. We will see later that this feature is important when using the neural network to predict models from observations.
We inverted a training set for a neural network. The training set consisted of the full field-of-view for two non-consecutive time steps of the series, but in order to speed up the calculations we only inverted every seventh pixel of the FOV. Furthermore, in order to efficiently train the neural network it is more convenient to have a statistical picture of FOV, with many pixels spread across the region and many different observed spectra, rather than many pixels in a small area that could all present similar observed spectra.
We initialized the magnetic field vector of the initial input model using the spatially-regularized weak-field approximation proposed in Morosin et al. (2020). The presence of strong velocity gradients as a function of depth can greatly distort the line profile of very strong lines. In out tests, this was one of the main sources of degeneracy in the output models. In order to better sample the parameter space in line-of-sight velocity, we prescribed five initializations of the stratification that were added to the FAL-C model (Fontenla et al. 1993): three with constant values at 0, ±5 km s −1 , and two with strong upflowing and downflowing gradients (see Fig. 2). The number of nodes used to run the inversion are the same for all five models. We used 10 nodes in temperature, 4 in v LOS , 4 in v turb , 3 in B || , 2 in B ⊥ and 1 node in φ. For each pixel, we selected the model that yielded the best χ 2 value and then applied a mild horizontal smoothing. The inversions were re-started with an increased number of 11 nodes in temperature and 4 nodes in B . Once all cycles were finished, we switched on the NLTE equation of state and re-run the final cycle, that essentially accounts for the ionization of hydrogen in NLTE and the ionization of the rest of elements is calculated in LTE. The last cycle essentially affects the temperature stratification and the derived electron densities. The first two rows in Fig. 3 present respectively the T and v LOS for three different depths in the atmosphere, that correspond to lower photosphere, upper photosphere and upper chromosphere. For B || , B ⊥ and v turb just two depths have been chosen: they are represented in each column of the two bottom rows.
The inversion results reproduce many features of the solar atmosphere that are well known from past research. In the photospheric temperature map, for example, the granulation pattern is visible, while moving upwards in the atmosphere the hotter and elongated fibrilar structures extend toward the outside. Moreover, the magnetic field presents the typical plage structure with strong concentration of field and field-free gaps in the photosphere, while in the chromosphere it is more extended and less strong (Buehler et al. 2015;Morosin et al. 2020). We note that the Q and U signals in the chromosphere are very noisy and the inversion code struggles to reconstruct clean maps for the |B ⊥ | and the azimuthal components of the magnetic field.
Due to the S/N ratio of the observation and to the different sensitivity of the emerging intensities to the different parameters of the model, the depth resolution of the model is greatest in temperature, and line-of-sight velocity, and much more limited in microturbulence and B . Therefore the magnetic canopy does not appear as sharply in the reconstructed magnetic field stratification as in the temperature reconstruction and in both cases it is much smoother than it probably is in reality because of the limitations of a nodes-based inversion.

Neural network application
Since it would be extremely time consuming and very computationally expensive to invert all the 35 time frames of the observation with STiC, we suggest instead a easier and faster approach using neural networks. They have showed a good performance in terms of accuracy and speed for solving different problems, for example to identify and predict solar flares (Panos et al. 2018), to denoise solar observational images (Díaz , or to study the mapping between spectral lines and the solar atmosphere (Socas-Navarro 2005; Centeno et al. 2021). In this case, we have used the resulting model atmospheres from the inver-sion to train a neural network that can quickly predict the model atmosphere for the rest of the dataset (Asensio Ramos & Díaz Baso 2019). We refer the reader to Appendix A for a detailed explanation of the architecture, training process and validation of the neural network.

Results of the NN:
The results obtained with the training of the NN are shown in Fig. 4 in a format similar to Fig. 3. The represented time step is the same as in Fig. 3.
The temperature and line-of-sight velocity prediction from the NN are very well correlated with the results from the inversions (see 2D density plots in Fig. 5). Both components of the magnetic field vector are predicted smoother compared to those obtained with the inversions. The latter occurs because the noise from Stokes Q and U was dominating the inversion results and STiC was not able to fully reconstruct the signal. The introduction of a constant factor in the training helps the neural network to smooth the noise and to better estimate the values.
The NN seems to overestimate v turb very deep in the photosphere, while in the chromosphere this behaviour disappears. The neural network incorrectly correlates areas of strong photospheric B || with higher values of the microturbulence. We investigated ways of removing this degeneracy, but we could not find a solution. We decided to move forward regardless because the values of the microturbulence close to the continuum formation layer in the photosphere have no influence in the prediction of radiative losses through strong chromospheric lines.

Radiative cooling rates
The heating terms in the chromosphere can be approximately estimated by calculating the integrated radiative losses, since the energy required to sustain radiative losses must be sustained by chromospheric heating terms. To calculate the radiative cooling rates from the predicted model atmosphere from the NN, we first imposed hydrostatic equilibrium in order to derive a z-scale and the gas pressure scale. Most inversion codes operate in an optical-depth scale and therefore the gas pressure has to be calculated iteratively with the consequent re-calculation of the continuum opacity (see, e.g., Mihalas 1970). However, working in a column mass scale, simplifies the calculations and no iterations are needed (Hubeny & Mihalas 2014): where ξ is the column mass, known from the inversions, g is the solar gravity and p gas is the gas pressure.
We have calculated the electron densities in each depth point using a simple equation of state 1 proposed in Mihalas (1970) which includes hydrogen atoms bound in H − and H2 molecules. Once the electron densities were known, we estimated the total number of atoms for a given temperature assuming an ideal gas: where N a is the atoms number density, N e is the electron density, K B is the Boltzmann constant and T the temperature. In this way, it is possible to obtain the total number of atoms N a and, by multiplying it to the mean particle mass, we can obtain the density ρ at a certain depth in the atmosphere: where the mean particle mass < m > is given by: In this case, a i is the solar abundance of the i-th element and m i is its atomic mass.
Finally, we obtained the z-scale from the definition of column mass: therefore, in a discrete grid we can write: This hydrostatic z-scale is used to calculate the integrated radiative losses in the chromosphere. The latter are calculated integrating over a height interval, selected from the z-scale. Our final model contains all the quantities calculated under the assumption of hydrostatic equilibrium and the ones previously obtained after the inversion process.
We have used a modified version of STiC that evaluates directly Eq. (7) and also outputs the net radiative rates for all transitions and the atom population densities. In our case we considered Lyα, Ly continuum, Hα, the Ca ii H&K and the IR triplet lines and the Mg ii h&k and UV triplet lines. Our inversion setup includes lines that sample the solar atmosphere from the upper chromosphere to the photosphere. Although the Ca ii K line may be sensitive to the very lower part of the transition region, it is not sufficient to properly constrain its exact location and gradient. Therefore, the reconstructed transition region is very cold and extended in comparison to models from inversions including transition region diagnostics (see e.g., de la Cruz Rodríguez et al. 2016). The latter seem to have a large impact in the prediction of the Balmer continuum, which is unrealistically large in our calculations. We could not find an obvious solution to this problem, and therefore we did not include the Balmer continuum contribution in our study. The effect of this exclusion is that our radiative losses will potentially be even lower than the heating terms that we are approximating with them.
The inversions were calculated including the effect of NLTE hydrogen ionization in statistical equilibrium by imposing charge conservation (Leenaarts et al. 2007). But our NN does not predict the electron density. Therefore, we have first calculated a forward synthesis with the H atom, recovering not only the radiative losses but also the electron densities in NLTE. Then the electron densities in the original model have been replaced with the new ones in the input model. The synthesis was done for the Ca ii and Mg ii atoms with the updated electron densities.
The radiative cooling rates are calculated automatically inside the code by computing the divergence of the radiative flux (see e.g., Uitenbroek 2002;Rutten 2003): where α ν is the total absorption coefficient, S ν is the total source function and J ν is the mean intensity over solid angle. For boundbound transitions, Eq. (7) can be transformed into an expression that only depends on the net radiative rates and the level population densities: where h is the Planck's constant, ν 0 is the central frequency of the transition, n u/l are the population of the upper and lower level respectively and R ul/lu are the radiative rate coefficient from the upper to the lower level or vice versa. In order to obtain the integrated radiative losses, Q has to be integrated over the range of geometrical heights of the region of interest. In our case the integration limits are set in order to take into account only the chromosphere, so from the depth point after the temperature minimum to the depth point at which the temperature reaches T ∼ 10000K. We note that moving the integration limits can change the derived losses, and therefore some deviations can be expected when comparing numbers from different studies.

Single timestep
In order to study the spatial distribution of radiative losses, we have calculated the radiative losses for the initial time step of the observations. We have chosen this frame because it is one of those with better seeing conditions. We have divided our integration interval (in height) into three subregions (lower, middle, upper chromosphere) to understand how the energy deposition is taking place, which are shown in Fig. 6. The lower extreme of the interval has been chosen after the temperature minimum, that define the end of the photosphere and the beginning of the chromosphere. The upper extreme has been chosen in order not to include the transition region in the calculations. Since our inversions did not include lines that are strongly sensitive to the transition region, the steep temperature gradient and its exact location are not well constrained in our models. The derived radiative losses integrated over the entire interval and over the three subregions of the chromosphere, are shown in the top row of Fig. 7. The average integrated radiative losses over the whole plage region are ∼ −28 kW m −2 . The FOV shows very small-scale regions with peak values of ∼ −90 kW m −2 . In order to have a more precise insight about the calculated radiative losses, we have plotted in Fig. 8 the contributions of each atom in the tree sub layers of the chromosphere. The Ca ii contribution is dominating in the lower and middle chromosphere while the contributions from H i and Mg ii lines are negligible in this region. The middle chromosphere is the sub layer that presents the lower values associated with the radiative losses. In the upper chromosphere, Hydrogen is the main contributor to the radiative losses. The contribution from the Mg ii atom is approximately one half of that from Hydrogen.
The comparison of the integrated radiative losses with the line-of-sight component of the magnetic field provides great in-  Bottom row: Maps of the parallel magnetic field for four different heights in the solar atmosphere. The height is increasing from left to right. The first panels in the middle and bottom rows represent respectively the T and B || in the photosphere. The green contours indicate the area where Q < const * Q layer kW m −2 in the corresponding atmosphere layer, whereQ layer is the median value of Q in the corresponding layer. Const = 1.6 for the lower and middle chromosphere, const = 0.7 for the upper layer.
sight into the overall heating process (see Fig. 7). From left to right, the contours plotted in three of the four panels correspond to a fraction of the median radiative losses value (Q < const * Q layer kW m −2 , whereQ layer is the average net radiative loss in that layer) in the lower chromosphere, in the middle chromosphere and in the upper chromosphere. In the lower and middle chromosphere the bulk of the radiative losses are concentrated in the areas surrounding the strongest photospheric magnetic field concentrations but not inside the latter. The distribution of the largest temperatures is also greatly correlated with the magnetic canopy. The photospheric temperature panel shows that our FOV contains a number of small pores. We will discuss the effect of pores in §5.2. In the lower chromosphere, the peak values of the radiative losses reach ∼ −20 kW m −2 . The lowerleft panel in Fig. 8 shows that the Ca ii contribution dominates in the lower chromosphere, and the canopy shape is already visible there.
In the middle chromosphere the magnetic field becomes smoother and the magnetic canopy is clearly visible in the B image, suggesting that at this depth we are already sampling above the lower edge of the canopy. The Q middle shows a similar picture with slightly smaller radiative losses. The temperature image shows a nearly homogeneous value of approximately 6.5 kK, while most pores appear as cold holes in the canopy.
In the upper chromosphere, the integrated radiative losses are dominated by the H i contribution. In this layer the radiative losses, the enhanced chromospheric temperature and the magnetic field form a patch above the plage target with relatively constant values of < Q >≈ −22 kW m −2 , < T >≈ 8.5 kK and |B |≈ 370 G. In this layer only the strongest pores are visible in the temperature map and in the radiative losses map.

The effect of pores
The presence of pores in plage seem to have a clear imprint in the statistics of the derived physical parameters (see e.g., Chintzoglou et al. 2021). Our target contains several pores, that appear as colder and larger footpoints of the magnetic canopy than in the rest of the magnetic elements. But otherwise, the magnetic field is still strongly vertical. The imprint of pores is clear. Because they are much colder than the smaller bright flux tubes, the radiative losses are insignificant in comparison with the surroundings until we reach the upper chromosphere. Figure 10 shows a vertical reconstruction, marked with the white slit in Fig. 7. This slice cuts through two pores embedded in the plage region (at Y ∼ 8 arcsec and 8 < X < 22 arcsec), illustrating this effect. The figure also suggests that eventually, the imprint of pores is much smoother and weaker in the upper chromosphere.
Summarizing, having pores in the the FOV does affect the derived radiative losses, especially in the lower chromosphere, because the atmosphere is much colder than in regular flux tubes. Besides that, having spatially resolved maps greatly helps to separate their influence from the rest of the FOV.

Time-series analysis
The NN has made it possible to obtain the model atmosphere for all the time steps of the observation, and thereby, it has allowed us to estimate a time series of integrated radiative losses in the chromosphere over the entire FOV. We have focused on the red region highlighted in the top left panel of Fig. 7 and we have calculated the radiative energy balance for 21 time steps of the whole observations. Given the cadence of the CRISP instrument, the latter cover a total time of ∆t = 12.19 min. In §1 we mentioned that magnetoacoustic waves and shocks can have a significant contribution to the heating of plage, and their imprint should be periodic. Our aim is to separate the contribution of Ohmic heating from the contribution of waves/shocks (De Pontieu et al. 2007;Hasan & van Ballegooijen 2008) by analyzing a time series. In order to have a grasp on the dominant period (p), amplitude (A), phase (φ) and offset (C off ) of the oscillatory behaviour, we have fitted a very simple model y = C off + A · sin(φ + pt) to the v LOS and Q tot temporal curves. The results for two random pixels are shown in the left and in the middle panels of Fig. 9. The selected pixels are marked with green crosses in Fig. 7. The v LOS is estimated in the region at cmass = −3.8, corresponding to middle/upper chromosphere. The results are in line with previous studies about waves propagation in different targets of the solar atmosphere. Both curves have very similar periods, but there is a phase shift between Q tot and v LOS , which we discuss below.
In order to reduce the noise inherent to a single pixel measurement, we selected the green region of Fig. 7 and computed the average values of the parameters of the sinusoidal function. The area covers 20x20 pixels, considering that an average over a bigger region would lead to a loss of information as the temporal coherency of the curves is not achieved over large spatial scales. The two sinusoidal averaged functions are out of phase, as in the panels for the specific pixels, with a phase difference of ∆φ ∼ 138 • . The periods of the sinusoidal functions of Q tot and v LOS are respectively p Q = 5.4 min and p v LOS = 4.8 min. The offsets have been removed from the representation of the average functions for an easier comparison between the two results.
The fitting procedure has been extended to the red highlighted region in the left top panel of Fig. 7. We have plotted a map for each parameter of the sinusoidal function in Fig. 11. The displayed parameters represent the quantities characterizing the time evolution of Q, v LOS andv turb (from the top row). As a reference, average values of the parameters over the region are given in Table 1. In our model, the offset is the value of Q or v LOS not related to the periodic wave. In the case of the radiative losses, our interpretation is that it contains the contribution from other heating phenomena, such as Ohmic dissipation, ion-neutral collisions, etc. The top left panel of Fig. 11 shows a smooth distribution of the offset values. Although the smallest scales are not present in this plot, there is a large scale variation across the panel.
If we neglect the central-upper part of the maps, outside the boundary of the plage region, the average amplitudes become amp Q = 7.6 kW m −2 and amp v LOS = 3.2 km s −1 . The latter is in agreement with values reported by Centeno et al. (2009) in a facular region using the He i 10839 Å line. The relatively large value of the period in the chromosphere can be due to the propagation along an inclined magnetic field line which can extend the cut-off frequency in the chromosphere (e.g., Bloomfield et al. 2007).
The third column of Fig. 11 shows the period maps of both Q an v LOS . The yellow contours indicate those areas where the radiative losses are particularly strong (Q < −20 kW m −2 ). The average period of oscillation for Q is p Q = 5.5 min, while consid-  Fig. 11. Maps of the parameters (offset, amplitude, period and phase) of the sinusoidal functions obtained for the red region of Figure 7. The top row shows the quantities characterizing the sinusoidal time evolution of the radiative losses, while the second row represents the quantities characterizing the time evolution of v LOS at cmass = −3.8. In the maps of the period, yellow contours indicate the plage area corresponding to Q < −20 kWm −2 (see Fig. 7). The first column from the right shows the phase difference between Q and v LOS . The bottom row shows the quantities of the sinusoidal evolution for v turb at cmass = −3.8.
ering only the region inside the contour it drops to p Q = 5.2 min. These results found for v LOS are in line with previous works (de Pontieu 2004;Centeno et al. 2009), with a value of p v LOS = 5.5 min. We note that in the plage area (bottom half of the panel) the periods present a quite smooth distribution, between 4 an 5 minutes, whereas the upper area of the image, outside the plage region, shows elongated features with longer periods. The latter are located at the plage boundary, where the magnetic field is more horizontal. The phase difference between the radiative losses and the line-of-sight velocity is dominated by values close to π/2. Since the behaviour of the radiative losses is dominated by the temperature, we are essentially recovering the usual phase relation for running waves. Although we have not included plots showing the imprint of the periodic signal in each of the sublayers in which we divided the chromosphere, we have performed the fits individually for each layer. In the lower and middle chromosphere the modulation amplitude is of the order of 1 kW m −2 , whereas in the upper chromosphere we get values much closer to the integral over the entire chromosphere. Therefore, we can conclude that although the imprint of waves is present in the lower and middle chromosphere, their contribution is larger in the upper chromosphere.

Discussion and conclusions
Our temporal analysis has allowed us to quantify the contribution from a periodic component, which we associate with wave heating. This component is weaker than the offset (background) value of the heating terms and it has a mean modulation ampli-  Table 1. Average values of the parameters of the sinusoidal functions obtained over the red region of Figure 7. Each value of the table corresponds to the average of the respective panel of Figure 11. tude of ∼ 7.0 kW m −2 . This component is responsible for the very fine structure that we observe in the Q tot maps. The offset value, which we associate with a more static or very slowly evolving component, has a mean value of ∼ −26.1 kW m −2 . The map constructed from the offset value is relatively smooth, which could also point to a magnetic origin. On the Sun, the β = 2µP g /B 2 = 1 layer is usually located in the lower chromosphere, and therefore, the magnetic field becomes smooth and room-filling above that layer. Having a relatively smooth offset map signals a magnetic origin. The amplitude of the periodic component of the radiative losses is almost a factor four larger in the upper chromosphere than in the lower and middle regions. Given the spatial distribution of periods and the relatively homogeneous π/2 phase difference, we associate the periodic component with compressible acoustic waves (the slow-mode of magnetoacoustic waves when v A > c s ). This argument is further supported by Fig. 12, where we show the time-evolution of the 8542 Å line at three random pixels selected in the middle of the canopy areas that are located in the surroundings of the photo- spheric magnetic elements. In all cases, the classical saw-tooth pattern from acoustic shocks is clearly visible.
Our results show that in the lower and middle chromosphere the radiative losses are distributed in areas surrounding the photospheric footpoints of the magnetic canopy. We do not see enhanced radiative losses within the magnetic footpoints of the canopy. Furthermore, de la Cruz  showed that in those regions the Ca ii 8542 Å lines profiles have a peculiar shape that can be explained by the presence of a hot magnetic canopy in the chromosphere that extends over a relatively quiet photosphere. The magnetic canopy should have a relatively sharp lower edge, where current sheets should be found through the relation j = ∇ × B/µ. Although, our inverted models have a low depth resolution in the magnetic field reconstruction due to the S/N ratio of our observations, the temperature stratification was derived with more than twice number of nodes and we do appreciate a relatively sharp canopy boundary there. Therefore, we argue that in the lower chromosphere of plage, Ohmic current dissipation must be responsible for the bulk of the heating. We also argue that if wave heating was a dominant phenomena in the lower chromosphere, its imprint should also be visible inside the magnetic footpoints, which we do not observe in our results (see e.g., Hasan & van Ballegooijen 2008). Brandenburg & Rempel (2019) also showed that in Ohmic dissipation works efficiently in the photosphere / lower chromosphere.
In the upper chromosphere the radiative losses form a patch that covers the entire plage region, including most of the pores that are present in the photosphere. Within this patch, the derived chromospheric temperatures are relatively homogeneous and larger than in the surroundings, with temperatures the order of ∼ 7.5 − 8 kK and a magnetic field strength of the order of ∼ 400 G. The largest contribution to the integrated radiative losses is from the H i diagnostics in the upper chromosphere. The contribution of wave heating is, in our results, largest in this layer. From our results alone, we cannot claim to have direct evidence of turbulent Alfvén wave heating (van Ballegooijen et al. 2011) or from neutral-ion collisions (Khomenko & Collados 2012). Given the smooth nature of the magnetic field, we do not expect current sheets to be present in the upper chromosphere, making Ohmic dissipation of currents a less likely heating mechanism. The periodic behaviour that we observe has a relatively long period of 5.5 m. At a cadence of ∼ 30 s, wave patterns with periods lower than 1 m are not properly sampled in our observations, so we cannot resolve high frequency waves.
In order to investigate deeper the origin of the wave behaviour and the deposition of energy in the chromosphere, we briefly extended our study to to the microturbulence velocity v turb . However, the temporal analysis, represented in the last row of Fig. 11, doesn't show a clear pattern in the offset or in the amplitude. It is definitely different from the typical white noise, but doesn't correlate with the patterns shown in the panels of Q or v LOS . We could expect an imprint of the wave pattern in the microturbulence if the shocks were relatively unresolved in depth by the inversion code, forcing a larger microturbulence value to account for the extra broadening. We do not observe such behaviour.
Our results are somewhat different than those reported by Anan et al. (2021), but surprisingly also compatible. Their observations were based on arguably lower spatial-resolution slitspectrograph raster scans in the Mg ii h&k lines (IRIS) and in the He i 10830 Å line. Although their analysis was not based on the calculation of radiative losses, they used the integrated h&k line intensity as a proxy, similarly done with the Ca ii H&K lines by Leenaarts et al. (2018). We have shown that the radiative losses in the lower chromosphere are very small in the Mg ii lines. By not having the Ca ii deeper contribution, they would also miss the heating closer to the lower boundary of the magnetic canopy. As for their estimates of the chromospheric magnetic field, they are based on inversions of the He i 10830 Å line, and the latter usually samples the middle/upper chromosphere (see Fig. 1  The present study is unique in that we have calculated and studied the distribution of radiative losses as a function of depth and time from very high spatial-resolution spectra. The canonical value of Q AR ∼ 20 kW m −2 , derived in the 70's and 80's using spatially and temporally averaged spectra, could not catch the complexity and dynamic behaviour of the solar chromosphere as shown in our analysis. Our analysis is also different from previous studies in that we have estimated the contribution from waves directly from the periodic modulation of the radia-tive losses, and not by estimating from Doppler velocities the energy that is carried out by waves in the chromosphere (Abbasvand et al. 2020).
Ultimately, our results could not provide direct evidence that would allow discriminating what heating mechanism is dominating in the upper chromosphere. In our opinion, future studies of the same nature should analyze higher cadence time-series and higher spatial resolution observations in order to attempt finding observational signatures of high frequency waves or very smallspatial-scale variations that could point to turbulent Alfvén wave dissipation. In order to estimate ambipolar diffusion heating, an accurate estimate of the upper chromosphere density is also needed, and therefore inversion codes must be modified in order to include the Lorentz force (Pastor Yabar et al. 2019) and the support effects derived from velocity gradients.

Appendix A: Training of the neural network
In general, an artificial neural network is defined by the dimensionality of the input, the number of layers, the number of neurons per layer and the dimensionality of the output. The number of neurons in each layer does not have to be constant, and can vary depending on the complexity of the problem. The most used type of neural network is the fully connected network (FCN; Schmidhuber 2015), in which every input is connected to every neuron of the following layer. Figure A.1 shows, in a simplified way, the architecture and connections of a FCN. Each connection is described by a simple function that combines linearly the input x multiplied by a weight w and summed with a bias b and finally returns the value of a certain user-defined nonlinear function f (x). In mathematical notation, the information that will pass from the input neurons i to the neuron j of the next layer will be: This output will be the input for another neuron of the next layer.
Since the first operation is linear, the activation is the one that introduces the non-linear character of the FCNs. The optimization of a neural network is called training and it involves the iterative modification of the weights and biases so that a loss function that measures the ability of the network to predict the output from the input is minimized. In our case, we have trained the neural network to learn the mapping between the observed Stokes profiles and model atmosphere obtained from the inversion. Once the network is trained, we are able to reconstruct the temperature T , the line-of-sight velocity v LOS , the turbulence velocity v turb , the parallel component of the magnetic field B || , the perpendicular one B ⊥ and the azimuth angle φ for the entire time series. The dimensions of the input are defined by the total number of pixels and by the four Stokes parameters. On the other hand, the output dimensions are the total number of pixels and the number of obtained parameters multiplied by the number of grid points.
Architecture of the NN: For our purposes, we design a fullyconnected neural network with 5 layers and 200 neurons per layer. We have found this configuration to be optimal after many tests in terms of training time and accuracy. The activation function that we decided to use is the rectified linear unit or ReLU (Nair & Hinton 2010). It has a linear behavior for a positive input, otherwise, if the input is negative, it is equal to zero: It is applied after every layer of the NN, except for the last one, to avoid obtaining only positive outputs. During the design of the network architecture we detected that the noise in the polarization was amplified and propagated to some physical parameters. To avoid this problem we decided to split the model into two parts: only Stokes I was employed in the calculation of the temperature and v LOS , while all the four Stokes parameters are used for the other quantities of the model atmosphere. Although Stokes Q, U and V have information of the gradient of the source function and the line-of-sight velocity stratification, in cases where the profiles are very noisy they do not play an important role in the derivation of the temperature and line-of-sight velocity, so it is better to avoid propagating that noise. Training process and validation set: Optimization is routinely solved using simple first-order gradient descent algorithms, which modify the weights along the negative gradient of the loss function with respect to the model parameters. To scale the magnitude of our weight updates we have to use the parameter called learning rate. This parameter has to be adjusted to find a compromise between network accuracy and convergence speed. If this number is too small, it will take too long to reach the solution, while, if it is too large, there is a risk of overshooting the optimal solution. In our case we used a learning rate of 10 −4 . For the optimization method, we have used a gradient descent variant called Adam (Kingma & Ba 2017) which has been developed to automatically adjust the learning rate, making the solution convergence faster. Our goal is to optimize a user-defined loss function that evaluates how well our network models the data. The most common loss function is the mean squared error which measures the average squared difference between predictions and desirable outputs. However, to get an idea of the dispersion of the estimate we have used the quantile regression (Koenker & Bassett 1978), a loss function for estimating any percentile value: where y is the training value and f (x) is output value of the network. During training, q is randomly varied between 0 and 1 so that the network can learn all possible percentiles. This allows us not only to estimate the mean value (q=0.5) but also the dispersion (q=0.16 and q=0.84) which is equivalent to one standard deviation σ in the case of a normal distribution. The dispersion can give us an idea of the uncertainty of the inversion since similar Stokes parameters could have had different atmospheric models as solutions (Díaz Baso et al. 2021a). The gradient of the loss function with respect to the free parameters of the network is obtained using the backpropagation algorithm. Since the networks are defined as a stack of layers, the gradient of the loss function can be calculated by the chain rule as the product of the gradient of each module and ultimately of the last layer and the specific loss function. The main problem with some activation functions is that the gradient vanishes for very large values due to the derivative of this function, making it