Overdense Threads in the Solar Corona Induced by Torsional Alfvén Waves

High-resolution and high-cadence observations have shown that Alfvén waves are ubiquitous in the solar atmosphere. Theoretical works suggest their ability to transfer large energy fluxes from the photosphere to the corona and solar wind. In this proof-of-concept Letter we show that torsional Alfvén waves can induce the formation of filamentary plasma structures in the solar corona. We perform high-resolution 3D ideal MHD simulations in an initially uniform coronal plasma permeated by a line-tied twisted magnetic field. We find that torsional Alfvén waves develop Kelvin–Helmholtz instabilities as a result of the phase mixing process. The Kelvin–Helmholtz instability drives plasma compression that breaks the uniformity of density, creating elongated overdense threads aligned with the direction of the magnetic field. With synthetic modeling of SDO/AIA imaging we show that the overdense filaments could be seen in observations as fine strands that illuminate the underlying magnetic structure.


Introduction
Alfvén waves are magnetohydrodynamic (MHD) waves whose restoring force is magnetic tension (Alfvén 1942). They are polarized transverselly to the magnetic field direction and are nearly incompressible. Torsional Alfvén waves (TAW) are a subclass of axisymmetric waves that drive periodic twisting (torsion) of the magnetic field lines (see, e.g., Spruit 1982). In the solar atmosphere, TAW have been detected in photospheric bright points (Jess et al. 2009), chromospheric spicules (De Pontieu et al. 2012), during solar flares (Aschwanden & Wang 2020), and in coronal structures (Kohutova et al. 2020). Theoretical works suggest that TAW are able to transport large energy fluxes from the photosphere to the corona and solar wind despite the strong chromospheric filtering (e.g., Soler et al. 2019Soler et al. , 2021. In the presence of a continuous spectrum, TAW undergo the process of phase mixing and develop shear flows that can trigger the Kelvin-Helmholtz instability (KHi; see Heyvaerts & Priest 1983;Browning & Priest 1984). In coronal loops, the nonlinear evolution of this process leads to turbulence (Díaz-Suárez & Soler 2021). Similar results were previously obtained in the case of kink MHD waves (see, e.g., Terradas et al. 2008;Antolin et al. 2014;Howson et al. 2017;Karampelas et al. 2017;Pagano et al. 2020, among other efforts). However, unlike kink MHD waves, TAW do not require a preexisting tube or density enhancement that acts as an established waveguide. Therefore, TAW can be considered as a more fundamental type of oscillation of the magnetic field.
This Letter is a proof-of-concept work whose purpose is to demonstrate that the nonlinear evolution of TAW can induce the spontaneous formation of filamentary structures in the solar coronal medium. The process relies on the excitation of the KHi via phase mixing shear flows and does not need the presence of a preexisting density enhancement. The KHi drives plasma compression in the form of elongated threads locally aligned with the direction of the magnetic field. To illustrate the fundamentals of this mechanism, we consider the paradigmatic scenario of a twisted magnetic field embedded in a uniform plasma. The evolution of standing TAW and the associated formation of overdense threads are studied with MHD numerical simulations, while synthetic modeling of extreme ultraviolet (EUV) emission is used to discuss possible observational signatures of the process.

Initial Configuration and Numerical Setup
In this numerical experiment, we aim to represent with a simple model the solar coronal plasma permeated by a twisted magnetic field. We assume a uniform background density, ρ = ρ 0 . Gravity is neglected. For the magnetic field, we consider the uniform twist model adapted from Gold & Hoyle (1960), namely with B r , B j , and B z the magnetic field components in a cylindrical coordinate system (r, j, z). B 0 is the magnetic field strength at r = 0, R is the radial length scale, and c is a parameter that controls the amount of twist. The magnetic field is invariant in the azimuthal, j, and longitudinal, z, directions. Since the magnetic field is force-free, we consider a uniform gas pressure, p = p 0 , whose value is set so that the plasma β at . Therefore, the simulations are conducted in a low-β scenario.
The Cartesian computational box extends from −3R to 3R in the x-and y-directions and from −L/2 to L/2 in the z-direction, with L = 10R. The magnetic field is assumed to be anchored (line-tied) at the bottom, z = −L/2, and top, z = L/2, boundaries, which represent the base of the corona. Hood & Priest (1979) studied the kink instability of the uniform twist model considering line-tying of the field lines and found that the maximum twist expressed with their parameter Φ = LB j /rB z is Φ crit. ≈ 3.3π. Their twist parameter Φ is related to our twist parameter by c = ΦR/L, so that the critical value is c crit. ≈ 1 in our model. Therefore, we shall consider c < 1.
At t = 0 we aim to excite the longitudinally fundamental mode of standing TAW, but we want to exclude other wave modes to be present initially. To this end, we consider an initial velocity field that is purely perpendicular to the magnetic field lines, since any perturbation parallel to the magnetic field would also excite slow MHD waves. We assume v 0, 4 is the modulus of the magnetic field and v ⊥ is the perpendicular velocity given by where v 0 is the maximum velocity amplitude and the radial profile, A r ( ), is the same as that used in Díaz-Suárez & Soler (2021). The velocity is zero at r = 0, its maximum is located at r ≈ 0.64R, and then it decreases as r increases. We is the Alfvén velocity at r = 0, with μ 0 the vacuum magnetic permeability.
We use the PLUTO code (Mignone et al. 2007(Mignone et al. , 2012 to numerically evolve the 3D nonlinear ideal MHD equations. The code setup is equivalent to that used in Díaz-Suárez & Soler (2021). The computational domain is initially discretized with a uniform mesh of 100 cells in each direction. We use adaptive mesh refinement (AMR) with five levels of refinement, so that the maximum effective resolution is 3200 cells in each direction. The refinement criterion is based on the second derivative error norm of a quantity related with the kinetic energy density. Regarding the boundary conditions, density, pressure, and the magnetic field components are fixed to their initial values while the velocity components are set to zero at all boundaries. This choice avoids unwanted problems like nonphysical energy injection from outside the domain, although reflections at the lateral boundaries may happen. However, the lateral boundaries are located sufficiently far away from the location where the important dynamics takes place and their influence on the simulations is minimal. Our aim is different from that of Díaz-Suárez & Soler (2021), where the onset of turbulence in a preexisting untwisted coronal loop was studied. Here, we investigate the formation of dense structures in an initially uniform corona with a twisted magnetic field.

Results
We study the evolution of density driven by the presence of TAW in the initially uniform plasma. Unless otherwise stated, we consider c = 0.6 in the simulations. Figure 1 shows the temporal evolution of density in a transverse plane (top panels) and in a longitudinal plane (mid panels). To focus on the relevant dynamics, Figure 1 only shows a smaller subdomain of the complete numerical box.

Formation of Overdense Threads
The plasma is initially uniform but, as time increases, the density evolution undergoes two differentiated phases. The first phase is caused by phase mixing of the TAW and the associated ponderomotive force. The second phase is driven by the KHi triggered when the shear flows generated by phase mixing are strong enough.
Besides the density, we also have explored the behavior of the perpendicular velocity component (not shown here). Before the KHi onset, v ⊥ shows the alternation of positive and negative velocities in adjacent radial positions, i.e., shear flows. This is an unequivocal signature of phase mixing (see, e.g., Heyvaerts & Priest 1983;Nocera et al. 1984;De Moortel et al. 2000;Smith et al. 2007;Díaz-Suárez & Soler 2021). The reason behind the presence of phase mixing in our model with uniform density is the magnetic field nonuniformity. According to Halberstadt & Goedbloed (1993), the line-tied Alfvén continuous spectrum is given by where n denotes the harmonic number. Even if the density is uniform, there is an Alfvén continuum owing to the radial variation of B z (see Equation (3)). The Alfvén continuum naturally leads to the development of phase mixing.
Density perturbations are produced during the evolution of the phase-mixed TAW. These perturbations appear in the form of concentric rings in the transverse cut (top left panel of Figure 1) and as periodic converging/diverging flows in the longitudinal cut (mid left panel of Figure 1). These periodic density perturbations are caused by the ponderomotive force, which is a well-known nonlinear effect that couples Alfvén waves with slow MHD waves (e.g., Hollweg 1971;Rankin et al. 1994;Tikhonchuk et al. 1995;Terradas & Ofman 2004). For phase-mixed TAW, the longitudinal flows by the ponderomotive force also get out of phase as time increases (Díaz-Suárez & Soler 2021). The development of density perturbations intensifies with the triggering of the KHi owing to the shear flows associated with phase mixing (Heyvaerts & Priest 1983;Browning & Priest 1984). Unlike in Díaz-Suárez & Soler (2021), where the KHi rolls are 2D structures formed in planes perpendicular to the straight magnetic field, here the KHi rolls are 3D structures because of twist (see, e.g., Howson et al. 2017;Terradas et al. 2018).
The KHi evolves nonlinearly and develops secondary instabilities that further mix and compress the plasma (see the bottom panel of Figure 1). During the evolution of the simulation, we found overdensities of ∼10% the initial background density. To study the 3D shape of the overdense structures, Figure 2 shows isosurfaces of the density for a value corresponding to ∼5% enhancement. No such density enhancements are initially found, as expected. Periodic enhancements show up during the phase dominated by phase mixing because of the compression driven by the ponderomotive force, but no persistent structures are formed. Overdense structures are formed after the KHi develops at later times. These structures appear as elongated filaments that follow the local direction of the twisted magnetic field. At the end of the simulation, fully formed long and thin threads that outline the magnetic field structure are observed. The high-resolution AMR scheme is crucial to correctly capture the formation of the finest threads.
The plasma compression that drives the thread formation increases the gas pressure as well. Pressures larger than 20% the background pressure are found inside the threads. As a consequence of this, the inner parts of the threads are slightly hotter than the background, but to a large extent the plasma can be considered to remain approximately isothermal (see again the bottom panel of Figure 1).

Effect of the Amount of Twist
Magnetic twist plays a key role because it allows the development of phase mixing in the initially uniform plasma. Here, we explore the effect of the twist parameter, c. Additional simulations have been performed with different twist parameters. In principle, the larger c, the steeper the radial profile of ω A . This should lead to the more rapid development of phase mixing and the quicker triggering of the KHi. Conversely, the larger c, the larger B j , which could have a stabilizing influence.
Following Díaz-Suárez & Soler (2021), to assess the KHi onset we use the vorticity squared integrated over the whole computational box, namely t t v r r , d .  However, a different behavior is found when c  0.75. For strong twist, Ω 2 saturates and decreases at larger times. The KHi does not seem to be excited for strong twist. Detailed analysis of the simulations reveals that the KHi is actually excited, but very locally, and the perturbations are not allowed to fully grow into the nonlinear regime. Such a behavior is shown in Figure 4. The larger B j can nonlinearly avoid the instability full development. As explained in, e.g., Galinsky & Sonnerup (1994), Ryu et al. (2000), Hillier (2019) magnetic tension is responsible for nonlinearly suppressing the KHi, preventing the growth of vortices, while the shear flows developed in the system can remain linearly unstable. Nevertheless, the failed KHi is able to disrupt the phase mixing process and to stop the increase of Ω 2 . Therefore, in our model the formation of density threads via full nonlinear KHi requires weak or moderate twist.

Forward Modeling
Here, we explore what kind of observational signature the overdense threads obtained in the simulations may leave. To this end, we used the FoMo code  to calculate the synthetic EUV emission for the SDO/    Figure 3), the KHi is locally excited, and vortices begin to form. The right panel corresponds to a latter time, when the vortices growth has been suppressed and Ω 2 has decreased.
AIA imaging telescope at 171 Å. The FoMo code requires input data, namely velocity components, temperature, and number density, in physical units. For this application, we take R = 5 Mm, B 0 = 10 G, and an initial temperature of 2 MK, as typical values in the corona (see, e.g., Priest 2014). This choice results in a reference Alfvén velocity of v A,0 ≈ 525 km s −1 and a number density of ≈2.88 × 10 9 cm −3 in the initially homogeneous background.
Following Van Doorsselaere et al. (2016), the specific intensity along the line of sight (LOS) is computed in the FoMo code as I x y x y z l , , , where κ α is the instrument response function for a bandpass α, Δl is the resolution along the LOS, x y z , , ¢ ¢ ¢ are the coordinates in a rotated frame of reference, and i, j, k are the indices in the rotated grid. The z¢-direction is the LOS direction whereas x¢ and y¢ are coordinates in the plane of sky (POS). With respect to the original simulation coordinates, the LOS is set along the y-direction, while z and x are the POS coordinates.
We obtained the specific intensity by feeding the FoMo code with the results of the reference simulation with c = 0.6. Since actual coronal magnetic fields anchored in the photosphere are curved, we mapped the specific intensity from the original straight grid onto a semicircular grid using the transformation Figure 5 shows the specific intensity in the original and transformed POS. As expected, the initially homogeneous plasma emits radiation uniformly. Early in the evolution, weak intensity variations appear due to the periodic plasma compression and expansion caused by the ponderomotive force. This effect drives longitudinal flows that are discernible in the intensity. Later, the KHi-induced threads become visible in the images. The threads appear as bright strands that illuminate and clearly reveal the underlying twisted magnetic structure. At the end of the simulation, the strands are about ∼5% brighter than the background plasma. The apparent thickness of the bright strands in the intensity images ranges from ∼150 km to ∼1500 km, approximately. The effective AMR resolution in this simulation is ∼10 km. The width of the finest strands is below the AIA 0 6 per pixel actual resolution (Lemen et al. 2012), but the thickest structures should be observable. The final simulation frame corresponds to a dimensional time of ∼20 minute after t = 0. Similar results are obtained for other AIA filters (not shown here), although the magnitude of the intensity variations and the apparent thickness of the strands may depend on the particular filter considered.

Discussion
In this Letter we report the results of 3D ideal MHD numerical simulations of standing TAW in a uniform plasma permeated by a line-tied magnetic field. Even with a uniform density, an Alfvén continuum is present because of the nonuniformity of the magnetic field that in our model is caused by twist, but magnetic expansion should have a similar effect (see Howson et al. 2019). The nonlinear evolution of phase mixing and the associated KHi break the uniformity of density, generating overdense structures in the form of elongated threads. These fine structures are locally aligned with the magnetic field direction. Forward modeling of SDO/ AIA imaging indicates that the formed threads could appear in observations as bright strands that illuminate the underlying coronal magnetic field.
The nonlinear evolution of TAW shares several features with that of kink MHD waves explored in the previous literature. The detailed studies of Antolin et al. (2014Antolin et al. ( , 2015Antolin et al. ( , 2016Antolin et al. ( , 2017 on kink oscillations of a preexisting denser tube embedded in a lighter environment show how KHi-induced rolls deform the cross-sectional area of the tube and mix the internal and external plasmas at the tube boundary. This dynamics produces strand-like structures in EUV intensity images as a result of complex plasma mixing of the rolls and their appearance depends on the LOS angle. The apparent strands last for timescales of a period of the kink oscillation (Antolin et al. 2014). The apparent strands produced during kink oscillations and the TAW-induced threads reported here share a common basic mechanism, i.e., the KHi, but this Letter introduces a new process. TAW do not require a preexisting tube, density enhancement, or waveguide. In a sense, TAW represent more fundamental oscillations of the magnetic structure. The threads obtained in the present simulations are actual density structures formed in an initially uniform plasma. We demonstrate that TAW evolving in a uniform corona are able to spontaneously generate field-aligned density enhancements that may further influence the structure and dynamics of the coronal medium.
We find that the threads are ∼10% denser than the background plasma, while the whole medium remains isothermal to a large extent. While 10% could be considered a small enhancement, it may importantly affect the subsequent evolution of the medium. The thermodynamic state in the corona is determined by a delicate balance between radiation losses, thermal conduction, and the unknown heating. Under the optically thin approximation, radiation losses in the coronal plasma are proportional to ∼ρ 2 , so that a 10% density increase would result in about ∼20% increase of energy losses. Equivalently, radiation losses should decrease in the regions where the density is depleted. If the background heating remains the same, the thread formation may evolve toward a situation of local thermal imbalance. We speculate that the overdense threads could behave as potential seeds of thermal instability, which is an essential ingredient of the coronal dynamics (see, e.g., Antolin 2020). This should be confirmed with nonideal simulations including thermal effects.
Finally, although we perform the simulations in a low-β scenario, the specific value of the plasma β might influence the strength of the KHi-driven compression that produces the overdense threads. A detailed analysis and parameter study of the role of the plasma β is needed. That is beyond the aim of this Letter and is left for forthcoming works.