Resonant Damping of Propagating Kink Waves in Non-stationary, Longitudinally Stratified, and Expanding Solar Waveguides

Propagating kink waves have been observed in many magnetic waveguides in the solar atmosphere, like coronal magnetic loops, spicules, and fine structures of prominences. There are also observational evidences that these waves are damped. At present resonant absorption is considered as the most likely candidate for explaining this damping. First the attenuation of propagating kink waves due to resonant absorption was studied using the simplest model with a straight magnetic tube and the density only varying in the radial direction. Later a more advanced model with the density also varying along the tube was used. It was shown that the variation of the wave amplitude along the tube is determined by the combined effect of resonant damping and the longitudinal density variation. In our article we extend the analysis of resonant damping of propagating kink waves to take into account the magnetic loop expansion. We also consider non-stationary magnetic tubes to model, for example, cooling coronal loops.


INTRODUCTION
Propagating kink waves have been observed in many magnetic waveguides in the solar atmosphere, like coronal magnetic loops (Tomczyk et al., 2007;Tomczyk and McIntosh, 2009), spicules (De Pontieu et al., 2007;He J.-S. et al., 2009), fine structures of prominences (Okamoto et al., 2007), and in filament threads (Lin et al., 2007(Lin et al., , 2009. It was also observed that these waves damp. At present resonant absorption is considered as the most likely candidate for explaining this damping. The theoretical modeling of the spatial damping of traveling kink waves due to resonant absorption has been carried out by Terradas et al. (2010) and Verth et al. (2010) analytically, and by Pascoe et al. (2010Pascoe et al. ( , 2011 numerically. In all these article the simplest model of a straight magnetic tube with the density only varying in the radial direction was used. Later mode sophisticate models were used. Soler et al. (2011a) took into account the effect of partial ionization in the single-fluid approximation. As a result, kink waves were damped by both resonant absorption and ion-neutral collisions. Soler et al. (2011b) studied the resonant absorption of propagating kink waves in the presence of flow. Soler et al. (2011c) investigated the propagation and resonant absorption of kink waves in a magnetic tube with the density varying both along and across the tube. They showed that the variation of the wave amplitude along the tube is determined by the combined effect of resonant damping and the longitudinal density variation. Ruderman et al. (2010) studied the effect of non-linearity on the resonant damping of propagating kink waves and showed that non-linearity can strongly enhance the damping efficiency.
All papers cited above used the theory of resonant damping that can be called classical. In this theory the wave amplitude decays exponentially with the distance from the place where it is driven. This result is based on the assumption that a propagating kink wave is an eigenmode of the linearized dissipative magnetohydrodynamics (MHD). In the case of standing kink oscillations (Ruderman and Roberts, 2002) showed that after the initial perturbation the kink oscillation of a perturbed magnetic tube is very well-described by an eigenmode of the linear dissipative MHD after a transitional time of order of the oscillation period everywhere but in a vicinity of the resonant surface. However, in the vicinity of the resonant surface phase mixing continues until it creates perturbations with so small spatial scale that viscosity and/or resistivity stops it. Only after that the perturbation is described by an eigenmode of the linear dissipative MHD everywhere. It was shown numerically that, for typical parameters of coronal magnetic loops, the time when the phase mixing stops is at least by an order of magnitude larger than the typical damping time of kink oscillations (Arregui, 2015).
Hence, we conclude that the main assumption of the classical theory of resonant damping is not satisfied. This problem was addressed numerically by Pascoe et al. (2012Pascoe et al. ( , 2013, and analytically by Hood et al. (2013). It was shown that at the initial stage the amplitude variation with the distance from the driver is described by the Gaussian profile. And only later the amplitude decays exponentially. As a result, the damping length of a kink wave is somewhat longer than that predicted by the classical theory of resonant damping. Hence, to correctly describe the spatial damping of propagating kink waves we need to use the advanced theory developed by Hood et al. (2013). The importance of Gaussian damping strongly depends on the thickness of the transitional layer where the density drops from large value inside the tube to low value in the surrounding plasma. It also depends on the ratio of densities inside and outside the tube. The distance where the transition from the Gaussian to exponential damping occurs reduces fast when the thickness of the layer decreases, and also when the ratio of densities increases. Ruderman and Terradas (2013) carried out an analytical analysis of resonant damping of standing kink waves similar to that made by Hood et al. (2013) for propagating kink waves. They, in particular, concluded that the classical theory of resonant damping underestimates the damping time. But, for typical values of coronal magnetic loop parameters, the error never exceeds 20%. Although a similar estimate was not obtained for propagating kink waves, on the basis of the analogy between the spatial damping of propagating waves and temporal damping of standing waves, we believe that, although the classical theory of resonant damping underestimates the damping length, the error is quite moderate. On the other hand, the advanced theory of resonant damping is much more complex than the classical theory. Hence, in this article we will use the classical theory of resonant damping.
Sometime it is observed that waveguides in the solar atmosphere are non-stationary. For example, Aschwanden and Terradas (2008) and Aschwanden and Schrijver (2011) reported observations of kink oscillations of cooling coronal loops. Inspired by these observations Ruderman (2011b) and  studied the resonant damping of kink oscillations of cooling coronal magnetic loops. Morton et al. (2010) and Ballester et al. (2018) investigated the propagation of magnetosonic waves in a cooling plasma. In our paper we extend their analysis to propagating kink waves in non-stationary magnetic flux tubes. In particular, we study the kink wave propagation in cooling coronal loops.
Our article is organized as follows. In the next section we formulate the problem and present the governing equations. In section 3 we derive the equation governing the evolution of kink waves propagating along non-stationary magnetic tubes. In section 4 we derive the expression determining resonant damping of kink waves. In section 5 we derive the equation for the wave amplitude. In section 6 we consider kink wave propagation in static and non-expanding magnetic tubes. In section 7 we study kink wave propagation in non-stationary and expanding magnetic tubes. Finally, in section 8 we summarize the results obtained in the paper.

PROBLEM FORMULATION AND GOVERNING EQUATIONS
We study propagating kink waves along a straight magnetic tube with the circular cross-section of variable radius (see Figure 1). The characteristic radius of the tube cross-section is R * . Since the tube expands the magnetic field is spatially dependent, but the scale of its spatial variation is L * ≫ R * . This assumption implies that we consider a thin magnetic tube. Below we use cylindrical coordinates r, ϕ, z with the z-axis coinciding with the tube axis. We consider an axisymmetric equilibrium meaning that all equilibrium quantities are independent of ϕ. The magnetic field is not twisted meaning that its azimuthal component is zero. Its radial and axial components are expressed in terms of the magnetic flux function ψ as Since h ′ (z) ∼ h(z)/L * , it follows that where ε = R * /L * and B = (B 2 r + B 2 z ) 1/2 is the magnetic field magnitude.
The plasma density varies both along and across the tube. The tube consists of a core region and a boundary layer where the density monotonically decreases from its value ρ i inside the tube to its value ρ e in the surrounding plasma. Here and below the indices "i" and "e" indicate that a quantity is calculated in the core region and in the surrounding plasma, respectively. The characteristic thickness of the boundary layer is lR * , where l ≪ 1. This implies that we use the thin boundary approximation. The transitional layer boundaries are defined by equations r = R(z)(1 − l/2), r = R(z)(1 + l/2). (2.5) Now it follows from the magnetic flux conservation that the magnetic field tube magnitude and the tube radius are related by the approximate equation It follows from Equations (2.2) to (2.6) that the transitional layer boundaries are magnetic field lines and their equations can be written in an alternative form as (2.7) In the core region and in the surrounding plasma the characteristic scale of density variation is L * both in the longitudinal and radial direction. In the transitional layer the characteristic scale of density variation in the longitudinal direction is also L * , but it is lR * in the radial direction. The density also can depend on time. The temporal variation of density can cause the plasma flow along the magnetic field lines. Again, the velocity only weakly varies in the radial direction in the core region and outside of the tube, but it varies on the scale lR * in the transitional layer. The density ρ and velocity U = (U r , 0, U z ) are related by the mass conservation equation Since the velocity is parallel to the magnetic field it follows from Equation (2.4) that where U = (U 2 r + U 2 z ) 1/2 is the velocity magnitude. We integrate Equation (2.8) over the area of the tube core cross-section, that is over a circle of radius R(z)(1 − l/2). As a result we obtain (2.11) Substituting this result in Equation (2.10) we obtain in the leading order approximation with respect to ε and l Next we integrate (Equation 2.8) over the ring region R(1+l/2) ≤ r ≤ ςR, where ς − 1 ≫ l and ς is of the order of unity. This yields (2.13) Similar to Equation (2.9) we obtain U r U r=ς R = ςR ′ . (2.14) Using Equation (2.11) with −l substituted for l and (2.14) we obtain from Equation (2.13) In the leading order approximation with respect to ε and l we neglect l in comparison with unity, and take U z ≈ U and in Equation (2.13). Then, dividing the obtained equation by R 2 (ς 2 − 1) yields It was shown in Paper I that long linear kink waves, which are waves with the wavelength much larger than R * , in an expanding and non-stationary magnetic tube is described by the equation (2.17) In this equations µ 0 is the magnetic permeability of free space, P is the perturbation of magnetic pressure, and ξ ⊥ is the plasma displacement in the ϕ = const plane and perpendicular to the magnetic field lines. In the case of nonexpanding tube ξ ⊥ = ξ r , where ξ r is the radial component of the plasma displacement. In a slowly expanding tube ξ ⊥ = ξ r [1 + O(ε)]. In the thin tube approximation ξ ⊥i is independent of r. This property is same as that first obtained in the case of nonexpanding magnetic tubes (e.g., Ruderman and Erdélyi, 2009). The quantities δη and δP are the jumps across the transitional layer defined by (2.19) Equation (2.17) with the right-hand side defined by Equation (2.19) is used below to study the propagation of kink waves.
In the thin tube approximation the radial and azimuthal components of both the plasma displacement and magnetic field perturbation are independent of r inside the tube and proportional to r −2 outside the tube. The wave energy density is equal to the sum of the kinetic and magnetic energy density. The kinetic energy density is proportional to the sum of the squares of the radial and azimuthal components of the plasma displacement, and the magnetic energy density is proportional to the sum of the squares of the radial and azimuthal components of the magnetic field perturbation. Hence, the wave energy density is independent of r inside the tube and proportional to r −4 outside the tube. The energy behavior in the boundary layer strongly depends on the dissipative coefficients and can be either monotonic or oscillatory. At distances from the place where the wave is driven that are much smaller than the damping distance the wave energy density in the transitional layer is quite small, of the order of l/R * ≪ 1. However, at distances comparable with the damping distance almost all energy is concentrated in the transition layer due to resonant absorption.

DERIVATION OF THE EVOLUTIONARY EQUATION
In this section we consider propagation of kink waves along an expanding and non-stationary magnetic tube. Using Equations (2.12) and (2.15) we transform Equation (2.11) to Now we consider waves with the length much longer than R * , but much shorter than L * . We denote the ratio of the characteristic wavelength to L * as ǫ ≪ 1. The condition the wavelength is much larger than R * implies that ǫ ≫ ε = R * /L * . We also assume that the wave period is much shorter than the characteristic scale of the temporal density variation. To study the wave propagation we look for the solution to Equation (3.1) in the form where θ is real and S is complex (Bender and Orszag, 1978). The presence of transitional layer results in the resonant damping of waves. We will see below that the damping length is of the order l −1 times the wavelength. On the other hand, the effect of inhomogeneity manifests itself on a distance from the place where the wave is driven that is of the order of ǫ −1 times the wavelength. We would like to derive the equation for the wave amplitude that takes both effects into account in the same order approximation. In accordance with this we put l = ǫ. When l ≫ ǫ the effect of resonant absorption strongly dominates the effect of Frontiers in Astronomy and Space Sciences | www.frontiersin.org the axial density variation that can be neglected. In the opposite case where l ≪ ǫ the effect of the axial density variation strongly dominates the effect of resonant absorption. We will see below that L is of the order of ǫ −1 when l = ǫ. This estimate inspires us to introduce L = ǫL. The characteristic time of wave damping due to resonant absorption is ǫ −1 times the wave period. If the characteristic time of density variation is much larger than that time then its effect can be neglected. On the other hand, if the characteristic time of density variation is much smaller then the damping time then the effect of density variation will strongly dominate the wave damping. We aim to study the competition of the two effects. In accordance with this we assume that the ratio of the wave period to the characteristic time of density variation is of the order of ǫ. Substituting η = S exp(iǫ −1 θ ) in Equation (3.1) and collecting terms proportional to ǫ −2 in the obtained equation yields This approximation is usually called the approximation of geometrical optics. Next, we collect terms proportional to ǫ −1 . This results in (3.5) Multiplying this equation by SR 4 and using Equation (2.6) we obtain (3.6) We further transform this equation to (3.7) We now transform the terms on the left-hand side of this equation that are not full derivatives. Using Equations (2.12) and (2.15) we obtain and using Equation (3.3) we rewrite Equation (3.9) as Here E is proportional to the wave energy density per unit length along the magnetic tube. Below we assume that the temporal variation of the density is very slow. To be specific, we consider as an example kink waves in cooling coronal loops. One observation of kink oscillation of a cooling coronal loop was reported by Aschwanden and Schrijver (2011). In this event the period of fundamental mode was 395 s, and the loop length was 163 Mm. Hence, the phase speed of the kink wave was 893 km/s. The cooling time was 2050 s. Taking this time as the characteristic time in Equation (2.12), and the loop length as the characteristic length, we obtain the estimate U i ∼ 80 km/s. Aschwanden and Schrijver (2011) did not give any information about the temperature of the plasma surrounding the loop. Cooling mainly occurs due to radiation with the intensity proportional to the plasma density squared. The plasma in the loop is much denser than that surrounding the loop. Hence, even if the external plasma cools, its cooling time is much larger than that of the plasma in the loop, and we can expect that U e ≪ U i . Consequently, we conclude that in the event reported by Aschwanden and Schrijver (2011) the speed of the flow induced by cooling is much smaller than the phase speed.
On the basis of this example we introduce the definition that the temporal density variation is very slow if the flow speed induced by this variation is much smaller than the phase speed.
The assumption that the temporal density variation is very slow enables us to neglect the terms containing U i and U e in Equation (3.3). In addition, we only consider waves propagating in the positive z-direction. Then Equation (3.3) reduces to .
( 3.12) We also can neglect the terms proportional to U i and U e in the expression for E and write it in the approximate form as (3.13)

DERIVATION OF EXPRESSION FOR L
The assumption that the temporal variation of density is very slow enables us to use the linear equations of static MHD to describe the plasma motion in the transitional layer. However, we take the dependence of the density on time. To remove the singularity at the resonant surface we take the viscosity into account. Below we use the system of equations derived by . In this system ψ is used as an independent variable instead of r. Since Shukhobodskiy and Ruderman (2018) considered a static problem with the density independent on time they took the perturbation of all variables proportional to e −iωt . However, it is easy to restore the time dependence. It suffices to substitute ∂/∂t for −iω. As a result, we obtain In these equations ξ ϕ is the ϕ-component of the plasma displacement, w = Bξ ⊥ , and ν is the kinematic viscosity. We note that r is the function of ψ and z. These equations are valid both in the transitional layer as well as in the core region and external plasma where we can neglect the terms proportional to ν. The characteristic scale of variation of perturbations with respect to z is lε −1 R * ≫R * . The characteristic time of variation of perturbations is lε −1 R * /V * , where V * is the characteristic value of the phase speed. We can take V 2 * = B 2 * /(µ 0 ρ * ), where B * and ρ * are the characteristic values of the magnetic field and density, respectively. Using these estimates we obtain that P ∼ l −2 ε 2 V 2 * R −1 * ξ ϕ . Then the ratio of the left-hand side of Equation (4.1) to the second term in the brackets in this equation is of the order of l −2 ε 2 ≪ 1, which implies that the left-hand side of Equation (4.1) can be neglected. The ratio of the third term in the brackets on the right-hand side of Equation (4.1) to the fourth term is l −1 ε 2 ≪ 1, so the third term also can be neglected. As a result, Equation (4.1) reduces to Finally, since in the core region the dependence of both B and ξ ⊥ on ψ can be neglected, in this region we also can drop the first term in Equation (4.4).
In the WKB method all dependent variables must have the same functional form. In accordance with this, recalling that η = S exp(iǫ −1 θ ) and l = ǫ, we put (4.5) Using the relations w = Bξ ⊥ , η = ξ ⊥i /R and η = S exp(iǫ −1 θ ), and taking into account that the dependence of ξ ⊥ on ψ in the core region can be neglected we obtain the expression valid in the core region,ŵ Substituting Equations (4.5) and (4.6) in Equation (4.1) with the small terms neglected, and in Equation (4.3) with ν = 0, and collecting the leading terms with respect to l we obtain the following equations valid in the core region, Eliminatingξ ϕ from these equations yieldŝ Now we proceed to calculating L. First we further simplify Equations (4.2)-(4.4). We note that we can take r(ψ, z) ≈ R(z) in the transitional layer. Using Equations (2.3) and (2.4) we also take B z ≈ B. We can disregard the dependence of B on ψ. Finally, the characteristic scale of both the density and perturbations of all quantities in the transitional layer is ψ e − ψ i = lBR 2 . Using these estimates we can easily show that the second term in the square brackets in Equation (4.2) strongly dominates the first term, and the second terms in brackets in terms proportional to ν are much smaller than the first terms. Finally, the ratio of the third and first terms in Equation (4.4) is of the order of l. Hence, the third term can be dropped. Then, using Equation (2.6) we reduce the system of Equations (4.2)-(4.4) to (4.12)

Solution Outside of the Dissipative Layer
To obtain the solution in the transitional layer we use the method of matched asymptotic expansions (e.g., Bender and Orszag, 1978). In accordance with this method we split the transitional layer in the dissipative layer and two layers sandwiching this layer where we can neglect viscosity. We look for the solution to the linear dissipative MHD in the dissipative layer and to the linear ideal MHD outside of this layer. Then we match the two solutions in the two overlap layer where the both solutions are valid. The solution in the dissipative layer is called internal, and the solution outside of the dissipative layer external. We start from looking for the solution to the linear ideal MHD outside of the dissipative layer embracing the resonance surface defined by the equation ψ = ψ A , where ψ A is determined by the condition V A (ψ A ) = C k . Since the variation of P in the transitional layer is of the order of lP i , we can substitute P i (ψ = ψ i ) for P in Equation (4.12). Now, we substitute Equations (4.5) and (4.6) in Equations (4.10)-(4.12). Since we need to calculate L in the leading order approximation with respect to l, we only keep leading terms. As a result, we obtain where V A = B(µ 0 ρ) −1/2 is the Alfvén speed and ρ A = ρ(ψ = ψ A ). When deriving Equations (4.14) and (4.15) we used the relation ω = C k k. It follows from Equation (4.15) that . (4.16) We see that there is a singularity ofξ ϕ at ψ = ψ A (t, z). Substituting Equation (4.16) in Equation (4.13) and integrating the obtained equation yieldŝ (4.17) Whileξ ϕ has a singularity of the form (ψ − ψ A ) −1 ,ŵ only has a logarithmic singularity. Finally, substituting this result in Equation (4.14) and integrating the obtained equation results in It is easy to see thatP is continuous at ψ = ψ A .

Solution Inside the Dissipative Layer
We now look for the solution in the dissipative layer embracing the resonant surface. Ruderman et al. (1995) was the first to show that solution character in the dissipative layer depends on the value of viscosity (see also Ruderman and Roberts, 2002;Goossens et al., 2011). The spatial dependence of variable perturbations in the dissipative layer is monotonic when the viscosity is not very small, while it is oscillatory for very small values of viscosity. Ruderman et al. (1995) studied a planar problem where the transition from monotonic to oscillatory behavior is determined by the relative values of two small parameters, the ratio of the thickness of transitional layer to the wavelength, and the inverse Reynolds number. He also studied the temporal damping of kink waves. However, the results that he obtained is easily translated to the spatial damping and cylindrical geometry. In this case, the variable spatial dependence is determined by the relative values of three small parameters, l, ε = R * /L * , and the inverse Reynolds number Re −1 , where Re = R * V * /ν. The parameter determining the character of the spatial variation of variable perturbations in the dissipative layer is l(εRe) 1/3 . When l(εRe) 1/3 ≪ 1 the spatial dependence of variable perturbations in the dissipative layer is monotonic, while it is oscillatory when l(εRe) 1/3 1. We mainly aim to apply the results of this study to the solar atmosphere, where the typical value is l 0.1, ε 0.01, while Re ≫ 10 6 , so that l(εRe) 1/3 > 1, which implies that the behavior of perturbations in the dissipative layer is oscillatory. However, in this case the equations describing the motion in the dissipative layer are very complex and, at present, it is not clear how to solve them. On the other hand, we only need to calculate the jumps of w and P across the dissipative layer. Ruderman et al. (1995) found that these jumps are independent of the value of viscosity. The only condition that must be satisfied is that Re is sufficiently large, so that the thickness of the dissipative layer is much smaller than the thickness of the transitional layer. This result was later confirmed in subsequent studies (see, e.g., the review by Goossens et al., 2011). The solution of equations describing the plasma motion in the dissipative layer in the case when l(εRe) 1/3 ≪ 1 is relatively simple. The result that the jumps of w and P across the dissipative layer are independent of l(εRe) 1/3 was obtained for a static magnetic tube with a constant cross-section radius. However, it looks like a viable conjecture to assume that this result remains correct even for a non-stationary and expanding tube. We note that the following derivation is similar to that in the case of a non-expanding tube. The only difference that in the case of a non-expanding tube we use the variable r, while in our derivation we use the variable ψ instead.
Hence, we assume that l(εRe) 1/3 ≪ 1. Since the thickness of the dissipative layer is much smaller than the thickness of the transitional layer we can approximate any equilibrium quantity in the dissipative layer by its first non-zero term of Taylor expansion with respect to ψ −ψ A . In particular, we can substitute ρ A = ρ(ψ = ψ A ) for ρ and take (4.20) Since we assume that the density monotonically decreases in the radial direction in the transitional layer, it follows that > 0. Now, substituting Equation (4.5) in Equation (4.12), collecting terms of the order of l −2 , using Equations (4.9) and (4.20), and substituting P(ψ = ψ i ) for P we obtain (4.21) When deriving this equation we took into account that ω = C k k.
The thickness of the dissipative layer is defined by the condition that the two terms on the left-hand side of Equation (4.21) are of the same order. Using Equations (2.2) and (2.4) we easily obtain that this thickness is Then the condition that the thickness of the dissipative layer is much smaller than the thickness of the transitional layer reduces to εRe ≫ 1. For typical conditions in the solar atmosphere this inequality is definitely satisfied. Together with the condition that the spatial behavior of variable perturbations in the dissipative layer is non-oscillatory this gives 1 ≪ εRe ≪ l −3 . (4.23) The solution in the dissipative layer has to match the solution outside of this layer in the overlap layer defined by δ A ≪ |r − r A | ≪ lR * . Using Equation (4.16) we obtain that the solution in the overlap layer has the form Hence, the solution to Equation (4.21) must have this form for |ψ − ψ A | ≫ R * B * δ A . The solution to Equation (4.21) satisfying this condition is obtained in Appendix A. It is given by Equation (A6). Using Equation (A1) we rewrite it aŝ Using Equations (4.25) and (4.26) we obtain from Equation (4.13) The functions F( ) and G( ) were introduced by Goossens et al. (1995). Finally, substituting Equation (4.5) in Equation (4.11), collecting terms of the order of l −2 , and using the relation ω = C k k, and Equations (4.20) and (4.22), and (4.26)-(4.28) we obtain (4.29)

Matching Solutions
The matching procedure is the following. First we find the asymptotic expansion of the internal solution valid for ≫ 1. Next we find the expansion of the external solution valid for |ψ − ψ A | ≪ BR 2 . Then we substitute ψ − ψ A = RBδ A in this expansion. The matching condition is that the leading terms of the two expansions must coincide. We found that it is more convenient to compare not the expansions but the jumps across the dissipative layer. The jump of w across the dissipative layer is given by w( ) − w(− ) with ≫ 1. We obtain It is obvious that the second term in this expression tends to zero as → ∞. Using this result we obtain from Equation (4.27) that the jump ofŵ across the dissipative layer iŝ (4.31) Using Equation (4.17) we obtain another asymptotic expression for the jump ofŵ across the dissipative layer, (4.32) where δŵ =ŵ(ψ = ψ e ) −ŵ(ψ = ψ i ) and P indicates the principal Cauchy part of an integral. This asymptotic expression is valid for |ψ − ψ A | ≪ 1. The leading terms of the two asymptotic expressions, one given by Equation (4.31) and the other by Equation (4.32), must coincide. It follows from this condition that Using Equation (4.29) we obtain that the jump of P across the dissipative layer is given bŷ The integral on the right-hand side of this equation is evaluated in Appendix B. Using Equation (B8) we obtain This result and the matching condition imply that the expansion with respect to ψ − ψ A of the jump of P across the dissipative layer calculated using the external solution must start from the term proportional to (ψ − ψ A ) 2 . In particular, it follows from this condition that the term in this expansion proportional to unity must be zero. Using Equation (4.18) we write this condition as (4.36) where δˆ =ˆ (ψ = ψ e ) −ˆ (ψ = ψ i ). Since δŵ/ŵ = O(l) and we only need to calculate δˆ in the leading order approximation with respect to l, we can substituteŵ(ψ = ψ e ) ≈ŵ(ψ = ψ i ) = RBS. Then, noticing that the only quantity that depends on ψ in Equation (4.36) is V A , each single integral is of the order of l, and each double integral is of the order of l 2 , we reduce this equation in the leading order approximation with respect to l to We note that we would obtain exactly the same expression for δˆ if we assume from the very beginning that we can neglect the jump of the pressure perturbation across the dissipative layer. This assumption was first made ad hog by Hollweg and Yang (1988). Later it was rigorously proved in 1D plasma equilibrium by Goossens et al. (1995).
Since the jump of across the dissipative layer is zero, it follows that the expression for obtained using the ideal MHD equations is a continuous function in the whole transitional layer. Then, using Equation (4.37) we obtain from Equation (4.18) the expression valid in the whole transitional layer in the leading order approximation with respect to l, (4.38) When deriving this expression we neglected the second terms in the brackets in Equation (4.18) because their ratios to the first terms are of the order of l. Using Equation (4.19) we obtain (4.39) Now we proceed to calculating L. We substitute η = Se iǫ −1 θ , δη = (δŵ/RB)e iǫ −1 θ , and δP = δP e iǫ −1 θ in Equation (2.17). Then, using Equation (3.12) and the condition of very slow temporal density variation implying that U e ≪ C k we obtain (4.40) Finally, using Equations (4.33) and (4.37) we arrive at

DERIVATION OF GOVERNING EQUATION FOR THE WAVE AMPLITUDE
The wave evolution is described by Equation (3.11) with E and L given by Equations (3.13) and (4.41), respectively. We write S = Ae iχ . Then, substituting Equations (3.13) and (4.41) in Equation (3.11), multiplying the obtained equation by e −2iχ , and separating the real and imaginary parts yields where Q = (ρ i + ρ e )ωR 4 A 2 . Equation (5.1) determines the temporal and spatial dependence of the wave amplitude, while Equation (5.2) describes a small phase shift. We are mainly interested in the variation of the wave amplitude in space and time, so we do not use Equation (5.2) below.

WAVE PROPAGATION ALONG A STATIC AND NON-EXPANDING WAVEGUIDE
In this section we reproduce the results previously obtained for static and non-expanding waveguides. Hence, we assume that the tube radius is constant and equal to R.

Waveguide Homogeneous in the Axial Direction
Here we consider the same problem as that studied by Terradas et al. (2010), which is the resonant damping of kink waves propagating along a magnetic tube homogeneous in the axial direction. We now assume that the density only varies in the radial direction. We assume that a harmonic wave is driven at z = 0 and propagates in the region z > 0. In that case A is independent of time and it follows from Equation (5.1) that A = A 0 e −γz , where A 0 is the amplitude at z = 0. Using the relation ψ = 1 2 Br 2 we obtain from Equation (4.20) In this case both k and ω are constant. Hence, θ = kz − ωt, which implies that the wavenumber is k * = l −1 k. Now, using Equation (6.1) and the relation C k = V A (r A ) yields This expression coincides with that obtained by Terradas et al. (2010) (see their Equation (10) with m = 1).

Waveguide With the Density Varying in the Axial Direction
Now we study the resonant damping of kink waves propagating along a magnetic tube with the density varying along the tube. This problem was first addressed by Soler et al. (2011c). We aim to reproduce their results. We assume that ρ i (z)/ρ e (z) = ζ = const and ρ(r, z)/ρ e (z) = f (r). Previously these assumptions were made by Dymova and Ruderman (2006) when studying resonant damping of standing kink waves, and by Soler et al. (2011c) when studying resonant damping of propagating kink waves. We again assume that the kink wave with the amplitude A 0 and the constant frequency ω is driven at z = 0. Since now C k is a function of z, the same is true for the wavenumber: k(z) = ω/C k (z). Note that in non-scaled variables the wavenumber is k * (z) = l −1 k(z).
Since Q is again independent of t it immediately follows from Equation (5.1) that When deriving this equation we used the relation ρ i C 2 k = ρ f C 2 f . Equation (6.1) remains valid. Then, using the relation ρ(r, z) = f (r)ρ e (z) we obtain from Equation (4.42) (6.5) G = 4/π 2 for the linear density profile, and 2/π for the sinusoidal density profile. After substituting Equation (

WAVE PROPAGATION ALONG AN EXPANDING AND NON-STATIONARY WAVEGUIDE
As an example of application of the general theory we consider a generalization of the same problem that was studied by Soler et al. (2011c), and take the loop expansion and cooling into account. We first describe the general method for studying the wave propagation, and then apply it to a particular loop with given cross-section radius and density variation along the tube, and the temporal density variation.

General Theory
We assume that a kink wave is driven at one of the loop footpoints and impose the boundary condition ω = ω 0 , A = A 0 at z = 0. (7.1) Driving starts at t = 0. Before that the loop is at rest, so we also have the initial condition The equations describing the wave propagation are solved for t > 0 and z > 0. We start from calculating θ (t, z). It follows from dispersion Equation (3.12), ω = C k k, and Equation (3.4) that θ (t, z) satisfies the equation Since θ (t, z) is defined with the accuracy up to an additive constant we can take θ (0, 0) = 0. Then it follows from Equation (7.1) that Since the loop is at rest at t = 0 we can take θ = 0 at t = 0. (7.5) The equation of characteristics of Equation (7.3) is It follows from Equation (7.3) that θ = const along a characteristic. We consider the characteristic that starts at the coordinate origin. Let its equation be z = z b (t), where z b (0) = 0. This characteristic separates the perturbed and unperturbed regions in the (t, z)-plane, so we will call it the boundary characteristic. Let us consider a point with coordinates (t 1 , z 1 ) satisfying the condition z 1 > z b (t 1 ). This implies that this point is above the boundary characteristic. Since the characteristic containing the point (t 1 , z 1 ) cannot intersect the boundary characteristic, it follows that it starts at the z-axis. Then, using Equation (7.5) we obtain θ (t 1 , z 1 ) = 0. Now we consider a point (t 1 , z 1 ) that is below the boundary characteristic meaning that z 1 < z b (t 1 ). Let the characteristic containing this point starts at t = τ (t 1 , z 1 ) on the t-axis. Then θ (t 1 , z 1 ) = −ω 0 τ (t 1 , z 1 ). As a result, we determine θ (t, z) in the whole region t > 0, z > 0. Differentiating θ (t, z) with respect to t we calculate ω. Then k = ω/C k .
Next we proceed to solving Equation (5.1). The equation of characteristics of this equation is also Equation (7.6). The variation of Q along the characteristic is defined by After substituting in this equation a solution to Equation (7.6) found when calculating θ (t, z) , Equation (7.7) becomes the equation determining the variation of Q along a characteristic. The solution to this equation must satisfy the initial condition Q = (ρ i + ρ e )ω 0 R 4 A 2 0 at t = τ (t 1 , z 1 ). (7.8) In this equation the equilibrium quantities are calculated at t = τ and z = 0. We now consider a point (t 1 , z 1 ) with z 1 > z b (t 1 ), which implies that it is above the boundary characteristic. In that case the characteristic that contains this point starts at the z-axis where A = 0. Then it follows that A(t 1 , z 1 ) = 0, that is the tube is at rest for z > z b (t). Hence, the equation z = z b (t) describes the propagation of the wave front along the magnetic tube. Below we apply the general theory to particular cases.

Wave Propagation in Cooling and Expanding Coronal Loop
We now consider the kink wave propagation in a coronal loop of half-circle shape immersed in an isothermal atmosphere. We assume that the loop is in a vertical plane. Cooling of coronal plasma mainly occurs due to radiation. The radiation intensity is proportional to the plasma density squared. Since the plasma density inside the loop is substantially higher than that of the surrounding plasma, the plasma inside the loop cools much faster than that outside the loop. This observation inspires us to make a viable assumption that cooling only occurs inside the loop, while its temperature outside the loop does not change. Then the density inside and outside the loop is given by (7.9) where L is the length of the loop. Following to Aschwanden and Terradas (2008) and Ruderman (2011a,b) we assume that the plasma density inside the loop decreases exponentially, so that (7.10) Here we do not discuss the mechanisms of coronal loop cooling, although the main cause of cooling of moderately hot coronal loops with the temperature of the order of or less than 1.5 MK is the radiative cooling. Magyar et al. (2015) studied transverse oscillations of radiatively cooling coronal loops numerically. They did not present the dependence of temperature on time. However, in general their results related to the time dependence of the oscillation amplitude are in good agreement with that obtained by Ruderman (2011a,b) who assumed the exponential temperature decay. Hence, it seems that the exponential dependence of temperature on time is a reasonable approximation. The pressure inside the loop is assumed to be in equilibrium with the outside medium during the cooling. Since the plasma beta in the corona is very low this condition does not impose any serious restriction on the plasma parameters. We adopt a model of expanding coronal loop first introduced by Ruderman et al. (2008), and later also used by  and . In this model the cross-section radius of the magnetic tube is given by (7.11) Here l c is a free parameter with the dimension of length, and λ is the expansion factor equal to the ratio of the cross-section radius at the loop apex and footpoints, that is λ = R(L/2)/R f . In our numerical study we took L/l c = 6 and l = 0.2. Here it is worth making a comment. In section 2 we imposed the condition that L * ≫ R * , where L * is the characteristic scale of the tube radius variation along the tube, and R * is the characteristic tube radius. In the model that we adopted here L * = l c = L/6. Since typically L ∼ 50R * , it follows that L * /R * ∼ 8. Hence, the condition L * ≫ R * is satisfied. The condition that the speed of the flow caused by cooling is much smaller than the phase speed is N = L(t c C f ) −1 ≪ 1. (7.12) In our analysis we neglect the effect of the tube curvature and consider it as straight. To our knowledge the effect of tube curvature on the propagation and damping of kink waves has not been studied. However, it was studied in the case of standing waves. Van Doorsselaere et al. (2004) analytically and Terradas et al. (2006) numerically showed that the coronal loop curvature has very minor effect on the frequency and damping of kink oscillations. It looks like a viable assumption that the same is true for propagating waves when the curvature radius is much larger than the tube radius.
The main aim of our study is to investigate the effect of cooling on the kink wave propagation. Since cooling decreases ρ i it increases C k . Hence, the stronger the cooling is the faster the wave perturbation launched at one footpoint at t = 0 reaches the other footpoint. We studied the wave propagation for N = 0 (no cooling), N = 0.1 (slow cooling), N = 0.2 (moderate cooling), and N = 1/3 (strong cooling). In the case of strong cooling the wave front arrives at the second footpoint at t = t end . In all other cases the wave front arrives at the second footpoint at t > t end . We calculated the spatial dependence of the wave frequency, wavenumber, and the amplitude at t = t end . At z = 0 the wave frequency is l −1 ω 0 and the wavenumber is l −1 ω 0 /C f . In our calculation we took the wavelength at z = 0 equation to one fifth of L, that is L = 10πlC f /ω 0 = 2πC f /ω 0 . FIGURE 3 | Dependence of the frequency on the distance along the loop for T = T end . The left panels correspond to λ = 1 and the right to λ = 1.5. The upper, middle, and lower panels correspond to κ = 0.5, κ = 1, and κ = 2, respectively. The solid, dashed, dotted, and dash-dotted curves correspond to N = 1/3, N = 0.2, N = 0.1, and N = 0, respectively.
(7.15) The quantity T end is defined by the equation Z b (T end ) = ω 0 L(lC f ) −1 , where Z b (T) is the solution to Equation (7.15) with α = 1/30π satisfying the initial condition Z b (0) = 0. Using Equation (7.15) we calculated the dependence of the wave front position on time for various values of κ and λ. This dependence is shown in Figure 2. We see that the stronger the cooling is the faster the wave front moves. This results is not surprising because FIGURE 4 | Dependence of the wave number on the distance along the loop for t = t end . The left panels correspond to λ = 1 and the right to λ = 1.5. The upper, middle, and lower panels correspond to κ = 0.5, κ = 1, and κ = 2, respectively. The solid, dashed, dotted, and dash-dotted curves correspond to N = 1/3, N = 0.2, N = 0.1, and N = 0, respectively. the cooling causes the enhancement of the phase speed C k . We can also see that the effect of cooling is rather weak, although the tube expansion makes it more pronounced.
We also calculated the dimensionless frequency and wavenumber K for various values of λ, , and κ. The results are presented in Figures 3, 4. First of all, we note the frequency is constant when there is no cooling as it must be. When there is cooling, in a non-expanding loop the frequency increases with the distance from the footpoint where the wave is driven. The stronger the cooling is the more pronounced this effect is. However, this effect is quite weak. The dependence of frequency on the cooling rate is much stronger in an expanding loop. We see that it is especially strong when the tube expands and the loop height substantially exceeds the atmospheric scale height. The situation with the wavenumber is quite similar. Again cooling almost does affect it in non-expanding loops, while in expanding loops the effect of cooling is quite noticeable.
Finally, Figure 5 displays the variation of the amplitude along the loop. It is worth noticing that, in most cases, the amplitude first increases and then starts to decay. The amplitude increase is related with the stratification, while the decay is cased by the resonant damping. We can see that cooling always results FIGURE 5 | Dependence of the amplitude on the distance along the loop for t = t end . The left panels correspond to λ = 1 and the right to λ = 1.5. The upper, middle, and lower panels correspond to κ = 0.5, κ = 1, and κ = 2, respectively. The solid, dashed, dotted, and dash-dotted curves correspond to N = 1/3, N = 0.2, N = 0.1, and N = 0, respectively. in the amplification of waves. This result is similar to that found by Ruderman (2011a,b); Ruderman et al. (2017), and  in the case of standing kink waves.

SUMMARY AND CONCLUSIONS
In this article we studied the kink wave propagation along an expanding magnetic tube with the density varying along the tube and in time. The tube consists of the core region where the density is almost independent of the radial coordinate, and the boundary layer where the density decreasing fast from its value inside the core region to its value in the surrounding plasma. This value is assumed to be thin meaning that its thickness is much smaller than the tube radius. We used the cold plasma approximation. We also used the thin tube approximation meaning that the wave length is much larger than the tube radius, and the short wave approximation meaning that the wavelength is much smaller than the characteristic scale of the density and tube radius variation along the tube. Using the WKB method we derived the equation describing the dependence of the wave amplitude on time and on the distance along the tube.
First we studied the kink wave propagation in a magnetic tube homogeneous in the axial direction. In this case the only effect affecting the wave propagation is the wave damping due to resonant absorption. We reproduced the results previously obtained by Terradas et al. (2010).
We then proceeded to studying the kink wave propagation in a magnetic tube with the density varying in the axial direction. In this case the wave propagation are affected both by resonance absorption and the axial inhomogeneity. We reproduced the analysis by Soler et al. (2011c).
Finally, we studied the kink wave propagation in an expanding and non-stationary magnetic tube. We obtained the general expressions determining the spatial and temporal dependence of wave frequency, wavenumber, and amplitude. We then applied the general theory to a particular case of kink wave propagating along a cooling coronal loop. We assumed that the loop has a half-circle shape and immersed in an isothermal atmosphere, the temperature of plasma inside the loop decays exponentially, while the temperature of the surrounding plasma does not change. We adopted the dependence of the loop cross-section on the distance along the loop previously used by Ruderman et al. (2008Ruderman et al. ( , 2017 and . The equations governing the wave propagation were solved numerically. We assumed that the wave was started to be driven at one of the footpoints at the same time when the plasma inside the loop started to cool. Our main aim was to study the dependence of the wave properties on the intensity of cooling. First we studied the dependence of the distance that the wavefront travels on time. We found that the stronger the cooling is the larger the distance that the wave front travel at a given time. This is an expected result because cooling enhances the phase speed thus accelerating the wavefront. We also calculated the dependence of the wave frequency and wave number on the distance along the tube. When doing so we chose the moment of time when the wavefront arrives at the second footpoint in the case of strongest cooling. The general conclusion is that cooling results in the increase of the wave frequency. In contrast, it is difficult to make any definite conclusion about the effect of cooling on the wavenumber. Finally, we investigated the dependence of the wave amplitude on the distance along the tube. In most cases the amplitude first growths due to the equilibrium quantity variation along the tube, and then its starts to decay due to resonant damping. We found that cooling enhances the wave amplitude. This result is similar to one previously obtained for standing kink waves (Ruderman, 2011a,b;.

A. Solution to Equation (4.21)
In this section we obtain the solution to Equation (4.21) satisfying the asymptotic condition Equation (4.24). To simplify calculations we introduce the notation Using this notation we rewrite Equation (4.21) as To solve this equation we use the Fourier transform with respect to defined by (A3) Applying this transform to Equation (A2) we obtain where δ(σ ) is the delta-function. The solution to this equation decaying as |σ | → ∞ is given by where H(σ ) is the Heaviside step function. Calculating the inverse Fourier transform we obtain Using the integration by parts we obtain the asymptotic expression valid for | | ≫ 1, Using Equation (A1) it is straightforward to verify that Equation (A7) coincides with Equation (4.24).

B. Evaluation of Integral in Equation (4.34)
In this section we evaluate the integral on the right-hand side of Equation (4.34) for | ≫ 1. We immediately obtain