Nonlinear wave damping by Kelvin-Helmholtz instability induced turbulence

Magnetohydrodynamic kink waves naturally form as a consequence of perturbations to a structured medium, for example transverse oscillations of coronal loops. Linear theory has provided many insights in the evolution of linear oscillations, and results from these models are often applied to infer information about the solar corona from observed wave periods and damping times. However, simulations show that nonlinear kink waves can host the Kelvin-Helmholtz instability (KHi) which subsequently creates turbulence in the loop, dynamics which are beyond linear models. In this paper we investigate the evolution of KHi-induced turbulence on the surface of a flux tube where a non-linear fundamental kink-mode has been excited. We control our numerical experiment so that we induce the KHi without exciting resonant absorption. We find two stages in the KHi turbulence dynamics. In the first stage, we show that the classic model of a KHi turbulent layer growing $\propto t$is applicable. We adapt this model to make accurate predictions for damping of the oscillation and turbulent heating as a consequence of the KHi dynamics. In the second stage, the now dominant turbulent motions are undergoing decay. We find that the classic model of energy decay proportional to $t^{-2}$ approximately holds and provides an accurate prediction of the heating in this phase. Our results show that we can develop simple models for the turbulent evolution of a non-linear kink wave, but the damping profiles produced are distinct from those of linear theory that are commonly used to confront theory and observations.


INTRODUCTION
Observations show that the solar atmosphere is filled with highly structured plasma forming loops and threads of material tracing the magnetic field.A number of different phenomena lead to oscillations of these structures in the direction transverse to the field.A few examples are the flare-generated transverse coronal loop oscillations (Aschwanden et al. 1999;Nakariakov et al. 1999); the more prevalent small-scale disturbances in active region loops (McIntosh et al. 2011); Dopplershift disturbances in extended regions of the solar corona (Tomczyk et al. 2007); propagating transverse waves in prominences (Lin et al. 2007;Schmieder et al. 2013); or the occurrence of transverse waves generated from colliding plasma flows (Antolin et al. 2018).A common interpretation has been given to these oscillations, in terms of propagating or standing transverse magnetohydrodynamic (MHD) kink waves (see e.g., Ruderman & Erdélyi 2009;Goossens et al. 2011;Nakariakov et al. 2021, for comprehensive reviews).
The usual paradigm under which theoretical studies on MHD kink waves have been carried out is that of waves modelled as linear perturbations to an ini-tial background state (Roberts 1983(Roberts , 2000)).Theoretical models for the damping of kink waves have also focused mainly on the linear regime (Goossens et al. 1992(Goossens et al. , 2002;;Ruderman & Roberts 2002;Goossens et al. 2006).However, as evidenced by the catalogues compiled by Goddard & Nakariakov (2016) and Nechaeva et al. (2019), many of the observed kink oscillations have amplitudes that are large and their damping seems to depend on the oscillation amplitude.For at least some of the observed oscillations, the nonlinearities associated with the perturbations to the system are non-negligible.This can lead to damping of the oscillations through nonlinearities (e.g.Chen & Schuck 2007;Van Doorsselaere et al. 2021;Arregui 2021).
A well-established result from analytical analysis and numerical simulations is that plasma motions in a flux tube undergoing a nonlinear kink oscillation can lead to the development of the Kelvin-Helmholtz instability (KHi) (see e.g.Terradas et al. 2008;Antolin et al. 2014Antolin et al. , 2015;;Magyar & Van Doorsselaere 2016) and the subsequent turbulence the instability can induce (e.g.Hillier & Arregui 2019;Hillier et al. 2020).In this case, the Kelvin-Helmholtz instability is a parasitic instability that grows on the shear-flow that exists on the flanks of the oscillating flux tube.If the instability can grow, it is then able to develop nonlinearities and if there is sufficient energy in the flow it can then develop turbulence (Hillier 2019;Hillier & Arregui 2019).This turbulence can extract energy from an oscillation, either providing a saturation mechanism for the amplitude of the wave for a driven oscillation (Hillier et al. 2020) or lead to damping of an impulsively excited mode.The turbulence excited by the Kelvin-Helmholtz instability is one that is fundamentally different from MHD wave turbulence (e.g.van Ballegooijen et al. 2011), where nonlinear MHD waves interact to create a daughter wave of higher frequency.The key difference being that even though both of these mechanisms use the large-scale oscillation as the energy source, in the wave turbulence model the large-scale oscillation is also involved in creating the energy cascade process.For Kelvin-Helmholtz turbulence this is not the case, with the energy cascade being related to the nonlinearities of an instability.
In this paper we perform a detailed analysis of the evolution of the Kelvin-Helmholtz instability on the surface of an oscillating flux tube, identifying how the turbulent dynamics results in two different phases of the evolution of the oscillation amplitude of the tube.We then develop an analytic model for the first stage of the evolution of the flux tube.Here we focus on the large scale response of the oscillating tube to the Kelvin-Helmholtz turbulence, in particular the damping of the oscillations and the heating this creates.During the second phase of the dynamics we find that the turbulent energy dominates that of the wave energy, so we develop a model only for this based on classic models of decaying turbulence, using these to predict the heating rate in the latter stage.

SIMULATIONS OF KINK-WAVE-DRIVEN KELVIN-HELMHOLTZ INSTABILITY AND RELATED WAVE DAMPING
To build a model of wave damping through Kelvin-Helmholtz turbulence, it is necessary to first identify the fundamental processes that are occurring.To do this we perform a 3D ideal MHD simulation of a nonlinear kink oscillation that develops the Kelvin-Helmholtz instability.Through this we look to identify what it means for the Kelvin-Helmholtz instability to damp a kink wave, and look at the fundamental processes involved in the damping.

Simulation setup
We perform this ideal MHD simulation of an impulsively excited MHD kink wave using the MHD routines of the (PIP) code (Hillier et al. 2016).We solve the evolution in 3D in a Cartesian reference frame using the non-dimensionalised ideal MHD equations in conservative form, namely: which allows us to calculate the evolution of the primitive variables, i.e the density (ρ), velocity field (v = (v x , v y , v z ) T ), pressure (P ) and magnetic field (B), through the evolution of the relevant conserved quantities.Note that in our normalisation we have taken the magnetic permeability of a vacuum inside B meaning that the local Alfvén speed is given by V A = |B|/ √ ρ.
These equations are solved using a fourth-order central difference approximation for spatial derivatives (calculated on a uniform mesh) with a four-step Runge-Kutta time integration.We have non-dimensionalised the equations using a characteristic coronal density (ρ c ), sound speed (C s ) and lengthscale (L c ).
Here we do not include any explicit resistive or viscous terms in the equations.Dissipation at some level is inherent in the simulation due to the finite grid size and the use of flux limiters to smooth sharp structures.As the code is written in conservative form, any extraction of energy from the flow or magnetic field results in a corresponding increase in internal energy (i.e.heating).This allows us to quantify the magnitude of any heating that occurs in the simulation whilst being able to run calculations in the least viscous, least resistive regime we can achieve for the resolution.
The initial conditions are of a dense tube, aligned with the direction of the magnetic field, which we then perturbed to excite the fundamental kink mode of the system.The initial density profile is given by with r = x 2 + y 2 , R the radius of the tube, and ρ i and ρ e the internal density of the tube and the external density.We take R = 0.3, ρ e = 1 and ρ i = 3.The initial pressure is uniform throughout the domain taking a value P (x, y, z) = 1/γ (which gives an initial sound speed of C s = 1) and a magnetic field of B = (B x , B y , B z ) T = 2/(γβ)(0, 0, 1) T , with β = 0.05.The calculations are performed in a domain of x ∈ [−L x , L x ], y ∈ [0, L y ] and z ∈ [0, L z ] with L x = L y = 1 and L z = 10.We use a grid size of ∆x = 0.005, ∆y = 0.0025 and ∆z = 0.1.The anisotropic resolution is chosen to give most resolution to the structures in the x-y plane, with greater emphasis on the y direction to help the growth of the initial instability.This anisotropic resolution in the x-y plane will result in an anisotropic numerical diffusion of turbulent structures, where the diffusion associated to the larger grid-scale will likely dominate the diffusion of the system as structures are rotated in that plane by the turbulent motions.This tube is perturbed with a lateral velocity perturbation to the tube in the form to excite a fundamental kink mode of the system.We use a value of V 0 = 0.2.We also set the y-and zcomponents of the velocity to zero (i.e.v y = v z = 0).
We set the boundaries of our domain to be either periodic or possessing symmetries.The boundaries at x = −L x and x = L x are set to be periodic.The boundaries at y = 0 and y = L y are symmetric boundaries such that the magnetic field is imposed to be in the (x, z)-plane, as is the flow.This implies that v y = B y = 0 on the boundary and ∂v x /∂y = ∂v z /∂y = ∂B x /∂y = ∂B z /∂y = ∂ρ/∂y = ∂P/∂y = 0.The z-boundaries are symmetric with the magnetic field penetrating the boundary with the imposed symmetries on the magnetic and flow field being of those of a node (v x = v y = v z = 0 and ∂B x /∂z = ∂B y /∂z = ∂B z /∂z = ∂ρ/∂z = ∂P/∂z = 0 at z = 0) and anti-node (B x = B y = v z = 0 and ∂v x /∂z = ∂v y /∂z = ∂B z /∂z = ∂ρ/∂z = ∂P/∂z = 0 at z = 0 at z = L z ).

The evolution of the kink wave and their damping through KHi-induced turbulence
Having been subject to its initial kick, as explained in the previous subsection, the tube begins to oscillate.The wave excited is (or at least dominated by) a fundamental kink mode with frequency ω KINK = 0.53 (calculated from linear theory in the long-wavelength limit), resulting in the dense tube oscillating back and forth.For our initial conditions we have set a sharp boundary in density between the tube and the external medium (see Equation 7).Therefore, these oscillations are not subject to resonant damping, which requires a continuous non-uniform variation of density such that the fun-damental kink mode has its frequency in the Alfvén continuum (Van Doorsselaere et al. 2004;Goossens et al. 2006).However, as the natural state for these oscillations is for a strong shear flow to develop between the tube and the external medium (Sakurai et al. 1991;Goossens et al. 1992), the boundary of the oscillating tube becomes unstable to the Kelvin-Helmholtz instability (Terradas et al. 2008;Antolin et al. 2014Antolin et al. , 2015;;Magyar & Van Doorsselaere 2016).
Figure 1 shows the density evolution of the tube crosssection at the wave apex (the x − y plane at z = L z ).These snapshots are taken every 3.2 time units to show the evolution of the oscillation at approximately every quarter period.As the oscillations proceed, we can clearly see the development of Kelvin-Helmholtz rollups on the top of the tube.After approximately the first quarter period (t=3.2) the instability has started to grow, though the scales associated with this (and with that the area of the cross-section of the tube that has been disturbed by the development of the instability) are small.As the wave continues to oscillate the Kelvin-Helmholtz vortices become larger, and with this the area of the cross section that has been disturbed by the turbulence keeps on increasing.Ultimately, the vast majority of the tube becomes part of the turbulent layer.
Looking at this in further detail, Figure 2 shows the temporal evolution of the area of cross section of the tube at the apex (normalised by the initial value) for different density thresholds.We plot the density ρ ≥ 2.7 (shown by the red line), ρ ≥ 1.89 (shown by the black line) and ρ ≥ 1.1 (shown by the blue line).The area shown by the red curve acts as a proxy for the area that has not been disturbed by the Kelvin Helmholtz induced turbulence.It is clear that once KHi develops there is a clear evolution of material that is initially part of the kink oscillation of the tube but then joins the mixing layer.Once the undisturbed area of the tube at the apex becomes about 20% of its initial area, we can see that the decrease in size of the tube core slows.The area of material with ρ ≥ 1.89 stays approximately constant over time.This is an important result as it connects to the model of Hillier & Arregui (2019) where for the density contrast used in this calculation the density value of ρ = 1.89 is predicted to be an area conserving threshold, so is predicted to show no evolution.
Looking back at Figure 1 we can see by eye that the total displacement of the tube, at least as measured by the displacement along the x-axis, does not appear to be showing a strong decay with time.We can see this somewhat clearer in Figure 3 which shows the spatial evolution of v x at y = 0 and z = L z over time.The blue  and red colours show the positive and negative velocities respectively.Until a time of t ≈ 44, the magnitude of the displacement and the magnitude and coherency of the velocity field do not show any drastic change.However, after that time, corresponding to where the decrease in area for the density threshold of ρ ≥ 2.7 is clearly reduced in Figure 2, there is a complete change in the dynamics with clear reductions in the lateral displacement of the tube and a velocity field which is less coherent and of smaller magnitude.
The temporal evolution of the velocity field in the plane of the apex provides further details on the evolution from coherent oscillations to turbulent motions.Figure 4 shows the evolution of v x in the plane at the apex of the oscillation (z = L z ).In this plot we have added two black lines to approximately denote the inner and outer edges of the turbulent layer using the density thresholds of ρ = 2.7 and ρ = 1.1 respectively.Though this is more pronounced initially, throughout the evolution the x velocity for the region inside the inner black line (i.e. the core of the tube) is relatively coherent at all times.However, in the turbulent layer (i.e. the region between the two black lines) there is more significant fluctuation around the average motion.This includes local speed-ups in the flow which are a common feature of Kelvin-Helmholtz roll-ups (e.g.Hasegawa et al. 2006) and other turbulent structuring.The coherent component of motions in the turbulent layer clearly do not have to move with the core of the tube, where sometimes they both move in the same direction but sometimes they are moving in opposite directions.
We can quantify the difference in the velocity field between the core and the layer by looking at the RMS (Root-Mean-Squared) values of the x velocity and the y derivative of the x velocity in each of the regions.By dividing the former by the latter we can get an approximation of the lengthscale over which the flow is evolving.In the core this value is consistently around 10, implying the average variation of the flow happens on scales larger than the core.However, for the mixing layer this value is ∼ 0.02, consistent with a flow varying over small length scales.
There is one question that Figure 4 (along with Figures 1 and 3) elicits related to measuring how the wave damps.We can see that at least for the first few periods the core of the tube undergoes coherent oscillations, but the whole evolution of the tube takes into account the turbulent layer which is growing over time, has more incoherent motions and a coherent component of the motions that moves differently to the core.Therefore, it is necessary to ask: How do we measure the coherent motions of the loop?And are all these measures the same?
Clearly the way we measure the bulk motion is determining how we view the change in the oscillation amplitude with time.We can see this very aspect in Figure 5 where the centre-of-mass displacement in the x-direction has been calculated in three slightly different ways resulting in very different displacements measured.For the solid black line, the displacement in the 2D plane at the apex for all the material above a density threshold of ρ ≥ 1.62 is displayed.The solid red line is a similar calculation, but for a density threshold of ρ ≥ 1.1 and the blue line for a density threshold of ρ ≥ 2.7.We can see that initially the amplitude of the oscillation, as given by the maximum and minimum displacements, decreases monotonically for all cases.However, for the different density thresholds the rate at which the amplitude reduces (and the form of the envelope giving the damping of the oscillating profile) clearly differs.We can also see that there is a clear period drift for the cases with lower thresholds compared to the ρ ≥ 2.7 curve.After t ≈ 44 the oscillations are difficult to distinguish in the lower threshold curves, but are clear and consistent in the ρ ≥ 2.7 curve that captures the motion of the core of the tube.There are three important pieces of information that can be taken from this.Firstly, the density thresholds are reasonably accurate at capturing the extent of the turbulent layer at the apex, but not at the foot point.This can be understood dynamically as the KHi grows at the apex, with flows that transport the density field.These flows, often vortical, will excite waves that travel along the magnetic field and stress the magnetic field at the foot point.However, the fixed boundary means that there are no flows to move the mass.So at the foot point the density thresholds there do not correspond to the width of the turbulent layer (it is more the magnetic mapping from the apex to the foot point that would provide this).Secondly, it can be seen that at this time, the energy is greater in the mixing layer than the core, and the energy fluctuations in the mixing layer are of small scale (as expected from the lengthscale analysis presented above).Thirdly, the area with turbulent kinetic energy seems to be slightly larger than the region of turbulent magnetic energy, but this is balanced by the turbulent magnetic energy generally having larger magnitude.
Ultimately this leads us to a fundamental question: How do we quantify damping of a kink oscillation when  the loop develops a growing turbulent layer?The standing kink wave is a coherent oscillation of a dense tube aligned with magnetic field.The damping of these oscillations, when considering linear wave theory, should be the same wherever it is measured, not heavily dependent on how the tube is being measured.Clearly something different is happening in this simulation where not only the amplitude but the period of the oscillations measured is a function of the density threshold used to calculate a centre-of-mass amplitude evolution.
Another key aspect of this simulation is that there are two clear phases of the dynamics.The initial period of damping oscillations connected to the growth of a turbulent layer, followed by the later stage when the tube cross-section is almost completely turbulent and the growth of the layer is almost completely arrested.In the next two Sections we will present models for these two phases by developing and combining aspects of the models in Hillier & Arregui (2019) and Hillier et al. (2023), and bench-marking the predictions of the model with key aspects of the simulated evolution.

MODELLING THE DEVELOPMENT OF A TURBULENT LAYER AROUND AN OSCILLATING TUBE
In this section we develop a formulation to describe the evolution of the kink oscillation into a turbulent tube, i.e. the first stage of the damping as described in the previous section.A model for key aspects of the second stage will be presented in Section 4.
The first stage of the damping process happens as the turbulent layer on the boundary grows.This process exchanges momentum from inside and outside of the tube, manifesting as the originally coherent oscillations of the tube developing a more incoherently moving outer layer, with a coherently moving inner region.This can be seen in Figure 5 where the oscillations of the centre-of-mass for the material with density greater than 1.1 has an amplitude decaying faster than just the material with density greater than 2.7.
In this regime, we can expect the velocity jump across the turbulent layer stays relatively constant because the velocity at the centre of the tube (as evidenced by Figure 3) remains at approximately the same magnitude until t ≈ 44.If we consider a simpler flow, i.e. a nonoscillating hydrodynamic shear layer, it is well known that for a constant velocity jump across the layer the thickness (h) of the turbulent layer grows as (e.g.Winant & Browand 1974) where β MIX is a constant for a given mixing layer but the particular value will depend on the density contrast of the layer (Baltzer & Livescu 2020).Using a dimensional analysis of the model proposed in Hillier & Arregui (2019) and developed in Hillier et al. (2023) we propose that the mixing layer grows as where C 1 is a constant of between ∼ 0.1 and ∼ 1.Through comparison with mixing in hydrodynamic shear flows, the value of C 1 is expected to be in the range of C 1 = 0.3 (Baltzer & Livescu 2020) to 0.5 (Brown & Roshko 1974) as discussed in Hillier et al. (2023).To determine an appropriate value for the mixing constant C 1 for this problem we will look more at the mixing in this section.By doing this we can then use this model of an expanding shear layer, combined with information on the structure of the shear layer taken from Hillier & Arregui (2019), to calculate the rate at which mass and momentum are transferred into or out from the oscillating loop and with that how the oscillation damps.Before making these comparisons, it is necessary to measure the magnitude of the shear flow across the turbulent layer as this is what drives the Kelvin-Helmholtz instability-induced turbulence.Figure 7 shows the shear flow at the surface of the tube measured at x = 0 and z = L z at t = 11.6 (blue line) and t = 18 (red line).The velocity is reformulated such that the value is positive at y = 0 for easy comparison.Two dashed horizontal lines are added at v x = 0.15 and v x = −0.05 to show the approximate values of the flow inside the dense tube and outside the tube.This implies the magnitude of the shear flow is ∆V ≈ 0.2.We can see that locally the shear may appear to be larger than this value, but this is a consequence of the KHi vorticies driving local speedups in the flow (e.g.Hasegawa et al. 2006), which can be seen in Figure 4, and not a change in the large-scale velocity shear that feeds the turbulence.
Some consideration has to be given to the fact that this is an oscillatory flow and as such will not have the same magnitude of shear flow at any given time.We hypothesise that as the turbulent motions drive the growth in layer width, the RMS shear flow speed is the appropriate value to use, therefore we redefine h with a factor of 1/ √ 2, Note that this is not strictly necessary as the constant C 1 will have to be calculated and the value calculated will just scale up or down based on the numerical constants included in the model.To make the value of C 1 as close to unity as possible (i.e.explicitly including as many of the physical processes that determine its value as possible), we include this factor here.

Mass evolution
To study the evolution of the tube oscillations in Section 2.2 we looked at the centre of mass motions of the system based on given density thresholds (e.g. Figure 5).However, these different thresholds result in changes in the mass over time, so a first important step is to determine how the mass above a given threshold evolves with time.
As for this model, the actual detailed wave dynamics are not so important for the construction of the model we have decided for simplicity of the model to take a rectangular cross-section of length 2R and initial height H = Rπ/4 to model the initial cross section of the tube to simplify the calculation of the evolution of the dynamics in the analytical model.The simple justification of this is that this maintains the cross-sectional area, mass, momentum and kinetic energy of the tube whilst making analytic progress easier.Figure 8 gives a diagrammatic representation of the model tube cross-section (panel b) and it evolution as a consequence of the growth of mixing layers (panel c).With the model geometry and following the self-similar argument, we can predict the evolution of the total mass as where m ρ>ρ T is the mass above a given density threshold and ∆m ρ>ρ T is the change in mass above that density threshold for a mixing layer of unit width as predicted by the model of Hillier & Arregui (2019).
Based on the model predictions of Hillier & Arregui (2019), ∆m ρ>ρ T is a constant in time, therefore our model becomes Using the model of Hillier & Arregui (2019), we can calculate the value of ∆m ρ>ρ T .This is done by first calculating the initial mass of the dense region before mixing by multiplying the density of the high density region (for our simulations ρ = 3) by the fraction of a unit length that would have been occupied by this dense phase pre-mixing (i.e. the amount the mixing layer has moved encroached into the the original height of the tube area).Then we integrate the density component of the model of Hillier & Arregui (2019) over a unit length from the density threshold value to the maximum density value (in this case 3).This is plotted in Figure 9 for different thresholds for the initial density contrast of our simulation.We can see for a threshold that is both small enough and greater than 1 the mass will be growing in time, but for larger thresholds there is a loss of mass.The value where the mass evolution is roughly neutral, i.e. ∆m ρ>ρ T = 0, is ρ T = 1.62.Equation 13can be integrated to give (14) To confirm that Equation 14gives the expected evolution of the mass, we compare the evolution of the mass for three different density thresholds from the simulation with comparisons to the predicted evolution given by Equation 14. Figure 10 shows the evolution of the mass at the apex for three density thresholds: ρ T = 1.1 (blue line), ρ T = 1.62 (black line) and ρ T = 2.7 (red line).The dashed lines are the predicted evolution curves for these three thresholds (during this initial regime of the evolution) using C 1 = 0.3, which was determined by fitting by eye to the ρ T = 1.1 lines.The value of C 1 found here (C 1 = 0.3) is in the range of the values found in hydrodynamic simulations/experiments (Brown & Roshko 1974;Baltzer & Livescu 2020).
Overall the simulated curves follow the predicted linear trend.It is clear that though there is some evolution in the total mass for the density threshold of ρ T = 1.62, this is a sufficiently accurate threshold to use to main- thresholds are employed.The reason the ρ T = 1.62 threshold shows some evolution could be due to errors in applying this Cartesian model for the mixing to this more complex geometry, with the extreme density thresholds (that capture either almost all or almost none of the mixing layer mass) likely to be less sensitive to this effect as the mapping between the mass associated with the density value in the two geometries being probably more accurate though further work is needed on this.Other, well-established reasons, unrelated to the mixing process, like the ponderomotive force driving mass accumulation at the loop apex (e.g.Terradas & Ofman 2004), exist as a potential alternative explanation.Though our calculations show that the mass accumulation at the apex by the ponderamotive force can only account for 0.3% of the total mass evolution, which is at a level where it can be treated as insignificant.It is worth noting that the density threshold of ρ T = 1.89, which is predicted to be area conserving by the model of Hillier et al. (2019), is approximately area conserving for the duration of the simulation (as evidenced in Figure 2) One consequence of having determined the value for C 1 is that it allows an upper bound for the time this model holds to be determined.That is to say we can calculate the time when the mixing layer has become large enough that all of the initial cross-section of the tube is predicted to be engulfed in the turbulent mixing layer.Using the asymmetry of the mixing layer predicted by Hillier et al. (2019), the whole cross-section of the tube should have become turbulent when where as stated previously H is half the initial height used for the rectangular approximation of the tube given as H = πR/4.This implies that at a time of the whole tube would have become turbulent and, therefore, the model of a growing turbulent layer is invalid.This predicts that a second stage to the dynamics must have been reached by this time, which is exactly what is seen in Figure 3.

Momentum evolution
Following on from the model of the mass evolution, we can develop a similar (though slightly more complex) model for the evolution of the momentum.To do this we split the momentum above a given density threshold into two components: 1) the momentum held in the material still corresponding to the original tube (M core ), and 2) the momentum injected into the mixing layer (M L ).This is the same concept as has been used for the modelling of the mass evolution, but with the added complexity that the momentum changes sign as the system undergoes its natural oscillations.The combination of a global oscillation at the kink frequency with the momentum extraction/injection process acts to make the momentum extraction (from the core) or injection (into the layer) oscillatory.Therefore, the forcing in the system driving the momentum change is oscillatory, meaning that it is natural to attempt to approximate this system by a forced linear oscillator.
Firstly, to model the evolution of the momentum in the core, we know that the core of the tube is oscillating at the frequency of the kink wave (ω KINK ) and the momentum extraction is driven by the shear-flow dynamics which extract momentum at the same frequency.Therefore, we need to look at the case where the forcing occurs at the same frequency (i.e.resonant frequency) as the natural oscillations of the system.This can be stated in equation form as i.e. a forced linear oscillator where the system is being forced at a resonant frequency.Here M core is the momentum of the core, we define I core such that dI core /dt = M core and a forcing ( Ḟcore ) given by with M 0 the initial momentum density given by where V i is the speed of the dense material, in this case V i = 0.15 as can be seen in Figure 7.That is to say, Equation 19gives the rate at which momentum is lost from the core as a consequence of the mixing layer growing at a constant rate.Equation 17 has the solution If we follow our initial assumption that the forcing term is in phase with the oscillations, this implies that B = − Ḟcore /ω 2 KINK , which leads to This leads to the momentum evolution of the core following (22) To model the momentum injected into the mixing layer requires more subtlety.Unlike the mass evolution, the momentum is not a positive definite quantity, it oscillates around zero, so the sign of the momentum injected into the mixing layer changes periodically.On top of this, the local field lines in the mixing layer have their own characteristic frequency of oscillation (i.e.Alfvén frequencies), which can be different from that of the frequency with which momentum is being injected.
The simple model we put forward to explain how the momentum in the mixing layer develops is based again on a forced linear oscillator, but this time the characteristic frequency of the oscillation of the layer and that of the core are not assumed to be the same.We consider a single, representative Alfvén frequency for the mixing layer (a composite Alfvén frequency or kink frequency of the layer) based on the turbulent motions coupling all the different regions of the layer allowing it to develop oscillations at a single representative frequency.As the model of Hillier & Arregui (2019) proposes, the representative density of the layer is √ ρ i ρ e .Therefore, the simple estimate of the non-dimensional frequency would be ω A = B/(ρ i ρ e ) 1/4 .However, as the inverse of the square-root of a mean is not the same as the mean of the inverse of a square-root, we calculate the latter from the model of the density distribution across the layer from Hillier & Arregui (2019) and use this to calculate the approximate Alfvén frequency of the mixing layer (ω A ≈ 0.6).Here we model the whole layer as a single forced nonlinear oscillator with momentum M L .Mathematically, the evolution of the momentum would be modelled as: where the forcing is occurring at the kink frequency (ω KINK ).We define I L such that dI L /dt = M L and ḞL is the forcing term relating to the rate at which momentum is added into the mixing layer from the external regions.This is given by with M as the momentum injection calculated from the model of Hillier & Arregui (2019) over a region of unit thickness.The precise calculation is given by The term in the bracket gives the mean velocity of the mixing layer in the rest frame of the simulation.This is then multiplied by the integral of the density between the two density limits of interest across the predicted average density profile of the mixing layer model of Hillier & Arregui (2019) for a layer with width of unit length.Equation 23 leads to the well know solution (in the case that the characteristic Alfvén frequency of the layer is not the same as the kink frequency of the tube) (26) If we take that at t = 0 we have I L = 0 this implies that Equation 26 can be differentiated to give the momentum evolution in the mixing layer as for the simulation (solid) and model (dashed).The red lines show the evolution for the quantities where 1.1 < ρ < 2.7 for the simulation (solid) and model (dashed).

We then have
Figure 11 shows the comparison between the predicted momentum and kinetic evolution of the core of the tube (M core ) and the mixing layer (M L ) compared with the simulation results.What is presented here is not a fit, but the solution from an initial value problem compared with the results from the 3D simulations.This model clearly captures the key features of the evolution.
Going further, we can simplify Equation 29 by taking M L (0) = 0 and looking at the case where ω A − ω KINK is small.By setting ω D = (ω A −ω KINK )/2 and ω S = (ω A + ω KINK )/2, which means ω A = ω S + ω D and ω KINK = ω S − ω D , we have This can be expanded out through double angle formulas to give Taking that ω D ≪ ω S this can again be simplified to This implies we expect our momentum in the mixing layer to oscillate with a frequency of ω S inside an envelope that evolves as sin(ω D t).
Our result above implies that we can refine our prediction for when this initial phase of the dynamics will end.The model we have put forward for this period of the dynamics is that of the oscillations in the core of the tube driving the motions.However, it is clear from panel (b) of Figure 11 (showing the kinetic energy evolution of the core, mixing layer bulk motions and the turbulent fluctuations) that we reach a time where assuming that the core of the tube holds the majority of the kinetic energy associated with coherent motions no longer holds.At a time of t ≈ 33 the energy of the coherent motions of the mixing layer becomes greater than that of the energy of the oscillations of the core.So even without including the energy of the turbulent motions we can see that our assumption that the core energetically dominates the dynamics only holds for a limited time.This switch from the core to the layer dominating the energetics might explain why the model is less accurate once we reach a time of t ≈ 40.

Predicted evolution of amplitude and velocity amplitude of the wave motions
With a prediction for both the momentum evolution and the mass evolution above a given density threshold of the oscillation, this can then be turned into an estimate of the centre-of-mass velocity of the oscillation.This can be formulated mathematically as with ρ T the threshold density used and a the area of the integration.That is the integral of the density weighted velocity divided by the integral of the density gives a centre-of-mass velocity, but this is just the momentum divided by the mass.As both m ρ≥ρT and M ρ≥ρT have direct predictions as shown in the previous subsections, we can use these to subsequently make a prediction of V CoM,ρ≥ρT .Panels (a) and (c) of Figure 12 show the temporal evolution of centre-of-mass velocity for just the core of the tube at its apex (panel a) and when the mixing layer is also included (panel c).Unsurprisingly, given the accuracy of the momentum and mass evolution up until a time of t ≈ 40 these panels show the model is a good representation of the simulation results.
A measure of the evolution of the velocity of the loop apex over time can be used to make a prediction for the evolution of the position of the loop apex (i.e. the wave amplitude) over time.Then by integrating V CoM,ρ≥ρT over time, we can make a subsequent prediction for the amplitude evolution (A ρ≥ρT ) of Due to the nonlinear system we are studying, though dimensionally the same, the A ρ≥ρT defined above is different from the centre-of-mass amplitude which is defined as This difference in formulation is very important to consider when comparing results for a nonlinear system, as differing formulations will lead to some differences in results.This is seen in panels (b) and (d) of Figure 12, where for looking only at the core of the tube (with its simpler evolution) the model prediction for the amplitude evolution (A ρ≥ρT ) is a reasonable representation of the centre-of-mass amplitude from the simulation (A CoM,ρ≥ρT ).However, in panel (d) there is a much greater decrease in the centre-of-mass amplitude (not seen in the centre-of-mass velocity) compared to the model.

Bounding the heating rate in the first phase of the dynamics
An important question, especially in terms of how the results of this paper may impact our understanding of wave heating of the solar corona, is: can we predict the heating rates in the simulation from the model of Kelvin-Helmholtz turbulence.As we have performed an ideal MHD simulation, we do not have any explicit heating terms which to calculate this, but the numerical dissipation will (due to the use of a conservative scheme) result in energy lost from the kinetic and magnetic energies being added to the internal energy.Combined with the choice of boundary conditions which mean energy cannot leave the simulation domain this gives us some measure of the dissipative heating from the turbulence.The black line in Figure 13 shows the increase from the initial value of the total internal energy over time.At the end of this phase (around t ∼ 40) we find a total internal energy increase in the simulation volume that is ≈ 0.1% of the initial internal energy of the simulation domain.
To make a prediction of the rate at which the internal energy increases (i.e. a heating rate), first we need to make an estimate of how the magnitude of the turbulent energy (both in the velocity and the magnetic fields) held in the mixing layer evolves over time.To do this, we need to account for all the energy that has not been previously accounted for in Section 3.2, i.e. energy not in the bulk motions of the layer.
This energy can be distributed into three categories: the turbulent energy, the energy in the variation of the mean and the energy lost by the forcing of the mixing layer not being resonant with the layer's natural frequency.The first of these is simple to understand, it is the energy that is predicted to be in turbulent fluctuations by the model of Hillier & Arregui (2019).The total energy of this component is calculated as (Hillier & Arregui 2019) where a factor of 1/2 has been introduced as a consequence of integrating along the loop assuming the velocity shear depends on z as sin(πz/(2L z )).The second of these is not seen as turbulent energy by Hillier & Arregui ( 2019), but has to do with the mean flow having a distribution across the mixing layer.This enters the calculations here as the mixing layer is no longer moving with the shear driving but is drifting out of phase.As such it is not clear how much of the energy we expect to be in the mean variation is kept there and how much might be released for further turbulence.The total en- ergy of this component is calculated as Finally, we have the energy that is lost by forcing the mixing layer at a frequency different to its natural oscillatory frequency.This is just the difference between the energy of the bulk motions of the mixing layer when they are forced at the kink frequency and forced at the natural frequency of the layer.This difference is more pronounced for larger density contrasts highlighting a greater potential for heating in those situations.As this calculation would already take into account the crosssectional area, this would be multiplied by L z to give a total energy E force .
The total possible energy in the turbulence in the mixing layer can then be calculated by summing these three: i.e.E tot = E turb + E meandist + E force .The evolution of E tot is shown with the dashed red line in Figure 13.This can be used as an expected upper limit from turbulent heating in the simulation.
We can add more nuance to these arguments.We should expect a time delay between energy being brought into the mixing layer and its dissipation after a turbulent cascade.To model this, we can use the approximate self-similarity of the turbulent energy cascade which implies that the velocity fluctuations at a given length scale (l) scale as ∝ (l/l 0 ) 1/3 where l 0 is the largest scale of the turbulent cascade or integral length scale (Kolmogorov 1941).Assuming that the turbulent cascade involves the standard nonlinear process where energy is passed from one scale to the scale of half that size and that the turbulence cascades to infinitesimal scales (the Reynolds and magnetic Reynolds number are tending to infinity) leads to the following estimate of the dissipation time τ DISS (e.g.Onsager 1949) ) where the subscript denotes the level in the cascade and τ EDDY = l 0 /v 0 is the large-scale vortex turnover time.This leads us to estimate that at any given time the energy that has been brought into the mixing layer and has formed part of the turbulent fluctuations, τ EDDY /τ DISS = 1/2.7 of the energy will have been dissipated leaving 1.7/2.7 as part of the turbulence.
This leads to an estimation of a lower bound for the heating rate of (E turb + E force )/2.7.This is shown with the red solid line in Figure 13.Overall the two bounds we derive do provide bounds for the growth of the internal energy of the simulation.The fact this is closest to the upper bound may imply that the turbulent energy is very efficiently dissipated (in much less than 2.7 turnover times) or more likely as the energy from the initial kick is not all trapped in the oscillation some extra heating is occurring.Figure 14.Kinetic energy of the x component of the velocity at the apex of the oscillation for the core of the tube (black line), the coherent motions of the mixing layer calculated from the mean density and density weighted velocity of the layer (red line), with the energy of everything left, i.e. the incoherent, turbulent motions (blue line).

MODELLING ENERGY DISSIPATION IN THE SECOND PHASE OF THE DYNAMICS
Having developed a model to explain the evolution of the first phase, we now turn our attention to the second phase of the dynamics.The key characteristic of this regime that we use to guide the model is that the energy of the dynamics is dominated by the turbulent component (see Figure 14 after t ≈ 40).It is also important to note that as no further energy is input into the system and without an energy source the turbulence will decay over time (again see Figure 14 after t ≈ 40).Therefore, we focus on developing a model for the decay of the turbulence.
To model the decay of the turbulence, we will no longer consider the energy of the oscillations and restrict our arguments solely to the evolution of the turbulence.This simplification makes the model we propose for understanding the rate at which the turbulence decays to be that of simple decaying turbulence, and independent of the oscillatory dynamics.The model we use is based on those first proposed by Taylor (1935) and Kolmogorov (1941) where the turbulent transport takes the turbulent energy held at large scales to smaller scales until it is dissipated.This can be understood simply by the nonlinear turbulent transport through spatial scales, which can be approximated by: where E turb is the energy held in the turbulent fluctuations of the velocity and magnetic field and v turb are the turbulent motions that transport the energy to smaller scales.
The RHS of Equation 39 can be approximated by using the root-mean-square (RMS) of the turbulent motions (v RMS ) and the largest lengthscale of the turbulence (also known as the integral lengthscale) and further by connecting to the energy of the turbulent fluctuations to get (e.g.Kolmogorov 1941;Onsager 1949) ) is an estimate of the connection between the root-mean-squared velocity of the turbulent motions and the turbulent energy with A(t) the area of the turbulent region and L(t) the integral length-scale of the turbulent motions.Figure 2 shows that at late times the area of the turbulent region is approximately constant implying that A(t) is constant (which we denote as A 0 ).If the area is not changing then we can also assume that the integral lengthscale L(t) is also constant, which we denote as l 0 following Section 3.4 Separating the variables in Equation 40 and then integrating leads to the following integral equation: Therefore, with Based on the area of the mixing layer shown in Figure 1 we take A 0 to be an annulus with outer radius of 0.45 and inner radius of 0.15.This leads to a lengthscale l 0 which we take to be the difference between the inner and outer radius of 0.3.For dynamical consistency we take that the unknown constant C has the value D/2.7 where the 2.7 comes from the cascade time for Kolmogorov turbulence (see Section 3.4), i.e.C = D/2.7 = 0.3/2.7.This then gives The simple implication from this model is that the turbulent energy should decay at late times ∝ t −2 .This is a classic result for decaying turbulence with a fixed lengthscale of the large scale motions (e.g.Taylor 1935;Oberlack 2002;Sagaut & Cambon 2018).Figure 15 shows the comparison between the theoretically predicted decay and that shown for the kinetic energy in the simulation at the apex of the tube using a value of D = 0.6, which is twice the value of the constant C 1 found in the first stage of the dynamics.Note that α has been multiplied by a factor of 4/3 because we predict that the total kinetic energy will be 4/3 greater than the kinetic energy in the x-direction.This can be inferred from the model of Hillier & Arregui (2019) due to the approximate equipartition between the energy found in the component of the flow varying about the mean layer velocity for a shear flow and the turbulent component of the flow.Overall the match shown in Figure 15 between the model and simulation results is good, though this only over a relatively short time range making stronger statements difficult.
If the energy of the turbulence is decaying, the reason for this is that the energy held in the turbulent fluctuations is being dissipated.Therefore, we can expect the decay of the turbulence to directly connect to heating.If the turbulent energy goes down, this should directly correspond to the same level of increase in thermal energy.To make progress, we assume the turbulent energy (taking the turbulent energy to be the sum of turbulent kinetic energy and turbulent magnetic energy) decay is uniform along the flux tube.Figure 16 shows the distributions of the turbulent kinetic and magnetic energies along the tube at two times.From this we can see that, roughly speaking, the magnitude of the total turbulent energy is approximately constant along the tube at both times implying an approximately constant decay along the tube.Making this assumption means we can integrate E(t 0 ) − E turb over the length of the tube and this gives the amount of energy dissipated.To take into account the velocity (and with that magnetic field) fluctuations in the y-direction we multiply E(t 0 ) − E turb by 4/3 and again take the same factor of 4/3 in α.The comparison between the relative increase in internal energy of the model and the simulation are shown in Figure 17 with the black line giving the simulation results and the red line that of the model (the model is only shown after a time of t = 42 for which it was derived).The model clearly provides a reasonable representation of the evolution implying that the key physical concepts that are needed to understand the increase in thermal energy of the simulation have been included in the model.

SUMMARY AND DISCUSSION
In this paper we have investigated the dynamics of a nonlinear kink wave excited in a flux tube.The numerical study was designed to capture the fundamentals of the system we are interested, i.e. kink oscillations of coronal loops, whilst maintaining sufficient simplicity and conserving energy in the domain to make analysis straight-forward.The oscillating tube develops the Kelvin-Helmholtz instability on its surface, developing a turbulent layer.It is the growth of this turbulent layer that characterises the first stage of the dynamics in the simulation.Once the layer becomes sufficiently large, and it became the dominant source of kinetic energy, the dynamics transition to a second stage characterised by the decay of the turbulent energy in an annulus around the much reduced core of the loop.We presented simplified analytic models to explain key aspects of these two stages.
The first stage of the dynamics is characterised by the growth of a turbulent layer that grows until most of the tube cross-section is turbulent.To model this evolution, we combined together the mixing models of Hillier & Arregui (2019) and Hillier et al. (2023) to create new models for how the mass and momentum (and as a consequence energy, wave velocity, wave amplitude and heating) evolve over time.This model made some clear predictions about when this stage has to end.As this transition is linearly proportional to the magnitude of the shear flow, we would naturally expect that this stage finishes sooner for highly nonlinear waves.
A fundamental aspect of the model for this stage of the evolution is the self-similar evolution of the turbulent layer.This model was developed for steady hydrodynamic shear flows, e.g.Winant & Browand (1974).
However, in this study we have a strong magnetic field, an oscillating flow and a curved boundary.In spite of these departures from the base model, the self-similar evolution can still be clearly distinguished.This implies that this is a robust physical mechanism that could occur in flows in many systems with strong magnetic fields.As the self-similar model comes from dimensional analysis of the system, it inherently has unknown constants associated with it.For the case presented in this paper, the constant we determined was consistent with previous hydrodynamic studies of shear layer mixing.However, a wider parameter survey is necessary to determine if this constant has some as yet unknown dependency on the system that is as yet not built into the model.
In the second stage of the evolution, as the energy of the turbulent motions dominates the energy of the coherent wave motions the model we propose to understand this phase focuses on the decay of this turbulent energy.This model provides a good prediction of the temporal decay of the turbulent energy of the simulation.This allows for a reasonable prediction of the heating in the simulation due to this turbulent decay to be estimated.One key aspect of this model is we assume that the characteristic length-scale of the turbulence remains fixed.However, the growth of the turbulent area shown by the difference between the blue and red lines in Figure 2 implies that even after t ≈ 50 the area does increase.This could imply that there is also some growth in the largest scale of the turbulence in the layer over time.This would change the exponent of the power law of the turbulent decay (e.g.Kolmogorov 1941;Oberlack 2002).Clearly, the simulated energy evolves at close to ∼ t −2 , but it would be interesting to explore how including an evolving length-scale improves the modelling of the system.
With both phases of dynamics, what is striking about the models proposed is that they both are taken almost directly from their hydrodynamic equivalents without major revision due to the presence of strong magnetic fields.We have implicitly used the magnetic fields in both phases, mainly through assuming that Alfvén waves are travelling back and forth along the magnetic field to connect turbulent fluctuations in one region in z to other regions in z.In the first phase this implies the influence of the mixing at the apex is felt the full length of the tube, creating velocity and magnetic field fluctuations all along the tube, and allowing the coherent oscillations of the mixing layer to develop.For the second phase, this connection means that we assume turbulence has developed along the full length of the tube and assuming the magnitude of the turbulent energy (whether it is in turbulent kinetic energy or turbulent magnetic energy) is the same along the tube due to the magnetic field connecting regions along z.This assumption allows us to calculate heating rates for the system.
Another consequence of the magnetic field is it makes the turbulent motions highly anisotropic (with KHi swirls being highly elongated along the z-direction (e.g.Antolin et al. 2018) giving quasi-2D like behaviour.This may imply that the value of the constant C 1 should be nearer to 0.5 than the value of ∼ 0.3 we have found for the first phase of the dynamics.However, we know that some of the energy of the turbulence is held in magnetic fluctuations, which will reduce the magnitude of the turbulent motions, slowing mixing.Assuming that there is equipartition between the energy of velocity and magnetic field fluctuations this would correspond to a factor of 1/ √ 2 in Equation 11 which would increase the value found for C 1 by √ 2. Taking that the velocity shear is maximal at the apex and zero at z = 0 because of the nature of the MHD kink wave under study, this would again put another factor of 1/ √ 2 in Equation 11 and another subsequent increase in the value found for C 1 by √ 2. Therefore, in some sense part of the effect of the magnetic field on the turbulence model is wrapped up in the value of the constant C 1 .

Resonant absorption and its lack of consequence in our simulations
For the simulations presented in this paper, as with many others of KHi turbulence developing from oscillations of tubes with thin boundaries (Terradas et al. 2008;Magyar & Van Doorsselaere 2016), the turbulent mixing broadens the tube boundary.This takes the tube density profile to be one where resonant layers (Goossens et al. 1992) would exist for a linear wave perturbation.Therefore, it may be expected that resonant absorption could develop in the tube once mixing has made a clear boundary layer between the inside and outside of the tube.However, we find no evidence of a resonant layer forming in the mixing layer at any point in our simulation and to explain the evolution of the turbulence in our simulation our model does not need to invoke resonant absorption at all.Beyond the visual lack of a resonant layer, both panels of Figure 11 show that the momentum and energy of the core of the tube are accurately accounted for by the model we propose.This can only be the case if the only change in the momentum and kinetic energy of the core is because the turbulent layer is eroding the core.If there where significant energy transfer from the core to the external layer through resonant absorption, these predictions would overestimate the momentum and kinetic energy of the core.Though we cannot rule out local field lines in the mixing layer resonantly absorbing some energy from the core of the tube, it means that if it is there this process is dynamically unimportant.
So why is resonant absorption not present (or if it is, dynamically unimportant)?The simple answer is that the development of turbulence naturally inhibits the formation of a resonant layer.There are two explanations as to why we should not expect it to develop.Firstly, the KHi turbulence is constantly driving changes in the local conditions of the mixing layer meaning that the resonant field lines are constantly changing and through the swirling motion of the turbulence field lines with different Alfvén speed are locked together and unable to freely move (which inhibits them developing their own oscillations).Secondly (or another way of looking at this) is that turbulence at small scales can act to have a diffusion effect at larger scales (e.g.Taylor 1922).We can infer that there would be a strongly enhanced magnetic diffusion felt by the tube in the mixing layer as a consequence of turbulent motions (c.f.turbulent magnetic diffusion in mean-field dynamo models Brandenburg et al. 2023).As the efficacy of resonant absorption is greatly reduced in the presence of strong diffusion, it can be understood that turbulence would work to suppress the development of a resonant layer.It is worth noting that the magnitude of the turbulent diffusion scales with the magnitude of the velocity shear.Therefore, smaller shear flows will produce less turbulent diffusion and at sufficiently small magnitude may allow the resonant absorption process to manifest.
To understand which will be dominate, resonant absorption or the Kelvin-Helmholtz instability, we can consider the timescales of the system.In the initial evolution, if the KHi timescale is much shorter than that of resonant absorption, e.g. with thin boundary layers or large velocity amplitudes of the waves, then the KHi will likely grow and at no time will resonant absorption be important.However, in the converse situation, e.g. with large boundary layers or small amplitudes of the waves, then resonant absorption will dominate leading to the development of a resonant layer.However, this resonant layer is KHi unstable (e.g.Terradas et al. 2008;Antolin et al. 2015) which will drive turbulence that can kill the resonant absorption process.
There is, however, one further point to consider.The oscillation of the mixing layer (where a resonant layer could develop) appears to no longer be that of the original kink oscillation.This is its own oscillation being forced by momentum injection into the layer by mixing.Its frequency is not the kink frequency of the system but a hybrid of that frequency and the kink frequency of the layer itself.Therefore it is likely any resonant process would be associated with this new oscillation of the mixing layer and not that of the initial kink oscillation.If this is the case, then the growth of the KHi would completely remove the possibility of resonant absorption of the global kink mode from the system, only allowing energy to be absorbed from large-scale motions of the mixing layer to local resonances, which would preclude resonant absorption as a damping mechanism once the KHi mixing layer has developed.5.2.Why is this model different from the other "KHi" damping model The model presented by Van Doorsselaere et al. ( 2021) is labelled in their work as damping by the Kelvin-Helmholtz instability.However, the model presented in this paper and their model have vastly different formulations, with very different dependency on the density contrast (and in the case of the model in this paper two distinct regimes).This leads to the question: How can two different formulations be modelling the same physical phenomenon?The model presented in this paper is developed from direct analysis of the turbulence created by the Kelvin-Helmholtz instability and, as shown by the comparison in Hillier et al. (2020) and Section 2, gives results that are consistent with those found in simulations of loop oscillations that develop the Kelvin-Helmholtz instability.
The model presented in Van Doorsselaere et al. (2021) used the linear eigenfunction for a kink wave and then prescribed a nonlinear amplitude to the oscillation.By integrating the wave oscillation over one period they could calculate the energy extraction rate through these nonlinearities.As they only study the particular wave mode associated with the kink wave, there are no other modes that are excited in the system.This means that there is no perturbation to seed the linear growth of the Kelvin-Helmholtz instability, and as a consequence there can be no Kelvin-Helmholtz turbulence.That is to say, the model they calculated does not allow the growth of the Kelvin-Helmholtz instability so cannot have any damping through Kelvin-Helmholtz turbulence.
This leads to the question, what is the nature of the nonlinear damping investigated in Van Doorsselaere et al. (2021).The initial conditions of their model are not an exact solution of the nonlinear equations (unlike a linear Alfvén wave in incompressible nonlinear MHD).Therefore, the kink wave (characterised by a wavenumber k 0 and azimuthal wavemode m = 1) does work on the neighbouring Fourier mode (in this case with wavenumber k = 2k 0 and azimuthal wave number m = 2).As with any forcing problem, the energy will be trapped in that wave mode efficiently if the forc-ing (which occurs at twice the kink frequency), and the wave mode being forced are resonant (this is the basis of any wave turbulence argument, e.g.Goldreich & Sridhar 1995).This clearly is a different mechanism to KHi damping of a kink wave as presented in this paper, highlighting that when considering nonlinear mechanisms for wave damping there are many different ways these nonlinearities can manifest.

Heating rates and the implications for coronal heating
One of the important predictions of the models present here are the heating rates from both stages of the evolution.To make this comparison, we need to compare the characteristic heating rate of the simulation with the cooling time for comparable coronal plasma.The cooling time for the solar corona (based on a number density of 10 9 cm −3 ) is τ cool ∼ 10 3 s.In the time units of our simulation (i.e.sound crossing times), if we assume a loop length of 10 10 cm and a sound speed of 1.2 × 10 7 cm s −1 the dimensional cooling time becomes 24τ .As this is an exponential decay, it is the equivalent of the thermal energy reducing by a factor of 1/ exp(1).Note this ignores the losses through thermal conduction to the chromosphere and so should be taken merely as a crudely representative value.
For the heating, we can see from Figure 13, at t ≈ 40 the internal energy had increased by a value of ∼ 0.018.This gives a rate of increase of internal energy of 0.00045.By dividing the total internal energy at the start of the simulation divided by exp(1) by this number we get a heating time of ∼ 10 4 τ , which is significantly larger than the cooling time of 24τ .Therefore, even though the kink wave in this simulation is a relatively nonlinear perturbation (with the initial kick of 20% of the sound speed), the heating rates found are significantly smaller than those needed to balance radiative losses.The small heating rates we find in our simulation are consistent with those found by Howson et al. (2017) where explicit viscosity and magnetic resistivity were included.

A corollary on coronal seismology
The work presented in this paper leads to a very important corollary about the damping of kink waves in the solar atmosphere.It is standard to fit either exponential or, in some cases, Gaussian damping envelopes to the observed decay of loop oscillations (e.g.Nechaeva et al. 2019).However, these damping envelopes are fundamentally connected to linear wave theory.
The work presented in Sections 3.2 and 3.3 highlights how KHi turbulence produces non-linear wave damping.The wave damping envelope for this nonlinear mechanism has a different functional form to the linear damping profiles.Even more importantly, exactly how the tube motions are measured changes the damping envelope that is measured.This implies that when measuring damping of kink waves from observations of the solar corona, it is important to consider more than just linear damping envelopes to understand the evolution.The active branch of the (PIP) code is available on GitHub (https://github.com/AstroSnow/PIP).Data required to reproduce the figures is freely available to download from Zenodo doi:10.5281/zenodo.10655009.All simulation data is available upon reasonable request.

Figure 1 .Figure 2 .
Figure 1.Contour plots of the density distribution in the (x, y)-plane taken at z= Lz showing the time evolution of the KHi on the surface of the oscillating flux tube.Time is given in units of the sound crossing time.Colour contours show densities between a lower value of 0.99 and an upper value of 3.01.

Figure 3 .
Figure 3. Contour plot of the x velocity in the x-t plane at y = 0 and z = Lz.blue (red) colours show positive (negative) values.Colour contour is only given for regions with ρ ≥ 1.62.
Figure 6 panel (a) shows the kinetic energy distribution at the loop apex (i.e.z = 10) and panel (b) shows the magnetic energy at the loop foot point (z = 0).

Figure 4
Figure 4. x-component of the velocity in the x-y plane at z = Lz.Colours correspond to same velocities as used in Figure 3. Black lines show the ρ = 1.1 and ρ = 2.7 transition.snapshots are taken at the same time as those shown in figure 1

Figure 5 .
Figure 5. Plot of the evolution of the amplitude of the displacement of the centre-of-mass of the tube in the x direction over time.The three curves are shown for the centre-of-mass calculated for density thresholds of ρ ≥ 1.1 (red), ρ ≥ 1.62 (black) and ρ ≥ 2.7 (blue).

Figure 6 .
Figure 6.Kinetic energy at the loop apex (panel a) and magnetic energy at the foot point (panel b) at t = 35.2.White lines show the density contour for ρ = 1.1 and ρ = 2.7

Figure 7 .
Figure7.Variation of vx with y for t = 11.6 (blue line) and t = 18 (red line).The sign of vx is set such that it is positive at y = 0 for easy comparison.

Figure 8 .
Figure 8. Schematic to explain the mixing model.Panel (a) shows the cross-section of the tube with radius R, the purple arrows show the direction of the oscillations of the tube.Panel (b) shows the geometric manipulation of the tube cross-section used to formulate the analytic model with the width of the now rectangular cross-section kept at 2R and the height set as 2H with the value of H set to conserve the cross-sectional area.Panel (c) shows how the development of two Kelvin-Helmholtz mixing layers that are growing over time and each have height h (delimited by the red arrows).

Figure 11 .
Figure11.Comparison of the model and the simulation for (a) the momentum evolution and (b) the kinetic energy associated with the bulk motions of the oscillating tube.The black lines show the distribution for material with ρ > 2.7 for the simulation (solid) and model (dashed).The red lines show the evolution for the quantities where 1.1 < ρ < 2.7 for the simulation (solid) and model (dashed).

Figure 12 .
Figure 12.Comparison of the model and the simulation for the centre of mass velocity amplitude (panels (a) and (c)) and the centre-of-mass amplitude of the oscillatory motions (panels (b) and (d)).Panels (a) and (b) show the temporal evolution for the material with ρ ≥ 2.7 and panels (c) and (d) show the temporal evolution for the material with ρ ≥ 1.1.The solid lines are calculated from the simulation and the dashed lines are V CoM,ρ≥ρ T (panels (a) and (c)) or A ρ≥ρ T (panels (b) and (d)).

Figure 13 .
Figure13.Change in internal energy (solid black line) with predicted slope based on the model of turbulent heating as the Reynolds number tends to infinity (solid red line) and the predicted upper limit of the heating rate (dashed red line).

Figure 15 .
Figure 15.Plot of the normalised turbulent energy against modified time.The black line shows the results of the simuand the red line the predicted decay from the model.

Figure 16 .
Figure16.Plot of the normalised turbulent energy in the x and y components of the magnetic field (red line) and velocity field (blue line) and their total (black line) against z calculated in a tube aligned with z with radius 0.45.The solid lines are for t = 60 and the dashed lines for t = 80.

Figure 17 .
Figure17.Evolution of total thermal energy over time (black line).Model of the expected increase in thermal energy as a consequence of heating due to decay of turbulence is shown with the red line.Note that the red line is only plotted for t ≥ 42.
AH is supported by STFC Research Grant No. ST/V000659/1.IA and AH are supported by project PID2021-127487NB-I00 from Ministerio de Ciencia e Innovación and FEDER funds.AH and TM are supported by JSPS KAKENHI Grant Number 19K03669.AH would like to acknowledge the discussions with members of ISSI Team 457 "The Role of Partial Ionization in the Formation, Dynamics and Stability of Solar Prominences", which have helped improve the ideas in this manuscript.AH and TM would like to acknowledge support by the Research Institute for Mathematical Sciences, an International Joint Usage/Research Center located in Kyoto University.For the purpose of open access, the author has applied a 'Creative Commons Attribution (CC BY) licence to any Author Accepted Manuscript version arising.