Magnetohydrodynamic Poynting Flux Vortices in the Solar Atmosphere and Their Role in Concentrating Energy

The nature of energy generation, transport, and effective dissipation responsible for maintaining a hot solar upper atmosphere is still elusive. The Poynting flux is a vital parameter for describing the direction and magnitude of the energy flow, which is mainly used in solar physics for estimating the upward energy generated by photospheric plasma motion. This study presents a pioneering 3D mapping of the magnetic energy transport within a numerically simulated solar atmosphere. By calculating the Finite Time Lyapunov Exponent of the energy velocity, defined as the ratio of the Poynting flux to the magnetic energy density, we precisely identify the sources and destinations of the magnetic energy flow throughout the solar atmosphere. This energy mapping reveals the presence of transport barriers in the lower atmosphere, restricting the amount of magnetic energy from the photosphere reaching the chromosphere and corona. Interacting kinematic and magnetic vortices create energy channels, breaking through these barriers and allowing three times more energy input from photospheric motions to reach the upper atmosphere than before the vortices formed. The vortex system also substantially alters the energy mapping, acting as a source and deposition of energy, leading to localized energy concentration. Furthermore, our results show that the energy is transported following a vortical motion: the Poynting flux vortex. In regions where these vortices coexist, they favor conditions for energy dissipation through ohmic and viscous heating, since they naturally create large gradients in the magnetic and velocity fields over small spatial scales. Hence, the vortex system promotes local plasma heating, leading to temperatures around a million Kelvins.


Introduction
Despite considerable progress in the last few decades, the mechanisms responsible for maintaining the high coronal temperature are still not fully understood, including some critical aspects of energy transport from the lower solar atmosphere and the conversion of this energy into heat.The magnetic nature of plasma dynamics and its coupling across the solar atmosphere through magnetic field lines indicate that magnetic energy flux is the main mechanism in the energy transport (Nordlund et al. 2009;Shelyag et al. 2012).The turbulent motion in the intergranular convective zone can be considered a reservoir of energy for this purpose (Pontin & Hornig 2020).The shuffling of the magnetic field lines in the solar photosphere can stress the magnetic field quasistatically, generating direct currents (DC) or alternate currents (AC), depending on the timescale of the motion compared to the Alfvén time (Klimchuk 2006).Estimates have established that the photospheric motion produces enough magnetic energy that could, in principle, be transported to the upper atmosphere, where it is converted into heat.However, transporting that energy and transforming it into heat are aspects of both AC and DC heating that still need assessing (Klimchuk 2006;Reale 2010;Pontin & Hornig 2020).
The rate of inflow of magnetic energy flux is described by the Poynting flux, which is perpendicular to the magnetic field.Most of the Poynting flux generated by photospheric motion flows parallel to the solar surface, as shown by Silva et al. (2022).Nevertheless, the amount of upflow of magnetic energy generated by photospheric motions, although varying in time (Shelyag et al. 2012;Hansteen et al. 2015), is enough to justify the plasma temperatures in the chromosphere.The averaged X-ray brightness in active regions correlates with Poynting flux distributions, estimated to be in the range of 5 × 10 6 -4 × 10 7 erg cm −2 s −1 or 5-40 kW m −2 , providing enough energy to coronal plasma to compensate for losses and maintain the plasma temperature at observed values (Fisher et al. 1998;Tan et al. 2007).Additionally, observational estimates for solar plage regions show that the upward-produced Poynting flux also carries the required energy input (Yeates et al. 2014;Welsch 2015).The inflow of magnetic energy into the upper solar atmosphere can be produced by horizontal magnetic fields transported by vertical plasma motions or the work done by horizontal motions on the vertical magnetic field.The primary process to create upflows of energy will depend on the local magnetic configuration of the solar atmosphere.Shelyag et al. (2012) and Yadav et al. (2021) showed that the vertical Poynting flux is mainly created by plasma motion parallel to the solar surface doing work on the vertical component of the magnetic field.Based on realistic MHD simulations, Steiner et al. (2008) concluded that the expulsion of magnetic fields from the upper photosphere can also provide the vertical component of the photosphere's Poynting flux.
Flow or kinetic vortices (K-vortices), which are defined as vertical rotating columns of plasma flow (see, e.g., Silva et al. 2020), are also believed to be an essential mechanism to create magnetic energy input to compensate for chromospheric losses (Shelyag et al. 2012;Wedemeyer-Böhm et al. 2012;Yadav et al. 2021;Tziotziou et al. 2023).In K-vortices, the vertical Poynting flux is still mainly created by work done by the horizontal plasma motion on the vertical magnetic field encompassed by the vortex tubes (Yadav et al. 2021) and intensified around the vortex tube boundaries (Kitiashvili et al. 2012).Even in the chromosphere, swirling motions can also generate an upward-directed Poynting flux (Battaglia et al. 2021).The solar atmosphere also supports the magnetic vortices (M-vortices), coherent magnetic flux tubes where the magnetic twist is present for a considerable part of the tube, usually rooted along intergranular lanes (Silva et al. 2021).The M-vortices can be interpreted as torsional Alfvén perturbations propagating along the magnetic field (Battaglia et al. 2021).The M-and K-vortices are different structures, as their centers and boundaries are localized at different positions and they appear and disappear at different time frames.This is clear when looking at Figure 18 from Tziotziou et al. (2023), for example.Their coexistence leads to a structure defined as a magnetic tornado, described by Wedemeyer-Böhm et al. (2012).The K-vortices can also excite and act as a waveguide to Alfvén waves, as found in MHD simulations (see, e.g., Battaglia et al. 2021;Finley et al. 2022) and observations (Jess et al. 2009;Tziotziou et al. 2020).By interacting with magnetic flux tubes, the vortical motions create enough energy to be transported upward and even more effectively produce Poynting flux than other horizontal motions (Yadav et al. 2021;Chian et al. 2023).Coherently rotating magnetic field structures, the magnetic tornadoes, are found to produce enough Poynting flux to justify chromospheric and coronal temperatures (Wedemeyer-Böhm et al. 2012), and their presence leads to an enhancement of the energy flux by four times compared to time intervals in which there are no magnetic tornadoes (Kuniyoshi et al. 2023).
Most studies concerning the heating processes of the upper solar atmosphere have focused on restricted parts of the whole process, such as estimating the kind of heating (AC or DC), determining the amount of energy generated due to the shuffling of magnetic field lines by photospheric motions, and assessing the plasma response to heating (see, e.g., Klimchuk 2006;Reale 2010).In contrast, our study concentrates on mapping the properties of the flow of energy as it travels through the solar atmosphere.
A way to adequately describe the transport properties of the flow in fluid mechanics is to apply the methodology of Lagrangian Coherent Structures (LCSs; Haller 2015).Recently, this method has been successfully implemented for solar plasma environments and helped to shed light on some essential properties of flows.Chian et al. (2019) showed that the magnetic field distribution is organized by the attracting and repelling regions of the flow.Roudier et al. (2021) proposed the use of combined Lagrangian methodologies for investigating large-scale events.Importantly, LCS theory provides the necessary framework for tracking energy transport during the lifetime of coherent plasma structures, even for plasmas in a highly dynamic and nonlinear state, such as intergranular or supergranular turbulence, or, regarding large structures, solar flares or coronal mass ejections.Our study presents an LCS analysis applied to the Poynting vector field to show how magnetic energy is transported in the solar atmosphere.Based on the realistic quiet-Sun simulation obtained by the Bifrost code, our analysis reports the actual energy source and deposition locations by identifying the LCS in the energy flow and the existing barriers to magnetic energy transport in the solar atmosphere.In addition, we analyze how the rotation of the plasma flow can produce considerably steep gradients at very small spatial scales, a necessary condition for significant heating in a plasma of small dissipation coefficients.

Bifrost Data
The numerical data were obtained using the Bifrost code (Gudiksen et al. 2011), which realistically simulates the upper and lower convection zones of the solar photosphere, extended into the entire atmosphere reaching the low corona.The simulated domain consists of 768 × 768 × 768 points, covering 24 Mm in the horizontal directions and 16.8 Mm in the vertical one.There is a slow inflow of the horizontal field at the lower numerical boundary, located 2.5 Mm below the simulated visible solar surface.The left panel in Figure 1 shows the whole solar surface, z = 0.0 Mm (corresponding to where the optical depth is unity at a wavelength of 500 nm), at t = 4200 s, colored by the vertical component of the velocity field.The numerical domain extends vertically about 14.3 Mm above the simulated surface, comprising the lower corona.The model used to represent the corona is a representation of a region with open fields in the photosphere.To compute radiative losses in the transition region and chromosphere, scattering and non-local thermodynamic equilibrium conditions are applied (Carlsson & Leenaarts 2012).Further information on this particular simulation data set, identity number ch024031_by200bz005, is presented by De Pontieu et al. (2021).
We analyzed a part of the domain that covers 9.4 × 9.4 Mm in the horizontal plane and is delimited by a green square in Figure 1(a).The close view of the analyzed region is shown in Figure 1(b), where we show the simulated surface colored by B z and the xy-plane at z = 7.0 Mm colored by v z .The velocity field lines are shown in blue and their seed points were randomly selected from inside the region where a large-scale K-vortex is detected.We define a large-scale vortex as a coherent structure presenting both a sizable radius (greater than 1 Mm) and height (ranging from 10 Mm).A discussion on how we identify the vortex tubes is presented in Appendix A. In the vertical direction, part of our analysis is limited to z 5.0 Mm (marked in purple in the right panel), covering the photosphere, chromosphere, and low corona.

Mapping the Energy
Our study focuses on describing the magnetic energy transport based on the Poynting flux, S. In the MHD approximation, the Poynting flux can be written as Vector S describes the direction in which magnetic energy is flowing and its magnitude provides the rate per unit area at which the energy crosses a surface.However, just analyzing the Poynting flux at each time frame does not provide relevant information, such as where the barriers to energy transport are and how the plasma dynamic is limiting or even producing more outflow of magnetic energy.To construct the energy map, i.e., identify the transport barriers and the outflow and inflow regions, we will compute the forward/backward Finite Time Lyapunov Exponent (FTLE) field of the energy velocity, v e , defined by Heaviside (1922) as: where u m = B 2 /8π is the magnetic energy density.For a discussion on the energy velocity, see Schantz (2018).The calculation of the FTLE fields in 3D is carried out by advecting energy parcels initially placed along the original mesh of horizontal and vertical slices of the numerical domain.The advection of the energy parcels is solved by the following equation: Equation (3) is solved over a grid of initial positions r 0 , until the final positions r(t 0 + τ) are reached after a finite-time duration τ.To solve Equation (3), we used a fourth-order Runge-Kutta integrator and applied cubic splines for space and time interpolations using a discrete set of energy velocity field frames.The FTLEs of the particle trajectories for a 3D flow are calculated at each initial position r 0 as where l max is the maximum eigenvalue of the finite-time right Cauchy-Green deformation tensor Δ = J  J, ( ) is the deformation gradient,  denotes the transpose of a tensor, and ) is the flow map for Equation (3).In other words, to compute the FTLE, the flow map is obtained by integrating Equation (3), forward or backward.Then, for every grid point i, j, k, we compute the Jacobian of the flow map, i.e., the deformation gradient, J, by calculating how much the initial distances from the grid point at i, j, k to its local neighboring points have changed after the grid and neighboring points were advected during τ time units.Based on that, the right Cauchy-Green deformation tensor, Δ, is computed for each grid point.Then, the largest eigenvalue of Δ is obtained and the FTLE is computed from Equation (4), resulting in the mean amount of stretching experienced during τ time units by a fluid parcel moving in the direction of maximum stretching from that grid point.The SI unit of the FTLE is s −1 .Regions presenting contraction instead of stretching display negative values for the FTLE field.
The lines of local maxima in the forward-time FTLE (f-FTLE) field represent regions, or ridges, of maximum local stretching in the flow.In stationary flows, the stable manifolds of saddle points will produce such ridges, for two particles in the close vicinity of the stable manifold, but on opposite sides of it, they will approach the saddle for some time.When they are near the saddle point, they will suddenly diverge from each other, following different branches of the unstable manifold (see Figure 4 of Shadden et al. 2005); in time-dependent flows, the manifolds are dynamic structures with finite time and correspond to the LCSs.In the following sections, the ridges of the f-FTLE are called repelling regions, since they act as sources or outflow regions in the flow.The ridges of the backward-time FTLE (b-FTLE) field represent lines of maximum local contraction and, therefore, are called attracting regions.They are analogous to the unstable manifolds of steady dynamical systems.In our context, the ridges of f-FTLE for the energy velocity can be interpreted as sources of magnetic energy flux and the attracting regions can be interpreted as the deposit locations of magnetic energy flux.More information on how to interpret the f-FTLE and b-FTLE ridges as LCSs can be found in the studies by Haller & Yuan (2000) and Shadden et al. (2005).
Figures 2 and 3 show the f-FTLE and b-FTLE fields calculated using the 3D energy velocity superposing the Line Integral Convolution (LIC; Cabral & Leedom 1993) of the magnetic field (first column) and the velocity field (second column) at the initial time of the integration, t 0 , for parcels initially placed at different height levels along xy-slices.LIC is a technique that allows us to visualize, through darker lines, the points tangent to the vector field along the plane.The FTLE fields were normalized by their maxima and the nonzero values are not shown.To more clearly visualize the ridges, the FTLE was made transparent for values below 60% of their maximum.The value chosen for normalization was the maximum value of FTLE (forward or backward) among the three selected initial times t 0 .The FTLE fields in the horizontal xy-planes were calculated for a portion of the 9.4 Mm ×9.4 Mm domain, depicted by the green square in Figure 1.For our analysis, we choose the xy-planes placed at z = 1.28 Mm (Figure 2) and z = 8.0 Mm (Figure 3).Other height levels above the photospheric range present either no special features or results similar to the selected xy-planes (see Figures 13 and 14 in Appendix B).The integration time, τ, was chosen to be 50 s, as it was the time interval that provided the best results for the FTLE analysis, i.e., a time smaller than the necessary time for the particles to leave the domain, yet large enough to provide well-localized LCSs, i.e., thin structures.The red square indicates the region where two large-scale vortices, a kinetic (K) and a magnetic (M) one, will be created by the plasma dynamics.Those vortices appear, become coherent, and later decay during the time of the analysis, as can be noticed from the LICs of both the magnetic and the velocity fields.The panels in the first row show the FTLE fields for integration starting at t 0 = 3800 s, when the large-scale K-vortex and M-vortex are present in the lower part of the domain, as indicated by the LIC pattern inside the red square.The middle row displays the panels of the FTLE fields for t 0 = 4200 s; the K-vortex and M-vortex tubes are well developed, i.e., the tubes are connecting the simulated surface to the top boundary of the domain and the vortex pattern can be observed at all height levels.The third row depicts the FTLE fields for t 0 = 4450 s, when those vortices start to divide into two smaller vortical structures at the upper boundary and disappear or lose strength closer to the surface, z 1.28 Mm.
The intersection areas of backward and forward FTLE are a complex network of homoclinic, heteroclinic, and saddle points, and they are the least attracting regions for magnetic energy, indicating the regions where the transport can occur across the ridges of FTLE fields.Therefore, the region where the large-scale vortices appear presents a complex behavior in terms of magnetic energy transport: concentrating energy in some parts, displaying the outflow of magnetic energy at other points, and allowing the mixing of energy from different sources through the saddle points.Figure 2(b) shows the ridges of f-FTLE encompassing the K-vortex region (the spiraling LIC in the center of the red square), indicating that the region initially acts as an outflow of magnetic energy.Later, for the analysis starting at t = 4200 s (Figure 2(d)), the K-vortex is surrounded by a semicircular inflow/deposition location for magnetic energy transport.At t 0 = 4450 s, panels (e) and (f) show that once the K-vortex is not present inside the M-vortex at that height level (z = 1.28 Mm), the inner ring formed by the ridges of FTLE is replaced by repelling lines.Thereby, both in their beginning and decay, the vortices seems to act as sources of energy.In general, we can notice that the ridges of b-FTLE seem to be spatially connected, creating a network of energy.
Figure 3 shows that most of the FTLE ridges correspond to the outer parts of the K-vortices or regions between vortex flows, suggesting that the magnetic energy is generated by vortical motion and deposited along thin localized regions.The FTLE ridges also act as energy transport barriers and the fact that the sources and deposition locations are so close to each other indicates that the energy does not travel long distances across the horizontal plane.Moreover, once the K-and M-vortices are developed (Figures 3(c) and (d)), large circular structures are formed, acting as both depositories and outflows of the magnetic energy within the shown red square frame.Some of the FTLE ridges change from attracting (backward) to repelling (forward) for the same initial integration time.This change can be interpreted as the energy flowing to the location (the b-FTLE ridges) during the integration time up to t = t 0 and then part of that energy being transported from that area to another part of the domain, acting as an outflow region (f-FTLE ridges).By comparing panels (a) and (b) of Figure 3      parts of the vortex tubes presented hotter plasma than their surroundings; in contrast, the plasma inside the vortex tubes is colder in the higher region.
In order to analyze changes in the energy transport in the vertical direction and the role of the detected vortices in establishing an energetic channel, we selected three xz-planes: The FTLE provides a map of the energy distribution, but it lacks information on the direction of energy propagation.To establish the energy transport path followed during a time interval, we computed the pathlines of the energy velocity, Equation (2), calculated for τ = 50 s starting at t 0 = 4200 s.We select this initial time integration as this corresponds to the time when the M-and K-vortices are established.The pathlines are shown in blue in Figure 7

Discussions and Conclusion
Our results present an unprecedented map of magnetic energy sources and deposition regions in a quiet solar atmosphere simulated by the Bifrost code.To map and describe the properties of the energy flow, we expanded a cutting-edge methodology for flow analysis, the FTLE, to the energy flow velocity, which is calculated as the Poynting flux divided by the magnetic energy density (see, e.g., Schantz 2018).The regions of the solar atmosphere acting as an outflow or deposition of the magnetic energy flow were identified by the ridges of the forward and backward FTLE fields, respectively, as displayed in Figures 2, 3,  and 6.It was thereby possible to assess the impacts on magnetic energy transport due to the presence of large-scale K-and M-vortices, showing that those vortical structures completely change the energy mapping.In particular, the area where significant vortices emerge exhibits intricate behavior regarding the transport of magnetic energy, showing that the vortex tubes act as both energy sources and depositions during their lifetimes.As the M-and K-vortices interact, the Poynting flux follows a vortical motion, cospatial with the K-vortex, as displayed in Figure 4.This creates a vortex of magnetic energy flow through the solar atmosphere, as indicated in Figure 7.The FTLE analysis shows that the Poynting flux vortex interior initially acts mostly as a source of magnetic energy, i.e., the energy coming from the lower part of the domain will be transported from the vortex interior to other parts of the domain, such as the vortex boundary regions, which attract part of the energy input being transported by that vortex.Later, the energy map is changed, and there are attracting regions to magnetic energy inside the vortex.Therefore, part of the energy remains in the vortex, leading to the higher Poynting flux distribution in the region of vortices than before, as illustrated by Figure 8.
The ridges of the FTLE fields are barriers to energy transport (Haller 2015); i.e., no energy flow can cross the sources and deposition locations indicated by the yellow/orange and blue filaments in Figures 2, 3, and 6.Therefore, the interaction of the K-and M-vortices breaks the existing energy transport barriers across the solar atmosphere, as depicted in Figure 6, which prevents most of the energy generated by surface motions from reaching the upper atmosphere.Once the transport barriers cease, a significant amount of energy flow can reach the upper part of the simulation domain, as illustrated in Figure 8, where the volume rendering of S log 10 | | in the region depicted by the purple box in Figure 1 is shown at different moments.The animation associated with Figure 8 shows this evolution.Initially, t 4000 s, most regions presenting intense Poynting flux are only found close to the simulated surface.After the energy channel is established, Figure 8(c), intense Poynting flux is observed for the full vertical extension of the vortex, with a stronger concentration around the M-vortex boundary region.The Poynting flux distribution starts to break and concentrate at the M-vortex center close to the vortex decay stage, Figure 8(d), confirming the findings provided by the FTLE fields presented in Figure 6(c).
To assess the role of the vortex tubes in changing the energy input and plasma heating in the upper part of the domain, we computed the average energy input S log 10(| |), temperature T log 10( ), and plasma density ( r log 10( )) for each height level at every 50 s for (a) the vortex region and (b) the central region of the domain with the same size, where no large-scale vortices are detected in the given time interval.For both Figures 9(a) and (b), most of the energy is generated in the photosphere or in the vertical range 0.5 z 1 Mm.In the region where the vortices appear, their presence leads to three times more Poynting flux reaching the upper part of the domain and also more Poynting flux being generated in the lower part of the domain.The average temperature distribution is displayed in the middle panel of Figure 9.As the vortices are established, around t = 4100 s, the plasma inside the vortex area is cooled down, due to input of denser plasma in the upper region, as indicated in the last panel of Figure 9(a).The upper part of the vortex tube, at z 12 Mm, is on average heated up to more than 6 × 10 5 K.As the vortices start their decay around t = 4300, the plasma inside that region becomes less dense and hotter.No considerable changes are observed for the region where no large-scale vortices are present, as illustrated by Figure 9(b).
One necessary condition for the plasma heating is the conversion of the excess magnetic energy indicated in Figure 9 into heat.As the coronal plasma has considerably small dissipation coefficients, larger gradients of magnetic and velocity fields are required.In order to evaluate the changes due to the M-and K-vortex dynamics in creating conditions for plasma heating, we computed the squared current density, j 2 , and the square of the norm of the velocity gradient, ||∇v|| 2 , at different times.These quantities are proportional to the ohmic and viscous heating rates.Figure 10 shows xz-planes for the selected domain from Figure 1 placed at y = 10.94Mm, crossing the central region of the vortices at three different times: (a) t = 3850 s (vortex tubes are not established), (b) t = 4200 s (vortex tubes are fully developed), and (c) t = 4450 s (start of the vortex tube decay).The animation associated with Figure 10 depicts this evolution.In each row ((a), (b), and (c)), the first panel is colored by the line-of-sight (LOS) component of the velocity field, v y , which is superposed by the 3D vortex core of the Poynting flux computed according to the λ 2 -criterion (see Appendix A).To depict the K-vortex region, we use the velocity field crossing the panel, v y .A vortical motion will be crossing the xz-plane by "entering" and "outing" that panel, creating a pattern presenting coherent strong positive and negative v y in that xz-slice, thereby helping to identify the region in which K-vortices are present.From   10(a) show no significant differences for the analyzed quantities in the vortex region for the upper part of the domain, z 4 Mm.Once all the vortex tubes are established, we see from Figure 10(b) that in the upper part of the simulation, z 10 Mm, there are localized areas with higher temperatures in the region where the vortices coexist.For the height range 4 < z 7 Mm, the central region of the vortex presents higher temperatures than its boundary and cooler temperatures compared to the plasma surrounding the vortex tubes.That tendency is reversed for the lower region of the vortices, in the range z 4 Mm.Thus, depending on the atmospheric layer, the vortex tubes can concentrate either hotter or cooler plasma temperatures compared to their surroundings.The existence of large-scale vortex tubes also leads to high-velocity gradients and intense magnetic energy fluxes around the vortex boundary.In the vortex region, the current density is more substantial than in other parts of the upper atmosphere.In the last frame, corresponding to t = 4450 s, the interior of the vortex was heated up to even higher temperatures (around one million Kelvins), and in the upper part of the domain, the contributions from |S| and ||∇v|| 2 are more significant than before, concentrated along the interior vortex region.A comparison between frames at t = 3850 s and t = 4450 s indicates that the vortices can create intense currents and velocity gradients in the form of filamentary structures.These structures present short transversal scales relative to the vertical lengths where stresses build up that may lead to heating.This is confirmed by the localized high-temperature regions that appear in the second panel.
Although previous studies (see, e.g., Wedemeyer-Böhm et al. 2012;Yadav et al. 2021;Kuniyoshi et al. 2023) had already demonstrated the potential role of K-vortices for energy transport, our analysis has shown for the first time how the magnetic energy is flowing in the upper atmosphere and how the presence of vortices shapes the energy distribution.Our results strongly underline the effectiveness of magnetic energy mapping through FTLE in elucidating energy transport characteristics within the solar atmosphere.In summary, through this mapping, we have identified, for the first time, distinct energy transport barriers that suppress most of the flow of magnetic energy from the photosphere to the upper atmosphere.This indicates that observational estimates on the upflow of magnetic energy generated by photospheric motion might overestimate the energy that reaches the upper atmosphere.Our findings also reinforce the role of vortices as energy channels, as they breach these energy transport barriers, resulting in a significant augmentation of the magnetic energy availability in the presence of the vortices, around a three times higher Poynting flux than before, as quantified by Figures 9 and 10 and similarly obtained by Kuniyoshi et al. (2023) for a magnetic tornado.After the vortices were established, the input of the magnetic energy into the upper atmosphere by the Poynting flux vortex was around 10 22 erg per second, which is the same amount of energy required by the nanoflare heating scenario (Pontin & Hornig 2020).Moreover, our investigation shows that the magnetic energy generated by the M-and K-vortices is continuously transported upward as a Poynting flux vortex.The magnetic energy mapping also reveals that the transported energy follows a dual behavior, with a portion becoming concentrated in the vortex region and another portion being transported upward into the higher layers of the solar atmosphere.Finally, our analysis shows that the vortical dynamics imposed by K-and M-vortices results in large gradients of magnetic and velocity fields over small scales, conditions that favor the dissipation of magnetic energy and, thereby, plasma heating (Klimchuk 2006).
The Q-criterion is defined as ( ) ] ¯is the rate-of-strain tensor, and ∇v is the gradient of the velocity field.Therefore, vortex regions are defined by Q crit > 0.0 and regions with negative Q crit tend to be dominated by strain.From Figures 11(c) and (d), we see that the lower part of the vortex is indeed dominated by vorticity as Q crit > 0.0, but the upper part of the vortex has regions where Q crit is close to zero.In those regions, the strong asymmetry of the rotating flow velocity causes weaker vorticity, and therefore the square norms of spin and strain tensor present similar values, leading to Q crit presenting negative or close to zero values in the vortex region.The weaker vorticity in a flow presenting strong strain makes it challenging to identify the vortex using IVD and other methodologies that require vorticity-dominated flows.
In order to identify the region of vortex tubes, we proceed to use LIC in combination with the vortex core lines.When horizontal sections are taken across a vertical vortex tube, circular patterns of LIC will be present in each section.The core line, situated at the vortex's center, will consequently be connecting the circular LIC regions observed on the xy-planes at different height levels.Therefore, in this study, we define the vortex tube region as the location where a vertical core line, identified by applying a threshold to the λ 2 -criterion (Jeong & Hussain 1995), intersects regions with circular patterns in the LIC field in xy-planes for a given height range.Although in this simulation the effects of compressibility are important, it does not seem to affect the identification of the vortex in our analysis.The λ 2 -criterion defines the vortex core as the region where the second largest eigenvalue of [Ω 2 + E 2 ], λ 2 , is negative.The core lines of K-vortices were identified using the threshold λ 2 − 0.001.The identification of the vortex tubes' region at t = 4200 s is depicted in Figure 12(a), where four xy- planes displaying the LIC of the velocity field are placed at z = 2, 5, 8 and 11 Mm.The black line is the core line identified by the λ 2 -criterion and the blue lines were traced from 100 random seed points within the vortex region that is visible from the LIC.For the M-vortices, we followed the same procedure to keep consistency.The identified M-vortex tube is depicted in Figure 12(b), and in this case the LIC depicts the magnetic field orientation along the four xy-planes.The magnetic field lines are shown in red and were traced from 20 random points along an LIC line within the vortex region at z = 11 Mm.As the magnetic field lines are less twisted in the lower part of the domain, the core line of the M-vortex is not as precise as in the K-vortex.It is important to notice that IVD only failed to identify the large-scale vortex at z = 8.0 Mm; the other small vortices present at the same height level are well defined by IVD in Figure 11.Therefore, the strain rate is likely influenced by the interaction with the M-vortex that is intersecting the K-vortex, as indicated by Figure 12, and it is more twisted in the upper part of the domain.

Appendix B FTLE Superposing LIC for Magnetic, Velocity, and
Poynting Flux Fields In this section, we present the complementary images to the analysis shown in Figures 2 and 3.The regions of outflow for magnetic energy transport are identified by the ridges of f-FTLE (orange/yellow lines) and the deposition locations are provided by the ridges of b-FTLE (blue lines).Figure 13 displays the FTLE superposed by the LIC of the magnetic field at z = 0.65 Mm (left column) and z = 3.14 Mm (right column).Figure 14 shows the same results, but for the LIC of the velocity field.

Figure 1 .
Figure 1.Simulation domain of the Bifrost identity number ch024031_by200bz005 at t = 4200 s.(a) The simulated solar surface of the Bifrost numerical domain covers 768 × 768 pixels or 24 × 24 Mm 2 .The xy-plane is colored by the vertical component of the velocity field and the green square shows the selected 300 × 300 pixels region (≈9.4 × 9.4 Mm 2 ) used in our analysis.(b) Close view of the selected region.The simulated surface is showing the z-component of the magnetic field in grayscale and the plane placed at 7 Mm above the surface is colored by the vertical velocity.The blue lines show the velocity field lines inside the K-vortex.The purple cuboid depicts a region analyzed in other parts of the paper.
with the panels where the vortex is present, (c) and (d), we can notice that the full development of K-and M-vortices considerably changes the mapping of the energy transport.The b-FTLE ridges within the red-framed areas in Figures 2 and 3 reveal that energy deposition follows a somewhat circular trajectory within the region containing both an M-and a K-vortex.The f-FTLE ridges also exhibit a similar geometric pattern.Consequently, the FTLE fields suggest that the Poynting flux moves along a vortical path, since the Poynting flux converges to attracting b-FTLE ridges in the long run.A portion of this flux is deposited along this path, while the remaining energy continues to trace the vortical motion of the plasma in an ascending or descending manner.This phenomenon results in the formation of the observed semicircular patterns in the FTLE fields along the horizontal slice cut crossing the vertical M-and K-vortex tubes The presence of a Poynting flux vortex is confirmed in the left panel of Figure 4, which displays the magnetic (red field lines), kinetic velocity (cyan field lines), and Poynting flux (purple field lines) vortices at t = 4200 s.The field lines of the M-vortex encompass the kinetic and Poynting flux vortex all the way from z = 3 Mm onward.In the upper part of the domain (z 10 Mm), the energy and K-vortices disengage as the velocity fields lose their rotation and become more vertical.The white square indicates the close views shown in the right panels of Figure 4, depicting the LIC colored by the log(T) for the magnetic field (B), velocity (v), and Poynting flux (S), respectively.In the upper right panels ((b)-(d)), we show the LIC for z = 8 Mm, and in the bottom right panels ((e)-(g)), we display the LIC for the lower atmosphere at z = 1.28 Mm.A comparison between the

Figure 2 .
Figure 2. LICs of the magnetic field (first column) and velocity field (second column) superposed by the ridges of FTLE fields calculated at z = 1.28 Mm above the simulated surface for the initial integration times t 0 = 3850 s ((a) and (b)), 4200 s ((b) and (d)), and 4450 s ((e) and (f)).The regions of outflow for magnetic energy transport are identified by the ridges of f-FTLE (orange/yellow) and the deposition locations are provided by the ridges of b-FTLE (blue).The red square delimits the region of appearance of the M-and K-vortices.

Figure 3 .
Figure 3. LICs of the magnetic field (first column) and velocity field (second column) superposed by the ridges of FTLE fields calculated at z = 8 Mm above the simulated surface for the initial integration times t 0 = 3850 s ((a) and (b)), 4200 s ((c) and (d)), and 4450 s ((e) and (f)).The regions of outflow for magnetic energy transport are identified by the ridges of f-FTLE (orange/yellow) and the deposition locations are provided by the ridges of b-FTLE (blue).The red square delimits the region of appearance of the M-and K-vortices.

Figure 4 .
Figure 4. Identified vortex tubes at t = 4200 s in the analyzed region.(a) Close view of the domain, showing the field lines for the magnetic field (red), velocity (cyan), and Poynting flux (purple) traced from random seed points inside the vortex region.The simulated surface is colored by the vertical component of the velocity field and the planes at z = 1.28 Mm (8 Mm) are colored by T log 10( ).((b)-(g)) Close views of the xy-planes indicated by the white square in (a) colored by T log 10( ) superposed to the LIC of the magnetic field (first column), velocity field (second column), and Poynting flux (third column) at z = 8 Mm (upper panels (b), (c), and (d)) and at z = 1.28 Mm (bottom panels (e), (f), and (g)).

Figure 5 .
Figure 5.The selected domain from Figure 1 up to z = 5.5 Mm, colored by the vertical component of the velocity field.The black lines depict the selected xz-planes to compute the FTLE in vertical planes.
and the initial grid position is indicated by the black dots.In the upper and lower panels, we display the pathlines of energy parcels initially distributed in the xz-plane placed at y = 10.94Mm (upper panel) and in the xy-plane placed at z = 3.14 Mm (lower panel).Both panels indicate that the energy is mostly transported in the horizontal direction, except mainly for the vortex region, where the pathlines indicate a vertical component of the energy velocity.Therefore, as indicated previously in Figure 6, the vortex region facilitates a greater flow of energy to the upper part of the simulated atmosphere.An animation of the Poynting flux vector at t = 4200 is available with Figure 7.The animation shows the view from the top of the analyzed region.This animation of streamlines of the Poynting flux at time t = 4200 s indicates how the magnetic energy would flow if the Poynting flux was constant in time.The seeds for the streamlines are

Figure 7 .
Figure 7. Particles advected (in blue) using the energy velocity starting at the points shown in black.An animation showing the stationary flow along the Poynting flux vector at t = 4200 is available for the upper view.In the animation, random points in the domain create streamlines colored by H z , which is defined as the height divided by the unit Poynting flux vector in the z-direction.Warm and cold colors indicate upward and downward flow, with yellow/cyan for the near-surface flow and red/dark blue for the upper-domain flow.The x-and y-axes are labeled in Mm.The real-time duration of the animation is 43.6 s. (An animation of this figure is available.)

Figure 8 .
Figure 8. Poynting flux, S log 10 | |, in volume rendering, within the purple cuboid of Figure 1.(a) Time t = 3850 s, when no major vortices are present yet.(b) t = 4000 s, when M-and K-vortices start to appear.(c) t = 4200 s, when vortices are established throughout the atmosphere.(d) t = 4400, when the vortex starts to decay.An animation of this figure is available.It starts at t = 3850 s and ends at 4450 s, with 50 s intervals.The real-time duration of the animation is 6.5 s. (An animation of this figure is available.)

Figure 9 .
Figure 9. Averages of magnetic Poynting flux, temperature, and plasma density calculated at each height level and at every 50 s for the 3 × 3 × 14.3 Mm 3 region (a) indicated by the red square in Figure 2, where the large-scale vortex tubes are formed around t = 4000 s, and (b) in the central part of the domain, where no large vortices are present.

Figure 10 .
Figure 10.Vertical xz-plane of the selected domain from Figure 1 crossing the central part of the vortex region at y = 10.92Mm.Columns from left to right: LOS component of the velocity field, v y , with opacity=0.3superposed by the 3D S-vortex center; log 10 of temperature; the square norm of the velocity gradient, ||∇v|| 2 ; log10 of the squared current density; and log10 of the magnitude of the Poynting flux.The rows from top to bottom refer to the times t = 3850, 4200, and 4450 s.An animation of this figure is available.It starts at t = 3850 s and ends at 4450 s, in 50 s intervals.The real-time duration of this animation is 6.5 s. (An animation of this figure is available.)

Figure 11 .
Figure 11.Vortex identification methods.IVD ((a) and (b)) and Q crit ((c) and (d)) calculated at z = 1.28 Mm and 8 Mm for the region indicated in Figure 1(b).The white contours in panels (a) and (c) indicate the vortex boundary defined by the outermost convex contour of the IVD field.
Figures 15 and 16 show the FTLE fields superposing the LIC of the Poynting flux at z = 0.65 Mm, z = 1.28 Mm, z = 3.14 Mm, and z = 8.0 Mm.The red square indicates the region of the appearance of the large-scale vortex under investigation.

Figure 12 .
Figure 12.Illustration of the vortex tubes' detection using the LIC of the vector field (displayed in the xy-planes) and the core lines.The LIC was computed for (a) the velocity and (b) the magnetic field and is shown here for the xy-planes at z = 2, 5, 8, and 11 Mm.The core lines are shown in black and were identified by the λ 2 -criterion.The vector field lines were traced from regions within the vortex area defined by the LIC.

Figure 13 .
Figure 13.The LIC of the magnetic field superposed by the ridges of the FTLE fields calculated at z = 0.65 Mm (left column) and z = 3.14 Mm (right column) above the simulated surface for the initial integration times t 0 = 3850 s ((a) and (b)), 4200 s ((c) and (d)), and 4450 s ((e and (f)).The regions of outflow for magnetic energy transport are identified by the ridges of f-FTLE (orange/yellow) and the deposition locations are provided by the ridges of b-FTLE (blue).The red square delimits the region of appearance of the M-and K-vortices.

Figure 14 .
Figure 14.The LIC of the velocity field superposed by the ridges of the FTLE fields calculated at z = 0.65 Mm (left column) and z = 3.14 Mm (right column) above the simulated surface for the initial integration times t 0 = 3850 s ((a) and (b)), 4200 s ((c) and (d)), and 4450 s ((e) and (f)).The regions of outflow for magnetic energy transport are identified by the ridges of f-FTLE (orange/yellow) and the deposition locations are provided by the ridges of b-FTLE (blue).The red square delimits the region of appearance of the M-and K-vortices.

Figure 15 .
Figure 15.The LIC of the Poynting flux superposed by the ridges of the FTLE fields calculated at z = 0.65 Mm (left column) and z = 1.28 Mm (right column) above the simulated surface for the initial integration times t 0 = 3850 s ((a) and (b)), 4200 s ((c) and (d)), and 4450 s ((e) and (f)).The regions of outflow for magnetic energy transport are identified by the ridges of f-FTLE (orange/yellow lines) and the deposition locations are provided by the ridges of b-FTLE (blue lines).The red square delimits the region of appearance of the M-and K-vortices.

Figure 16 .
Figure 16.The LIC of the Poynting flux superposed by the ridges of the FTLE fields calculated at z = 3.14 Mm (left column) and z = 8.0 Mm (right column) above the simulated surface for the initial integration times t 0 = 3850 s ((a) and (b)), 4200 s ((c) and (d)), and 4450 s ((e) and (f)).The regions of outflow for magnetic energy transport are identified by the ridges of f-FTLE (orange/yellow) and the deposition locations are provided by the ridges of b-FTLE (blue).The red square delimits the region of appearance of the M-and K-vortices.