Magnetic shocks and substructures excited by torsional Alfv\'en wave interactions in merging expanding flux tubes

Vortex motions are frequently observed on the solar photosphere. These motions may play a key role in the transport of energy and momentum from the lower atmosphere into the upper solar atmosphere, contributing to coronal heating. The lower solar atmosphere also consists of complex networks of flux tubes that expand and merge throughout the chromosphere and upper atmosphere. We perform numerical simulations to investigate the behaviour of vortex driven waves propagating in a pair of such flux tubes in a non-force-free equilibrium with a realistically modelled solar atmosphere. The two flux tubes are independently perturbed at their footpoints by counter-rotating vortex motions. When the flux tubes merge, the vortex motions interact both linearly and nonlinearly. The linear interactions generate many small-scale transient magnetic substructures due to the magnetic stress imposed by the vortex motions. Thus, an initially monolithic tube is separated into a complex multi-threaded tube due to the photospheric vortex motions. The wave interactions also drive a superposition that increases in amplitude until it exceeds the local Mach number and produces shocks that propagate upwards with speeds of approximately $ 50$ km s$^{-1}$. The shocks act as conduits transporting momentum and energy upwards, and heating the local plasma by more than an order of magnitude, with peak temperature approximately $60,000$ K. Therefore, we present a new mechanism for the generation of magnetic waveguides from the lower solar atmosphere to the solar corona. This wave guide appears as the result of interacting perturbations in neighbouring flux tubes. Thus, the interactions of photospheric vortex motions is a potentially significant mechanism for energy transfer from the lower to upper solar atmosphere.


INTRODUCTION
Magnetic flux tubes (and networks of flux tubes) are frequently observed in the solar atmosphere in an environment of relative pressure equilibrium with lifetimes lasting minutes, hours or even days.These stable magnetic configurations may act as waveguides transporting motions from the lower solar atmosphere into the upper chromosphere and corona.Waves propagating along such stable tubes have been well studied from observational, numerical and analytical approaches (for example, Bogdan et al. 2003;Banerjee et al. 2007;Terradas 2009;de Moortel 2009;Wang 2011;Mathioudakis et al. 2013;Jess et al. 2015;Nakariakov et al. 2016).
Observationally, a wide range of MHD wavemodes have been detected in the lower solar atmosphere, e.g.sausage (Morton et al. 2012), kink (He et al. 2009;Kuridze et al. 2013;Morton et al. 2014), and torsional Alfvén waves (Jess et al. 2009;Sekse et al. 2013).Understanding the propagation and mode conversion of these waves as they progress through the lower solar atmosphere gives insight into the magentic structure of the chromosphere and the locations of magnetic waveguides.
Vortex motions are highly important in revealing the fundamental dynamics of the solar atmosphere since they can significantly stress the magnetic field, drive the dynamics of the upper solar atmosphere and may contribute towards the heating of the solar corona.Photospheric vortex motions have been shown to be an effective mechanism for supplying mass and energy to the upper atmosphere (Wedemeyer-Böhm et al. 2012;Park et al. 2016;Murawski et al. 2018).Recent advances in the automated detection of these photospheric vortex motions (e.g., Kato & Wedemeyer 2017;Giagkiozis et al. 2017) indicate that these small-scale swirls are far more populous than previously thought and are capable of supplying energy to the upper solar atmosphere.Many of these vortex motions also exist in close proximity to each other, allowing for potential interactions of these motions, and this is the key motivation of the current work.
Photospheric vortex motions have been studied numerically in an expanding flux tube, and are found to excite a wide range of wave modes, including fast and slow magnetoacoustic, and Alfvén waves (Fedun et al. 2011b,a;Shelyag et al. 2013;Mumford et al. 2015;Mumford & Erdélyi 2015).Theoretical investigations of torsional Alfvén waves indicate that torsional motions can be reflected by the transition region and damped, resulting in heating (Giagkiozis et al. 2016;Soler et al. 2017).
Extending investigations towards more complicated flux tube structures (e.g., multiple flux tubes, merging flux tubes, or multistranded loops) increases the complexity of the wave interactions.Studies of localised perturbations inside a larger loop show that the resultant dynamics cannot be modelled as a monolithic loop, i.e. the interactions of multiple interior waves greatly affect the overall tube motion (Luna et al. 2010).Similar behaviour is expected when localised flux tubes, and their interior perturbations, interact and the resultant wave motions in the merged tube are a superposition of the isolated perturbations.2D numerical studies of merging flux tubes have shown that shocks can occur in the chromosphere (Hasan et al. 2005), and wave dissipation can contribute towards heating (e.g., Hasan & van Ballegooijen 2008;Vigeesh et al. 2012).Networks of multiple merging flux tubes can be constructed and stabilised following the work of Gent et al. (2013Gent et al. ( , 2014) ) whereby analytically stable networks are constructed using a realistic (VAL IIIC Vernazza et al. 1981) temperature and gas density profile.
In this paper, we investigate numerically the interactions of vortex motions in a pair of expanding and merging flux tubes with a realistic (VAL IIIC) temperature and plasma pressure profile.The individual tubes are excited at their bases using torsional velocity drivers that generate perturbations that propagate up the tube.When the flux tubes merge, these vortex motions interact, rearranging the magnetic field into a multithreaded structure, and driving high-velocity shocks that propagate into the upper solar atmosphere with speeds 50 km s −1 .This is a novel mechanism for driving spicules (or spicule-like structures) in the chromosphere via the interactions of vortex motions in merging magnetic flux tubes.

NUMERICAL CONFIGURATION
The numerical simulations are performed using the Sheffield Advanced Code (SAC, Shelyag et al. 2008).SAC solves the ideal MHD equations in a multidimensional system.The code was devised to simulate the linear and non-linear interactions of arbitrary perturbations in a gravitationally stratified and magnetised plasma atmosphere.A key feature of SAC is the separation of background and perturbation variables to allow for macroscopic processes to be modelled as a perturbation in a stable background atmosphere.The full system of equations solved for the perturbations to the background state are detailed in Appendix A.
A pair of identical, axisymmetric flux tubes are constructed using the self-similar approach outlined by Gent et al. (2013Gent et al. ( , 2014)).The initial temperature and density profiles are constructed using the VAL IIIC model for the lower solar atmosphere.This enables expanding flux tubes that capture the overall observed properties of some flux tubes, and are stabilised using analytical background forcing terms (see Appendix B).The physical justification for these forcing terms in that stable networks of flux tubes are regularly observed in the solar atmosphere, existing in relative pressure balance for extended periods of time.Flux tubes are con- structed that match observed properties and are stabilised using forcing terms that cannot be balanced by the scalar pressure or density gradients.The forcing terms account for unknown small-scale forces that cannot be measured, for example, the temporary cumulative effect of small-scale turbulence in the chromosphere and/or external forces acting from below the photosphere or in the neighbourhood of the flux tubes.The forcing terms manifest as terms in the momentum and energy equations.The flux tubes are constructed to match observed models (e.g.Verth et al. 2011;Jeffrey & Kontar 2013) and stable multiple flux tubes are regularly observed in pressure equilibrium (e.g.Levine & Withbroe 1977;McGuire et al. 1977;Malherbe et al. 1983).Perturbations in such stable flux tubes are investigated in a large number of observational, numerical and theoretical studies (see, for example, reviews cited in the introduction).The constructed atmosphere allows us to study the wave interactions in stable networks of more realistic flux tubes.The magnetic structure of each flux tube expands radially outwards with height following an exponential-type profile.The temperature and density profile along the tube axis is given by the VAL IIIC profile (Vernazza et al. 1981).The horizontal temperature and density profile is specified using the 3D magnetic field profile, creating a non-force-free equilibrium (see Appendix A.2.2 in Gent et al. 2014).The analytic construction of the flux tube pair are outlined in Appendix B.
At the photospheric level the centres of the flux tubes are located at (x, y) = (−1, 0), (1, 0) Mm and both have a base magnetic field strength of 1000 Gauss.In the lower atmosphere (z < 0.8 Mm) the flux tubes are independent in the sense that there is no overlap or interactions of the magnetic field or density profiles in the xy-plane.In the chromosphere, the two separate tubes begin to merge into one tube at z 0.8 Mm.The local magnetic pressure increases as the tubes merge, resulting in a decrease in plasma pressure, and a decrease in plasma−β, as shown in Figure 1.
The model extends vertically upwards to z = 2 Mm, where z=0 represents the top of the photosphere.This is below the transition region and allows us to focus solely on the effect of the merging of the tubes and the interactions of previously isolated motions in the merged tube.
Since the background and perturbation variables in SAC are separated, the background forcing terms outlined in Gent et al. (2014) appear, in their applicable form, only in the energy equation.The additional energy resulting from these forcing terms is small compared to the total energy.Test simulations were performed with the domain specified in the perturbation variables in SAC such that the derived atmosphere is advected.The atmosphere was tested to be stable for at least 1000 s and the current simulation was performed up to 400 s.
The computational grid spans the range −2 ≤ x ≤ 2, −1.4 ≤ y ≤ 1.4, 0 ≤ z ≤ 2 Mm and is resolved using a cell count of (100,100,200) for the x, y, and z directions, respectively.There is no evidence of significant numerical reflections from the boundaries in the simulation.Note that as time advances, numerical asymmetries form due to the representation of the physical domain on the numerical grid.The asymmetries remain small compared to the overall dynamics.
Vortex velocity drivers are specified near the base of the flux tubes.These are of the same form as their counterparts in Fedun et al. (2011b), i.e., which defines the velocity components v x , v y in terms of radius r = x 2 + y 2 and height z, for amplitude A 0 = ±1000 m s −1 , radial centre of the driver r 0 , time t and period P = 30 seconds.∆r = 0.2 Mm is used to define the radial expansion of the driver, and ∆z = 0.3 Mm defines the vertical driver size.The vertical centre of the driver is located at z 0 = 0.06 Mm which prevents the maximum driver amplitude from occurring on the z = 0 boundary.Drivers of equal magnitude and opposite direction of rotation are centred at (x, y) = (−1, 0) Mm and (1, 0) Mm, i.e., at the flux tube axes.This prevents a shear layer developing between the two flux tubes on the z = 0 boundary.

Propagation of waves in separate tubes
Below 1 Mm, the flux tubes are distinct and vortex motions are free to propagate independently, expanding with the tube and driving a number of different wave modes.The vortex drivers stress the magnetic field and transport energy and momentum upwards.The propagation of the isolated vortex motions are not discussed in detail since it has been well-studied previously (e.g.Fedun et al. 2011b;Mumford et al. 2015;Soler et al. 2017).Instead, we focus on the new physics introduced by the interactions of the flux tubes.

Merging of the flux tubes
At height z ≈ 0.8 Mm, the initially independent tubes merge together into one larger tube, see Figure 1.When this happens, the radial extent of the tube rapidly increases, enabling the previously independent vortex motions to expand and interact, as shown in Figure 2.
The vortices pass through each other, stressing the magnetic field and generating a rotation in the v × B term in Ohm's law.This in turn generates an outward force that reduces the mass density in the centre of the domain.
During this phase, the vortex interactions also create a series of short-lived thin magnetic substructures, see Figure 3.These features exist for a relatively short time and are evidence of reorganisation of the magnetic field due to torsional motions.These features form and disappear sporadically as time progresses, generating a number of magnetic substructures, i.e., thin magnetic structures contained within the initial tube, see Figure 3.These magnetic substructures are co-located with peaks in vertical velocity on the order of 1 km s −1 and are effectively waveguides transporting energy upwards.The substructures are typically ≤ 0.4 Mm in width and drift horizontally away from the centre of the merged tube.Magnetic substructures also form further up the tube as time progresses (Figure 3).The structures are clearly evident in the local variations in magnetic field strength, as illustrated, but hardly distinguishable in the distribution of the plasma density.Observationally, this would mean that such fine structures are difficult to identify in intensity maps, despite having enhanced localised velocity and magnetic field perturbations.The formerly monolithic tube has become magnetically multi-stranded as a result of the applied photospheric vortex motions.
The formation of these magnetic structures is dominated by advection.SAC possesses a small (but nonzero) numerical diffusion that can allow the magnetic field to be reorganised.To assess the importance of diffusion in our simulation we compare the v × B and ηJ  terms from Ohm's law, using a conservative estimate of η.It is found that the advection term is several orders of magnitude above the diffusion term.Therefore, whilst some, non-zero diffusion is present in the simulation, the reorganisation of the magenetic field is dominated by advection, i.e. ideal (η ≈ 0) processes.
The perturbations begin to affect the magnetic field near the upper z boundary at t ≈ 254 s, as shown in Figure 3.There is no evidence of significant numerical reflections from the upper boundary.

Shock formation
The vortex interactions generate a superposition where the flux tubes merge.The v z velocity amplitude at this point increases until it exceeds the sound and Alfvén speeds (at the merge point, plasma-β ≈ 1), driving shocks.A time-series of this increasing amplitude wave developing into a shock is shown in Figure 4.Note that this figure is not at the lowest formation region of the shock.The increase in amplitude is a result of the continued stress created in this region from the torsional drivers.The 3D structure of the shock is approximately conical and there are no rotational or helical motions in the shock itself.The background atmospheric conditions change as the waves propagate upwards; the plasma-β and plasma density drop and the shock separates into magnetic and hydrodynamic components, resulting in two shock fronts propagating at different speeds.This is shown by the separation of sonic and Alfvénic Mach numbers at t 300 in Figure 4.

Energy transfer
Vortex motions have been shown to transport energy and mass upwards in the solar atmosphere (e.g., Soler et al. 2017) and it has also been shown that rotational structures such as tornadoes can contribute significantly to the heating of the solar corona (Wedemeyer-Böhm et al. 2012).From our numerical simulation, we quantify the energy transported to the upper chromosphere as a result of vortex motions, and the energy change from the plasma evacuation and shock.We define the different types of energy as follows: Magnetic Energy for mass density ρ, velocity v, magnetic field strength B, and thermal pressure P .The universal constants are the permeability of free space µ 0 = 4π × 10 −7 , and the specific gas ratio γ = 5/3.In particular, the integrated energy at each height in the simulation reveals the energy transport upwards.The normalised integrated energy at time t = 360 seconds is shown in Figure 6 for the three different energies, Equations ( 3)-( 5).Note that the internal energy is far larger than the other energies and as such the total energy ≈ .Energy is shown as a percentage increase from time t = 0 except for kinetic energy where k e (t = 0) = 0 by definition.Kinetic energy is normalised by its value at z = 0.1 Mm.
The bulk energy remains near the photosphere and horizontal attenuation limits the amount of kinetic energy that can propagate upwards.The vortex drivers supply velocity and stress the magnetic field hence the kinetic and magnetic energies have peak values close to the lower boundary.The velocity increases with z (see Figure 5), however, the density decreases with z resulting in little apparent increase in kinetic energy along the domain length despite the large increase in velocity.The integrated total energy (k e + m e + ) in the upper chromosphere occurs at z = 1.2 Mm and the total energy here increases by approximately 20%.Note that the shock creates a localised increase in energy which is averaged when we visualise the total energy along a slice.
At the core of the shock, there is a large increase in temperature (see Figure 7).At height z = 1.3 Mm and time t = 360 s, the temperature increases to 60, 000 K, an increase of an order of magnitude.The heating is localised to the shock and the temperature near the edge of the merged tube remains of a similar order of magnitude to the initial condition.There is a corresponding decrease in density at the heating location and the overall transverse density structure through the shock is shown in Figure 7.The interior of the shock is a reduced density region (with locally enhanced temperature) and an increased density at the shock edge.The high density on the edge of the shock reduces the temperature of the plasma through the ideal gas law.Note that the asymmetries present in Figure 7 are numerical in nature and originate from the representation of the physical domain on the numerical grid.

DISCUSSION
In this paper, we investigate the interactions of photospheric vortex motions in a pair of expanding and merg- ing flux tubes.Stable flux tubes are constructed that are independently perturbed by counter-rotating vortex motions at their footpoints.These perturbations interact both linearly and nonlinearly driving several features that reorganise magnetic field and transport energy throughout the system.
Formation of magnetic substructures -We have shown that vortex motions in a simple pair of magnetic flux tubes are capable of reorganising the magnetic field and forming a myriad of smaller magnetic substructures (Figure 3).Thus the initially monolithic tube becomes multi-threaded.The torsional motions stress the magnetic field causing this development of localised flux tubes inside the larger structure.These substructures can act as waveguides to transport energy and momentum to the upper solar atmosphere.In the solar chromosphere, spicules are observed to split and merge (e.g., Sterling et al. 2010a), however, the possible mechanism(s) that cause this are still not well understood.In our simulations, we demonstrate that torsional motions from adjacent and merging flux tubes interact creating smaller substructures.
Shock formation -The interacting vortex motions create a superposition in the centre of the domain, where the tubes merge.The continued driving increases the amplitude of the superposition until it exceeds the sound and Alfvén speeds and shocks.This shock propagates upwards with a speed of 50 km s −1 transporting energy into the solar corona.The development of this shock indicates a potential way of driving spicules and chromospheric jets via photospheric vortex motions.
Energy transfer -The presented model is a potentially efficient mechanism for transporting energy from the lower to upper solar atmosphere.The spatiallyintegrated energy over the (x, y)-plane in the upper chromosphere was found to increase by up to 20%, with the localised energy being even greater.The shocks heat the plasma to ≈ 60, 000 K in the upper chromosphere.This heated plasma propagates upwards and therefore can supply energy and momentum to the upper atmosphere.

Figure 1 .
Figure 1.(Left) 3D isosurfaces of plasma−β and arbitrary magnetic fieldlines (grey).At the photosphere, the two flux tubes are independent (i.e., there is no overlap or interaction in the density or magnetic field profiles) and have footpoints located at (x, y) = (1, 0), (−1.0)Mm.The flux tubes expand with height and merge into a single tube at z 0.8 Mm. (Right) Slice through the centre of the flux tubes showing the plasma−β (colormap) and magnetic field lines (blue lines).Selected contour lines of the plasma-β are overplotted to highlight the key plasma-β regimes throughout the system.

Figure 2 .
Figure 2. Colour plot of velocity component (vz) at z = 1 Mm at times t = 180 (left), 200 (centre), and 220 (right) seconds.Vector of the in-plane Lorentz force (v × B) is shown as arrows.

Figure 3 .
Figure 3. (Left column) Snapshots showing line-contours of the magnetic field strength at heights z = 0.18, 1.0 and 1.91 Mm at times 108, 152, 254, and 302 s, respectively.z = 0 colourmap shows the photospheric magnetic field.Left and rear colourmap show the mass density.(Right column) Corresponding slices at z = 1.0 Mm showing magnetic field line contours, colourmap of vz velocity, and vectors of in-plane velocity.

Figure 4 .
Figure 4. Alfvénic (black solid) and sonic (red dashed) Mach numbers at the centre of the domain (x = y = 0, z = 1.0 Mm) through time.This point is located above the merge point of the flux tubes.

Figure 5 .Figure 6 .
Figure 5. Vertical velocity (vz [km/s]) colour plot showing the high-velocity upflow region created in the centre of the domain at time t = 360 s.Arbitrary streamlines of magnetic field show the overall magnetic structure.

Figure 7 .
Figure 7. Mass density of the tube at height z = 1.3 Mm.Initial profile (t = 0) given by red dashed line.Black solid line is the mass density across the shock at time t = 360 s.Temperature at t = 360 s is shown by the green dash-dot line.