© ESO 2019

1. Introduction

Ever since the discovery of the high temperatures of the solar corona by Edlén (1943), scientists have been trying to find an explanation for this phenomenon. One possible coronal heating mechanism is alternating current (AC) heating, in which magnetic energy is dissipated by waves (Aschwanden 2006; Priest 2014; Parnell & De Moortel 2012; Arregui 2017). AC heating mechanisms were widely ignored for a long time until it was observed that waves are indeed ubiquitous in the solar corona (De Pontieu et al. 2007; Tomczyk et al. 2007; Krishna Prasad et al. 2012; Morton et al. 2012; Nisticò et al. 2013).

Tomczyk & McIntosh (2009) observed waves with the Coronal Multi-channel Polarimeter (CoMP) instrument and found that the spectrum of velocity perturbations peaks at the same frequency (∼3 mHz) as solar p-modes. Morton et al. (2016, 2019) confirmed that the power enhancement at this frequency is a global phenomenon. Therefore, it seems likely that the ubiquitous waves are at least partially driven by p-modes. However, it is not yet clearly understood how p-modes propagate into the higher atmosphere. Observations tracing p-modes in the chromosphere and above were done by, for example, Centeno et al. (2006), Marsh & Walsh (2006), de Wijn et al. (2009), Prasad et al. (2015), and Zhao et al. (2016), while numerical modeling was done by, for example, Khomenko et al. (2008), Fedun et al. (2011), Santamaria et al. (2015), and Griffiths et al. (2018). De Pontieu et al. (2004, 2005) found that, although p-modes with periods above the cutoff period are usually evanescent in the chromosphere, these p-modes can propagate upward along inclined magnetic flux tubes because the effective cutoff period is increased in a non-vertical magnetic field.

It was shown by Bogdan et al. (1996), Hindman & Jain (2008), and Gascoyne et al. (2014) that p-modes lose energy to magnetic tube waves, such as sausage waves and kink waves; the propagation of these excited waves into the solar atmosphere was not considered in those studies. Observations show that tube waves are indeed excited in the solar atmosphere. While propagating kink waves have already been observed many years ago Verwichte et al. (2005), and Grant et al. (2015) were more recently able to observe sausage modes propagating from the photosphere to the transition region for the first time. These authors find surprisingly strong damping for those sausage waves; this is not well understood and therefore indicates again our limited knowledge of wave propagation in that region.

We therefore stress the importance of understanding the propagation of p-modes through the chromosphere, as they might be the source of decayless transverse waves. Those waves could be connected to coronal heating (e.g., Karampelas et al. 2017). They have been observed by Wang et al. (2012), Nisticò et al. (2013), and Anfinogentov et al. (2013), for example. Transverse waves are currently modeled in the corona as loops with a horizontal driver in one or both footpoints (Karampelas & Van Doorsselaere 2018; Guo et al. 2019; Karampelas et al. 2019; Pagano & De Moortel 2019), but from where those horizontal plasma movements originate is usually not discussed. One possible mechanism could be a self-oscillatory process due to the interaction of the loops with quasi-steady flows, as discussed by Nakariakov et al. (2016). Another possibility, as mentioned above, could be photospheric p-modes, which are converted to kink waves.

In this work, we study the conversion of photospheric p-modes to sausage and kink waves in an atmosphere that is gravitationally stratified and has additionally structuring perpendicular to the magnetic field. The model we use is in magnetohydrostatic (MHS) equilibrium and contains four loops with different inclinations ranging from the photosphere to the lower corona. We perturb the equilibrium at the bottom with a p-mode driver in the form of an analytic solution for gravity-acoustic waves and simulate the propagation of the waves. The waves interact with the cylindrical structure of our model atmosphere and tube modes are excited. The MHS model is described in Sect. 2, while we explain the numerical setup in Sect. 3. We present assisting methods for data interpretation in Sect. 4 and present and discuss our results, including what kind of wave modes are excited in the corona and how their basic properties change, in Sect. 5. Finally, we present the summary and conclusions of our work in Sect. 6.

2. Magnetohydrostatic equilibrium model

We built a 3D MHS equilibrium atmosphere from the photosphere to the lower corona, which has to fulfill the condition

$$\begin{array}{c}\hfill \mathbf{\nabla }{p}_0-\frac{1}{{\mu }_0}(\mathbf{\nabla }\times {\mathit{B}}_{\mathbf{0}})\times {\mathit{B}}_{\mathbf{0}}-{\rho }_0\mathit{g}=0,\end{array}$$(1)

where p₀, B₀, and ρ₀ are the equilibrium pressure, equilibrium magnetic field, and equilibrium density, respectively, and μ₀ is the permeability of vacuum. For our models, the gravity vector g=(0,0,−g) is constant and points to the negative z-direction. In addition, we built the atmosphere as periodic in the horizontal directions. The domain has a size of n_x × n_y × n_z = 140 × 140 × 840 points with a resolution of Δx × Δy × Δz ≈ 14.3 × 14.3 × 6.0 km, which results in a domain with the approximate measurements of 2 × 2 × 5 Mm. The bottom seven planes of grid cells of the domain are located below the photosphere (z = 0) to make space for the driver (see Sect. 3).

In the first step we define a divergence-free magnetic field, where the total field strength forms several straight loops of reduced magnetic field with a gauss-shaped cross section

$$\begin{array}{c}\hfill B_{0,x}=B_{0,y}=0,\end{array}$$(2)

$$\begin{array}{c}\hfill B_{0,z}=a[1-\sum _i^{n}exp\{-\frac{{(x-x_{0,i})}^2+{(y-y_{0,i})}^2}{{\sigma }_i^2}\}]+b.\end{array}$$(3)

The sum corresponds to a sum over all loops, where x_0, i and y_0, i describe the coordinates of the loop centers and σ_i define the thickness of each loop. We place four loops evenly inside the domain with a distance of 1 Mm in the x- and the y-direction between the loops and use σ = 1/3 Mm for all i. The constants a = 309 G and b = 5 G define the strength of the magnetic field. The magnetic field has its minimum in the loop centers with value b, whereas the theoretical maximum outside the loops is a + b. However, owing to the tight structuring of our loops, this maximum is never reached. The resulting magnetic field ranges from 5 G to 300 G and is shown in Fig. 1 at the top left. Since we only have a magnetic field component in the z-direction and the magnetic field does not change with height, ∇ ⋅ B₀ = 0 is automatically fulfilled.

Fig. 1. Magnetohydrostatic model of four loops for θ = 0° (top row) and θ = 15° (bottom row). Left column: magnitude of the total magnetic field. Middle column: pressure (logarithmic scaling). Right column: temperature. All length scales are in Mm.

To study the effect of an inclined magnetic field as well, we rotate the vertical loop system from Eqs. (2) and (3) clockwise around the y-axis by an angle θ, which modifies the equations to

$$\begin{array}{cc}\hfill B_{0,x}& =\stackrel{\sim }{B}sin(\theta ),\hfill \end{array}$$(4)

$$\begin{array}{cc}\hfill B_{0,y}& =0,\hfill \end{array}$$(5)

$$\begin{array}{cc}\hfill B_{0,z}& =\stackrel{\sim }{B}cos(\theta ),\hfill \end{array}$$(6)

with

$$\begin{array}{c}\hfill \stackrel{\sim }{B}=a[1-\sum _i^{n}exp\{-{\varphi }_i\}]+b,\end{array}$$(7)

and

$$\begin{array}{c}\hfill {\varphi }_i=\frac{{(-sin\left(\theta \right)z+cos\left(\theta \right)x-x_{0,i})}^2}{{\sigma }_i^2}+\frac{{(y-y_{0,i})}^2}{{\sigma }_i^2}·\end{array}$$(8)

We note that to keep the domain periodic, we extend it in the x-direction and therefore also slightly change the corresponding resolution. That leads to a number of points in the x-direction of n_x = 145 for θ = 15°. The total magnetic field for the case of θ = 15° is shown in Fig. 1 at the bottom left.

In order to calculate the pressure and density we split Eq. (1) into its three components and slightly reorder these components as follows:

$$\begin{array}{c}\hfill \mathbf{I}:B_{0,y}{\partial }_y{B}_{0,x}+B_{0,z}{\partial }_z{B}_{0,x}-B_{0,y}{\partial }_x{B}_{0,y}-B_{0,z}{\partial }_x{B}_{0,z}={\mu }_0{\partial }_x{p}_0,\end{array}$$(9)

$$\begin{array}{c}\hfill \mathbf{II}:B_{0,x}{\partial }_x{B}_{0,y}+B_{0,z}{\partial }_z{B}_{0,y}-B_{0,x}{\partial }_y{B}_{0,x}-B_{0,z}{\partial }_y{B}_{0,z}={\mu }_0{\partial }_y{p}_0,\end{array}$$(10)

$$\begin{array}{c}\hfill \mathbf{III}:B_{0,x}{\partial }_x{B}_{0,z}+B_{0,y}{\partial }_y{B}_{0,z}-B_{0,x}{\partial }_z{B}_{0,x}-B_{0,y}{\partial }_z{B}_{0,y}={\mu }_0{\partial }_z{p}_0+{\rho }_{0}g{\mu }_0,\end{array}$$(11)

where ∂x_j = ∂/∂x_j, j = 1, 2, 3. Reforming and integrating the first two components gives us

$$\begin{array}{c}\hfill p_0={\stackrel{\sim }{p}}_{\mathrm{I}}+f_1(y,z)={\stackrel{\sim }{p}}_{\mathrm{II}}+f_2(x,z)={\stackrel{\sim }{p}}_{\mathrm{I}}+h(z),\end{array}$$(12)

where ${\stackrel{\sim }{p}}_{\mathrm{I}}$ and ${\stackrel{\sim }{p}}_{\mathrm{II}}$ are the pressures calculated from integrating Eqs. (9) and (10), respectively, and f₁ and f₂ are functions resulting from the integral that have to be determined. For our magnetic field model ${\stackrel{\sim }{p}}_{\mathrm{I}}={\stackrel{\sim }{p}}_{\mathrm{II}}$, which leads to the last equality of Eq. (12). The resulting expression for the pressure is

$$\begin{array}{c}\hfill p_0=-\frac{1}{2{\mu }_0}{\stackrel{\sim }{B}}^2+h(z)·\end{array}$$(13)

The function h(z) is arbitrary, as it has no influence on Eqs. (9) and (10), and represents the vertical pressure stratification, and the constant that has to be added to make Eq. (13) positive. However, we stress the significance of choosing h(z) wisely, as this term essentially determines our vertical density profile. Therefore, we define the vertical pressure stratification according to the VAL-C model and add an exponential term to modify the stratification

$$\begin{array}{cc}\hfill h(z)& =p_{\mathrm{bot}}exp\{-{\displaystyle \int \frac{1}{H(z)}\mathrm{d}z}\}\hfill \\ \hfill & +150[\mathrm{Pa}]exp\{-\frac{z}{600·6046[\mathrm{m}]}+0.015\}+422[\mathrm{Pa}]·\hfill \end{array}$$(14)

In this case, p_bot = 14 kPa is the pressure at the bottom of the domain for the first term of Eq. (14). The scale height H(z) is calculated with a temperature profile that follows the VAL-C model until the transition region and has a constant temperature for the corona. The two regions of the temperature profile are connected by a cosine-shaped transition region. Since this initial temperature profile is only used for calculating the vertical pressure stratification and is changed in the next step, we abstain from mentioning the exact expression. The constants in the exponential term of Eq. (14) are related to our practical implementation and their exact values have no deeper meaning. Figure 2 (left) shows the pressure profile as a function of height for θ = 0° for pressure according to the VAL-C model (dashed lines) and for pressure with added exponential term, as used in our model (solid lines). Both profiles have the same start and end points, but in the modified version the slope in the upper part of our model is larger, which leads to higher density values in that region. This is necessary because as a consequence of the constant magnetic field along the loops the density would be very low otherwise, which would lead to unreasonably high temperatures. However, the additional term also leads to a strong broadening of our transition region, which we deem a compromise that has to be made.

Fig. 2. Profiles and profile ranges of the atmospheric quanities for the vertical case as a function of z. The lines indicate the minimum and maximum values at each height, respectively. The striped region indicates the transition region. Left: pressure according to the VAL-C model (orange dashed lines) and with added exponential term (red solid lines). Right: temperature (green solid lines) and density (blue dashed line, constant in x and y).

As soon as the magnetic field and the pressure are known it is straight forward to get an expression for ρ from Eq. (11). Finally, we calculate the temperature from the ideal gas law. Figure 2 (right) shows the density (dashed) and temperature (solid) profiles as a function of z for the vertical case, whereas Fig. 1 shows a 3D plot of the pressure and temperature for both the vertical and the 15° inclined case. We note that the pressure is vertically much more stratified than horizontally structured.

Since the pressure is higher in the loop interior, where the magnetic field is lower, we get a sound speed profile that is much higher inside the loops than outside. On the contrary, the Alfvén speed has its maximum outside the loops, while being very low in their center. This would be expected for a coronal loop with higher density than its surroundings, however, to fulfill Eq. (1) the density is horizontally constant in our model. This is easily visible in Eq. (11) considering the vertical case with B_0, x = B_0, y = 0 and taking into account that ∂_zp₀ is independent of x and y for θ = 0°. The pressure and magnetic field distributions lead to a horizontally strongly structured plasma-β profile. The β = 1 contour is plotted in Fig. 3. As we go up in the atmosphere, the plasma-β decreases in the loop exterior already below the transition region to lower than unity, while it is always much higher than unity in the loop interior until the top of the domain.

Fig. 3. Contour of β = 1 for θ = 0° (left) and θ = 15° (right). At the bottom of the domain (photosphere) β ≫ 1.

3. Numerical setup

For our simulations we used the MANCHA3D code (Khomenko & Collados 2006; Felipe et al. 2010; Khomenko et al. 2018) developed at the Instituto de Astrofísica de Canarias in Tenerife, Spain. This code solves the fully nonlinear magnetohydrodynamics equations for perturbed variables, which is why we initially had to define a MHS equilibrium to perturb it in the simulations. For the present work we consider an adiabatic system and neglect changes of the mean molecular weight due to ionization.

On the vertical faces of the computational box we set periodic boundary conditions. Our system can therefore be viewed as an ensemble of thin loops (or loop strands) extending infinitely to all horizontal sides that reasonably represent groups of spicules. To allow waves to escape, we set a Neumann-type zero-gradient open boundary condition at the upper boundary. At the first seven layers of the lower boundary we applied a driver that follows an analytic solution for a vertical gravity-acoustic wave (Mihalas & Mihalas 1984). The driver is described in detail by Santamaria et al. (2015) and in this work we only repeat the general form,

$$\begin{array}{c}\hfill {\mathit{v}}_{z,1}=V_{0}exp\{\frac{z}{2H}+k_{\mathit{zi}}z\}sin(\omega t-k_{\mathit{zr}}z),\end{array}$$(15)

$$\begin{array}{c}\hfill \frac{p_1}{p_0}=V_0|P|exp\{\frac{z}{2H}+k_{\mathit{zi}}z\}sin(\omega t-k_{\mathit{zr}}z+{\varphi }_{\mathrm{P}}),\end{array}$$(16)

$$\begin{array}{c}\hfill \frac{{\rho }_1}{{\rho }_0}=V_0|R|exp\{\frac{z}{2H}+k_{\mathit{zi}}z\}sin(\omega t-k_{\mathit{zr}}z+{\varphi }_{\mathrm{R}}),\end{array}$$(17)

where v_z, 1, p₁, and ρ₁ are the velocity perturbation in the z-direction, pressure perturbation, and density perturbation, respectively. The values |P| and |R| are the relative amplitudes, while ϕ_P and ϕ_R are the phase-shifts of pressure and density perturbation compared to the velocity perturbation. The value H is the pressure scale height and V₀ is the amplitude of the velocity perturbation. The vertical wave number k_z is either complex or real, depending on the frequency ω compared to the isothermal acoustic cutoff frequency ω_c,

$$\begin{array}{c}\hfill k_z=k_{\mathit{zr}}+i{k}_{\mathit{zi}}=\frac{\sqrt{{\omega }^2-{\omega }_{\mathrm{c}}^2}}{c_{\mathrm{s}}}\end{array}$$(18)

with

$$\begin{array}{c}\hfill {\omega }_{\mathrm{c}}=\frac{\gamma g}{2c_{\mathrm{s}}},\end{array}$$(19)

where c_s is the sound speed and γ = 5/3 is the adiabatic index.

In this work we used a small perturbation of V₀ = 10⁻² m s⁻¹ to stay in the linear regime. In addition, we used a period of 100 s, which leads to a frequency of ω ≈ 0.063 rad s⁻¹. The main reason for choosing this small period is that the imperfect open boundary conditions at the top cause waves to reflect and propagate downward. The chosen period allows us to study at least half of a wave period before the reflected waves interfere. For this period, ω >  ω_c is valid for the whole domain, so the waves excited from the bottom of the domain never reach a cutoff region. To put it into a solar context, the p-modes in our simulations are no longer trapped within the solar interior and can propagate through the chromosphere. Since the waves are not trapped within a resonant cavity, their amplitude is not amplified by constructive interference, which validates our choice of a small V₀. We plan to use larger periods that allow a cutoff region in future work.

4. Methods for data interpretation

4.1. Decomposition into components

In order to distinguish the different wave modes, it is necessary to bring our simulation data into a form that allows us to visualize the characteristics of the expected modes. Tarr et al. (2017) decomposed their data into kinetic, acoustic, and magnetic energy densities, which allowed these authors to decouple fast from slow waves and magnetic from acoustic waves. Another decomposition method was carried out by Khomenko et al. (2018), who, following Cally (2017), constructed three quantities based on the physical properties of the waves: (1) f_alf for the incompressible perturbation propagating along the magnetic field, (2) f_long for the compressible perturbation propagating along the magnetic, and (3) f_fast for the compressible perturbation perpendicular to the magnetic field. While f_alf is a quantity describing the Alfvén waves for all β, f_long and f_fast decouple the slow and the fast magneto-acoustic waves only for β <  1.

However, since we have both β <  1 and β >  1 regions in our model and expect tube modes to be excited, we adopted the approach of Mumford et al. (2015), who split the velocities and fluxes into three orthogonal components defined by the magnetic flux surfaces. These components are defined by the unit vectors longitudinal (${\widehat{\mathit{e}}}_{\Vert }$), azimuthal (${\widehat{\mathit{e}}}_{\mathrm{a}}$), and normal (${\widehat{\mathit{e}}}_{\perp }$) to those surfaces. Since we only used small perturbations, magnetic flux surfaces and therefore also the resulting unit vectors for the three components are constant in time. The unit vector for the longitudinal component ${\widehat{\mathit{e}}}_{\Vert }$ points into the direction of the magnetic field and is therefore easy to compute. Because all field lines in our model are straight, this vector is the same for all points of the domain.

Calculating the azimuthal unit vector ${\widehat{\mathit{e}}}_{\mathrm{a}}$ proves to be more difficult. We solve it by calculating the 2D isocontour of the magnetic field for all horizontal layers and fitting a straight line to the contour for each pixel. The direction of the linear fit is then the direction of ${\widehat{\mathit{e}}}_{\mathrm{a}}$. For the inclined case ${\widehat{\mathit{e}}}_{\mathrm{a}}$ gets an appropriate vertical component to be normal to ${\widehat{\mathit{e}}}_{\Vert }$. However, with this method all azimuthal vectors point to the positive x-direction, which results in half of the vectors being clockwise (top half of the loops in respect to y), while the others are counterclockwise (bottom half of the loops in respect to y). Because of the regular structure of the loop system it is simple to distinguish those regions. We multiply all clockwise azimuthal unit vectors with −1 to have a consistent sense of direction. Finally, the normal unit vector ${\widehat{\mathit{e}}}_{\perp }$ is calculated by ${\widehat{\mathit{e}}}_{\perp }={\widehat{\mathit{e}}}_{\mathrm{a}}\times {\widehat{\mathit{e}}}_{\Vert }$. With this convention, ${\widehat{\mathit{e}}}_{\Vert }$ always points upward (and for the inclined case also into the positive x-direction), ${\widehat{\mathit{e}}}_{\mathrm{a}}$ is parallel to the flux surfaces and points to the counterclockwise direction, and ${\widehat{\mathit{e}}}_{\perp }$ points away from the loop centers. Figure 4 sketches the vector directions for vertical and inclined loops.

Fig. 4. Orthogonal decomposition vectors parallel (${\widehat{\mathit{e}}}_{\Vert }$), normal (${\widehat{\mathit{e}}}_{\perp }$), and azimuthal (${\widehat{\mathit{e}}}_{\mathrm{a}}$) to the magnetic iso-surfaces (cylindrical shapes) for a vertical (left) and inclined (right) loop.

4.2. Expected wave modes

Since our atmosphere in not horizontally uniform but has a cylindrical shape, we would also expect the wave modes excited in our simulations to be those of a plasma with cylindrical shape, such as the m = 0 sausage mode and m = 1 kink mode, where m is the azimuthal wave number. To approximately calculate what wave modes would appear in our setup, we assume a simple vertically constant cylinder with radius R and internal values f_i embedded in an external plasma with values f_e. Similar to our atmosphere, we assume that the magnetic field is parallel to the loop axis. We then followed the mathematical framework of Moreels & Van Doorsselaere (2013). As external and internal values for this simplified model we used the average external and average internal values of our atmosphere at a height of 2 Mm, where the loop boundary is defined by the β = 1 layer. If c_s, i is the internal sound speed, this leads to an external sound speed of c_s, e = 0.745c_s, i, an internal Alfvén speed of c_A, i = 0.618c_s, i, and an external Alfvén speed of c_A, e = 0.994c_s, i. Figure 5 shows the resulting phase speed diagram, where the internal and external sound and Alfvén speeds are indicated by horizontal gray lines. Also plotted are the internal and external tube speed

$$\begin{array}{c}\hfill c_{\mathrm{T},\mathrm{f}}=\frac{c_{\mathrm{s},\mathrm{f}}{c}_{\mathrm{A},\mathrm{f}}}{\sqrt{c_{\mathrm{s},\mathrm{f}}^2+c_{\mathrm{A},\mathrm{f}}^2}}\end{array}$$(20)

Fig. 5. Phase speed diagram of sausage modes and kink modes for conditions similar to the model atmosphere at a height of 2 Mm. The horizontal gray lines indicate various characteristic speeds.

and the kink speed

$$\begin{array}{c}\hfill c_{\mathrm{k}}=\sqrt{\frac{{\rho }_{0,\mathrm{i}}{c}_{\mathrm{A},\mathrm{i}}^2+{\rho }_{0,\mathrm{e}}{c}_{\mathrm{A},\mathrm{e}}^2}{{\rho }_{0,\mathrm{i}}+{\rho }_{0,\mathrm{e}}}}·\end{array}$$(21)

For both sausage and kink waves, we only get non-leaky solutions for fast surface waves and slow body waves. Slow surface waves are theoretically possible below c_T, i, but no solutions are found in this region. Fast body modes could occur above c_s, i, but they are leaky; we indeed find leaky solutions for a kink body mode there for higher kR. The phase speed line of the fast sausage surface mode (red dotted line) actually stops where kR is ≈0.8 as there are no solutions found for low kR, not even for leaky waves.

It is now possible to calculate the magnitude of the ratio of longitudinal displacements to perpendicular displacements for the modes shown in Fig. 5. Unlike the normal component in the decomposition we described in Sect. 4.1, perpendicular displacement describes plasma displacement in all directions perpendicular to the loop axis, so ${\xi }_{\mathrm{perp}}=\surd ({\xi }_{\mathrm{a}}^2+{\xi }_{\perp }^2)$. This ratio depends on the distance to the loop axis r and is shown in Fig. 6 at r = 0.69R (left) and at r = R (right). From the figures it is clear that, although the curves slightly change, the general behavior of the displacement ratios stays the same regardless of the chosen point inside the loop or the loop surface. While the slow body modes and the fast sausage surface mode have much higher parallel displacements than perpendicular displacements for small kR compared to larger kR, it is the opposite for the fast kink surface mode. This makes the fast sausage surface mode and the fast kink surface mode easily distinguishable from each other.

Fig. 6. Ratio of longitudinal displacement and perpendicular displacement of plasma for sausage and kink modes for conditions similar to the model atmosphere at a height of 2 Mm. Left: inside the loop at r = 0.69 R. Right: at the surface of the loop (r = R).

5. Results and discussion

5.1. Wave propagation for vertical flux tubes

We first study simulations with the vertical (θ = 0°) case. Our goal is to investigate the conversion of p-modes that arrive at the corona. For that purpose we look at the evolution of a horizontal cut through the domain at a height of 2 Mm. We later show that we have some issues with wave reflection from the upper boundary, so investigating the wave behavior at 2 Mm instead of higher up allows us to analyze a longer time sequence before the reflections from the upper boundary intervene. This is possible because the wave behavior does not change much above the transition region after 2 Mm, which is the case even though the high β regions inside the loops are still merged together at that height. The movie showing the time development of the horizontal cut before the reflections from the upper boundary arrive (t ≤ 170 s) is available online. A screenshot of the movie is shown in Fig. 7. It shows the three velocity components at t = 106 s together with the β = 1 contour. For this time series the longitudinal velocity perturbation is always larger inside the loop than outside. The waves arrive at 2 Mm at approximately the same time; the waves inside the loop arrive slightly earlier. The maximum perturbations of the normal component are more than one order of magnitude smaller than the maximum perturbations of the parallel component. In the first part of the time series, the normal component has a positive sign everywhere and thus shows an expansion of the whole loop cross section with some normal velocity components outside the loops as well. In its first maximum at t = 106 s (as indicated in Fig. 7) the expansion of the loop has similarities with a m = 4 fluting mode. However, this mode would also require some plasma to flow into the loops at the top, bottom, and sides of the loops. At around t = 114 s the normal component changes its sign and at t = 128 s it looks the same as in Fig. 7 but with changed sign (i.e., contraction instead of expansion). We therefore conclude that the wave the p-modes excited is in fact a m = 0 sausage mode that is deformed by the tight packing of the grid-like positioned loops. The deformed sausage modes are very similar to a superposition of a m = 0 sausage mode with a m = 4 fluting mode.

Fig. 7. Components of the velocity perturbation in a horizontal cut at 2 Mm for θ = 0° at a time of 106 s after the start of the simulation. The velocities are given in m s−1 and the spatial scales are in Mm. The black lines show the β = 1 contour, with β ≫ 1 inside and β <  1 outside the loops. The values of the first two pixels next to the margin are set to zero for the normal and azimuthal component, as the corresponding vectors were badly defined in that region. The temporal evolution is available as an online movie.

Immediately apparent in the azimuthal component of Fig. 7 is the ring-like structure of counter-flowing plasma velocities close to the β = 1 layer. However, this ring structure propagates from outside of the loop inward and is only coincidentally at the position of the β = 1 layer for this screenshot. This propagation inward is only apparent (Raes et al. 2017), since it results from a cone-shaped area of high c_A propagating upward. Those waves are probably Alfvén waves that are excited by the first impulse of the driver. Other than that, stationary counter-streaming regions appear around the β = 1 layer with the same periodicity as the normal component.

To help visualize the direction of the velocity perturbations we overplot the color scale for the parallel component with vectors showing the horizontal velocity perturbation (Fig. 8) for the same time and height as in Fig. 7. It is now obvious that the plasma expands from the loop centers toward the centers of the loop exteriors (and half of a period later vice versa). Where the loops are closest to each other, counter-streaming flows get deflected sideways toward the centers of the loop exteriors. This also solidifies our interpretation of a deformed sausage wave as mode identification.

Fig. 8. Vectors of the horizontal velocity perturbations at 2 Mm at time t = 106 s. The color scale shows the longitudinal velocity component in m s−1 and the black contours show the β = 1 border. The spatial scales are given in Mm.

In order to study the wave propagation through the system, we looked at the data of two vertical lines in the domain: one of the lines is located in the loop interior where there is a low magnetic field, while the other is close to the center of the loop exterior, which has the maximum magnetic field. The vertical line in the loop interior is at a distance of 0.69R from the loop center and therefore at the same position in our MHS atmosphere as the ratio of displacements in Fig. 6 (left) in the simplified model of Sect. 4.2. We refrain from using data from the exact center of the loops and loop exteriors because the azimuthal unity vector ${\widehat{\mathit{e}}}_{\mathrm{a}}$, and therefore also for the normal unity vector ${\widehat{\mathit{e}}}_{\perp }$, are not defined there. The positions of the vertical lines are shown in Fig. 9, while the velocity components in these lines are plotted as a function of time in Fig. 10. In the latter figure, four characteristic speeds are indicated: the local sound speed (dashed black line), local Alfvén speed (dotted black line), local tube speed (solid black line), and kink speed (dashed dotted line). As seen in Eq. (21), the kink speed is calculated by external and internal Alfvén speed and density. Those values are the mean of the values inside and outside the loop for each height with the border at β = 1. Since the loop width is very constant above the transition region, the kink speed is only plotted from the transition region upward.

Fig. 9. Location of the vertical lines of Fig. 10. The black lines show the β = 1 contour at a height of 2 Mm, where there is high β inside and low β outside the “circles”. The red and blue crosses denote the location of the vertical line in the loop interior (Fig. 10 left) and the loop exterior (Fig. 10 right), respectively.

Fig. 10. Time-distance plot for the longitudinal (top), normal (middle), and azimuthal (bottom) velocity component for a vertical line in the loop interior (left) and loop exterior (right). The dotted red lines denote the transition region and the solid red line shows the β = 1 height with high β beneath the line. Plots without that line have high β for all heights. The dashed black lines show the sound speed, dotted black lines show the Alfvén speed, dash dotted black lines show the kink speed, and solid black lines show the tube speed. The extreme values in the middle right plot are saturated to make the other structures visible as well.

The plots on the left side of Fig. 10 show the velocity components in the loop interior, where there is high β for all heights. For the longitudinal component inside the loop (top left) the waves propagate smoothly with approximately the sound speed or kink speed and do not seem to be disturbed by the transition region (region between red dotted lines), except for some slight reflection, which is visible between 110 and 130 s. No features propagating with slower speeds are visible. However, if we look onto the same component but in the loop exterior (Fig. 10 top right), we see a different picture. There, the waves have to travel through the β = 1 layer before going through the transition region. Below this border, where β ≫ 1, the waves travel again with the sound speed as for the loop interior, but as soon as the first waves pass the β = 1 layer there are suddenly wave features that travel with the Alfvén or tube speed. Following the definition for mode conversion and transition of Cally (2005), this could be interpreted as a conversion from fast acoustic waves to fast magnetic waves. However, we have to be cautious when calling a wave acoustic or magnetic above the transition region in our model because we are dealing with tube waves with high β inside the loops and low β outside the loops. In addition, we also see features propagating with the sound speed or tube speed above the β = 1 layer, which seems like a transition from fast acoustic waves to slow acoustic waves at first sight. Similar to the waves inside the loop, there is a sign of reflection from the transition region, which is best visible at a time of about 135 s. However, we also have unwanted reflections from the upper boundary, which distort the wave shapes coming from below.

Compared to the longitudinal velocity components, the normal velocity components (Fig. 10 middle) are about one order of magnitude smaller. This is no surprise, as only the longitudinal component is driven at the bottom of the domain, while the other components arise from mode conversion or coupling due to inhomogeneity. The general behavior is similar to the longitudinal component, where there are only waves propagating with the sound speed or kink speed in the loop interior and a combination of slow and fast waves in the loop exterior. However, apart from the much stronger reflections from the transition region compared to the general amplitude and the much less prominent reflections from the upper boundary, two striking effects appear. The first is the high velocity amplitude around the lower border of the transition region, which we do not investigate in this paper. From the first high amplitude wave, a wave is launched that travels faster than the fastest characteristic speed (i.e., sound speed for the loop interior and Alfvén speed for the loop exterior) until it vanishes. This could be a sign of a leaky sausage wave, as these waves can travel faster than the external Alfvén speed (Pascoe et al. 2007), which would be the Alfvén speed in the center of the loop exterior for our case. By looking closer at the simulation data we find indications that this is indeed the case. We did not find a leaky sausage mode with high phase speed for our simplified model in Sect. 4, however, it could still appear in our more complicated MHS model. The second effect are the vertical stripe patterns, which are more pronounced inside the loop than outside. These may either be interference patterns by partial reflection of the waves due to the temperature gradient (the vertical temperature gradient is stronger inside the loops than outside), or just artifacts due to the cylindrical structure within a Cartesian grid. Those patterns are not of interest for our study and are therefore not considered in the following.

The azimuthal velocity component (Fig. 10 bottom) closely resembles the normal component, except for the generally smaller amplitude and that we now also see waves with a phase speed of approximately the Alfvén speed inside the loop. In addition, this component is more strongly affected by the reflection from the upper boundary. The close relation between the normal and azimuthal component is no surprise, as Fig. 8 shows us that the azimuthal velocities arise from the deflection of expanding (normal) plasma movements toward the centers of the loop exteriors.

To roughly estimate how much of the wave energy is reflected at the transition region in Fig. 10, we determine the wave energy flux parallel to the magnetic field and look at the ratio of maximum (positive, upward) to minimum (negative, downward) amplitude of the first wave front. This shows that about 3% of the flux is reflected in the loop interior and about 4% in the loop exterior. Because of this crude approximation there might be an error of several percent, but the reflected flux is still very low. For a steeper transition region, more reflection would be expected. We refrain from giving an estimate about the amount of flux reflected from the upper boundary, since this is a purely numerical issue that holds no physical value.

From Fig. 7 we get an approximate loop radius of R ≈ 0.5 Mm and from Fig. 10 we can estimate that the wavelength lies between 5 and 10 Mm. From those values we determine that kR lies between 0.31 and 0.63, which allows us to compare our simulation data to the wave mode solutions of the simplified model in Figs. 5 and 6 for small kR. All our data suggests that the parallel displacement is larger than the perpendicular displacement, which immediately excludes the fast kink surface mode. We also exclude the slow kink body mode because there is no reason why a kink mode would be excited owing to the symmetry of our system. We therefore conclude that the fast waves we see in Fig. 10 are fast sausage surface waves, which is supported by the fact that the magnitude of the velocity perturbations is smaller in the loop center than in its surroundings, while the slow waves are slow sausage body waves. However, for small kR there is no solution of fast sausage surface waves (Fig. 5), so the simple model we used to calculate the phase speed diagram does not describe those wave modes well in that region. In addition, we might still have a plane-like wave traveling through the whole domain, which mostly ignores the horizontal structuring because of its relatively slow changes. It is, however, difficult to distinguish such waves from the fast surface sausage waves.

5.2. Wave propagation for inclined flux tubes

We now study the case with the inclined magnetic field for an inclination of θ = 15°. The time development of a horizontal cut of the simulation box with 15° inclined magnetic field at a height of 2 Mm is shown in the movie that is available online. A snapshot of this movie at t = 112 s is shown in Fig. 11. The parallel velocity component in Fig. 11 behaves very similarly to the vertical case (Fig. 7), where there are higher velocity perturbations within the loops than outside, which change signs over time. However, both the normal and azimuthal components show not only a similar stripe pattern, but also have the same magnitude. These effects arise from the whole plasma moving first to the left (like in the figure) and later to the right¹, which corresponds to a kink wave. During the transition from plasma moving from one direction to the other, we have a short time span in which plasma flows both to the left and right. Another kind of fast wave was also excited because we expect the magnitude of perpendicular displacement to be larger than the magnitude of parallel displacement for fast kink surface waves (see Fig. 6); however, we find the opposite in Fig. 11. This wave is either a fast sausage surface wave, such as in the case with vertical flux tubes or a plane wave that does not “feel” the loop structure.

Fig. 11. Components of the velocity perturbation in a horizontal cut at 2 Mm for θ = 15° at a time of 112 s after the start of the simulation. The velocities are given in m s−1 and the spatial scales are in Mm. The black pseudo-circles show the β = 1 contour, with β ≫ 1 inside and β <  1 outside the loops. The values of the first two pixels next to the margin are set to zero for the normal and azimuthal component, as the corresponding vectors were badly defined in that region. The temporal evolution is available as an online movie.

We checked the kink wave assumption by plotting the velocity disturbance perpendicular to the magnetic field in the loop center as a function of the height, where we find upward propagating waves for the inclined cases with θ = 15° and θ = 30°, whereas there is no kink wave for the vertical case (θ = 0°). These kink waves are excited in the simulations with the inclined loops, as the driver is still purely vertical and therefore gives the loops a push perpendicular to the magnetic field. Since the magnetic field is only inclined in the x-direction, there is no significant perturbation of the loop centers in the y-direction. A snapshot of these waves is shown in Fig. 12.

Fig. 12. Velocity perturbation perpendicular to the magnetic field in the loop center as a function of height at time t = 118 s for three different magnetic field inclinations.

Figure 13 shows the equivalent of Fig. 10 for the inclined case, for which we examine the wave propagation along two lines inclined 15° from the vertical (parallel to the magnetic field) that lie within the loop interior and loop exterior, respectively, at equivalent locations as for the vertical case. The black lines in Fig. 13 show the local sound speed (dashed), local Alfvén speed (dotted), local tube speed (solid), and kink speed (dashed dotted) for wave propagation along the loop, i.e., along the magnetic field. Similarly, the brown lines show these speeds for vertical wave propagation, i.e., the direction of the driver. The reason for the strange S-shape for the vertical propagation speeds is that for different heights z vertically propagating waves that reach that height also start from different horizontal locations x, where the characteristic speeds are different. Therefore, for some vertical lines, the waves would find more “favorable” conditions to propagate than for others, which causes them to arrive earlier at a larger height. The S-shape is more pronounced for the Alfvén speed (and therefore also the tube speed) because the Alfvén speed changes more between loop interior and exterior than the sound speed. Since vertically propagating waves would have to go through regions with very low Alfvén speed for the y location of the loop-interior-line (left column), the S-shape is extreme in those plots. However, we do not see any signs in Fig. 13 of waves propagating along these strange (brown) lines. Therefore, we can conclude that the waves mainly propagate along the magnetic field (black lines) from the transition region onward.

Fig. 13. Time-distance plot for the parallel (top), normal (middle), and azimuthal (bottom) velocity component for a line in the loop interior (left) and loop exterior (right) in a simulation with a magnetic field that is inclined 15°. The dotted red lines indicate the transition region and the solid red line shows the β = 1 height with high β beneath the line. The plots without that line have high β for all heights. The dashed lines show the sound speed, dotted lines show the Alfvén speed, dash dotted lines show the kink speed, and solid lines show the tube speed. The black lines show these speeds for wave propagation along the magnetic field, while the brown lines show speeds for vertical wave propagation. The extreme values (only in reflection zones) in the plots for the parallel and azimuthal component are saturated to make the other structures visible as well.

From Fig. 13 it is immediately apparent that the reflection problem from the upper boundary is much stronger for the inclined case. For this reason the plots for the parallel component (top row) and for the azimuthal component (bottom row) are saturated to allow a better visibility of the waves propagating upward. What is also noticeable is that the magnitudes of the normal and azimuthal components are much higher than before, as also seen in Fig. 11.

As in the vertical case, the parallel component propagates with the sound speed or kink speed in the loop interior for all heights and splits at the β = 1 layer into waves propagating with the sound speed or kink speed and waves propagating with the Alfvén speed or tube speed in the loop exterior. In the plots for the normal component we again see (small) wave signatures that propagate faster than the local fast speed in the transition region at about t = 180 s, which could be a sign of leaky sausage waves. These wave signatures are also present in the azimuthal component, but not visible in the displayed data range. In both the normal and azimuthal component we only see wave propagation with the sound speed or kink speed for high β and the Alfvén speed or kink speed for low β.

There is again some wave reflection from the transition region, which is about 3% of the energy flux in the loop interior and 6% in the loop exterior. Given that this is just an estimation, it is not possible to say if the reflection is dependent on the inclination angle, since these values are very similar to those obtained for the vertical case (3% and 4%).

Similar as before in the vertical case, we compare our simulation results with Figs. 5 and 6. We found from Figs. 11 and 12 that definitely a kink wave is excited, and from Fig. 13 that it propagates fast. Therefore, we identify it as a fast kink surface wave. To explain the ratios of velocity perturbation components in Fig. 11 there must also be another fast wave, which could be a fast sausage surface wave or a plane fast wave that ignores the cylindrical structuring. The slow waves appearing in Fig. 13 are symmetric around the center of the loop exterior and are therefore identified as slow sausage body waves.

To reiterate the results of this section, we inserted a fast (acoustic) wave at the bottom of our domain, which converted to different wave modes. This was also found by Cally (2017), who injected a fast (magnetic) wave into a model with gravitationally stratified Alfvén speed profile with discrete (and also “touching”), inclined flux tubes, using the cold plasma (β = 0) approximation. The initially m = n = 0 fast waves scattered in Fourier space into other modes, i.e., essentially the m = 0, n = −1 kink mode. In addition, there was also significant conversion to Alfvén waves, which decayed with height as they were also scattered into higher mode numbers as a consequence of mode mixing.

5.3. Mode conversion

We would like to estimate how much the initially acoustically dominated waves, as excited by the driver, take on magnetic properties. This can be described by mode conversion from acoustic to magnetic waves. When speaking about mode conversion in the following, we mean conversion from acoustic to magnetic behavior, not conversion from fast to slow waves or vice versa. An often used conversion coefficient was given by Eq. (26) of Schunker & Cally (2006) and we repeat for convenience, i.e.,

$$\begin{array}{c}\hfill C=1-T=1-exp{[-\frac{\pi k^2{k}_{\perp }^2}{|k_z|(k^2+k_{\perp }^2)}{\left(\frac{d(c_{\mathrm{A}}^2/c_{\mathrm{s}}^2)}{\mathrm{d}z}\right)}^{-1}]}_{c_{\mathrm{A}}=c_{\mathrm{s}}}·\end{array}$$(22)

The accuracy of this expression was analytically tested by Hansen & Cally (2009). It is defined by the portion of the wave energy flux that is converted from acoustic to magnetic, where C = 1 (transmission coefficient T = 0) describes full conversion from acoustic to magnetic waves (fast wave to fast wave) and C = 0 (T = 1) describes no conversion. In this equation, k is the wave number with its components in the z-direction k_z and the direction perpendicular to the magnetic field k_⊥. The equation is evaluated at the c_A = c_s layer, where mode conversion is supposed to happen. In the following, we check Eq. (22) in our simulation data.

We calculate the conversion coefficient at three different points at the c_A = c_s layer. All three points lie in the loop exterior below the transition region and are located at different arbitrary distances from the center of the loop exterior. Since the excited acoustic waves propagate vertically when below the transition region, we can assume that k_z = k and k_⊥ = k sin(θ). Furthermore, we assume that k = ω/c_eq, where ω is our driver frequency and c_eq is the sound or Alfvén speed at the equipartition layer. The resulting conversion coefficients are plotted in Fig. 14 (top) for various inclination angles in 5° intervals as solid red lines. Two of those lines are more similar because the points in which they were calculated lie closer to each other than to the third point. The mode conversion clearly increases with increasing angle, as we expected. At a field inclination of about 15° the curvature changes from positive to negative.

Fig. 14. Top: conversion coefficient describing mode conversion from acoustic to magnetic waves as a function of inclination angles. The full red lines are mode conversion according to Schunker & Cally (2006) (Eq. (22)) in three different points within the loop exterior at the equipartition layer. The dotted black line is calculated according to Eq. (25) at z = 2 Mm and averaged over the whole horizontal plane. The dashed black line is the same, but only considering fluxes in regions with β <  1 (see text). Middle: absolute error of the red curves (Schunker & Cally 2006) compared to the dashed curve (Eq. (25) for β <  1). Bottom: same, but relative error.

To compare the outcome of Eq. (22) with a reasonable quantity of our simulation results, we use the mean acoustic and mean magnetic energy flux defined by Bray & Loughhead (1974),

$$\begin{array}{c}\hfill {\mathit{F}}_{\mathrm{ac}}=⟨p_1{\mathit{v}}_1⟩,\end{array}$$(23)

$$\begin{array}{c}\hfill {\mathit{F}}_{\mathrm{mag}}=⟨{\mathit{B}}_1\times ({\mathit{v}}_1\times {\mathit{B}}_0)⟩/{\mu }_0,\end{array}$$(24)

where p₁, v₁, and B₁ are the perturbations of pressure, velocity, and magnetic field, respectively, and B₀ is the equilibrium magnetic field. The (outer) brackets denote the average over time, however, for the following we also average over space. We now define a quantity

$$\begin{array}{c}\hfill C_{\mathrm{f}}=\frac{|{\mathit{F}}_{\mathrm{mag}}|}{|{\mathit{F}}_{\mathrm{mag}}|+|{\mathit{F}}_{\mathrm{ac}}|},\end{array}$$(25)

which should tell us how much of the energy flux, which is initially fully acoustic, was converted into magnetic energy flux. Since we want to know how much our waves are dominated by which kind of energy flux, without taking into account the direction, we avoid positive and negative fluxes to be canceled out by taking the absolute value before averaging the fluxes. The value C_f has the same properties as C given that it is 1 for full conversion and 0 for no conversion. The dotted black line in Fig. 14 (top) shows C_f averaged in time over half of a period at a height of 2 Mm, starting from the time when the first wave reaches that height, and averaged in space over the whole horizontal plane. The dashed black line shows the same, but only averaged in space over the areas where β <  1. The latter line shows higher values, as there is more flux converted in the considered regions than in the rest of the plane because there the waves have to travel through the β = 1 (and c_A = c_s) layer. That line is therefore more comparable to mode conversion in the three chosen points than the dotted line. Like before, we see more conversion from acoustic to magnetic waves with increasing magnetic field inclination θ and a change of curvature at around θ = 15°.

By comparing C_f (Eq. (25)) with C (Schunker & Cally 2006, Eq. (22)) it is apparent that the curve of the latter conversion coefficient shows a higher inclination for increasing θ than the former. Figure 14 shows the absolute (middle) and relative (bottom) errors of C compared to C_f (for β <  1). Between inclination angles of around 10° and 30° the relative error stays within ±40%. Below 10° the relative error is due to the small values of C and C_f much higher, however, the absolute error stays within the interval [−0.1,0]. In general, both conversion coefficients show the same qualitative behavior and are much more similar than we expected. This result is interesting because in our simulations we do not just have simple plane waves, except perhaps right above the driver at low heights along with some additional plane-like waves higher up. Instead, we mainly have waves that are modified by the cylindrical structure of the atmosphere, in particular sausage waves and kink waves. The local analysis around a point at the c_A = c_s layer of Schunker & Cally (2006) allows the use of the simple analytic formula in Eq. (22) for those cases as well.

We note that another, less general conversion coefficient was given by Cally (2005). It yields the same results as in Fig. 14 for θ ≤ 10°, but deviates from Eq. (22) for higher inclination angles, as the curve does not change its curvature.

6. Summary and conclusions

In this paper we presented a simple method to calculate a 3D MHS equilibrium when a divergence-free magnetic field is given. We built an equilibrium model resembling the solar atmosphere from the photosphere to the lower corona, including four flux tubes of decreased Alfvén speed and increased sound speed with the inclination θ from the vertical. This led to a tube-like plasma-β = 1 layer with β >  1 inside the tubes and everywhere in the bottom layers, and β <  1 outside the tubes starting from a certain height. We then perturbed the plasma at the bottom with vertically polarized gravity-acoustic waves according to an analytical solution. We investigated the resulting waves by studying the behavior of the velocity perturbations parallel to the magnetic field, perpendicular to the magnetic iso-surfaces and azimuthal to these surfaces in a horizontal plane above the transition region. In addition, we studied the propagating velocity disturbances in two different lines (inside and outside the loop) along the magnetic field. By comparing our results with a simplified model of waves in a cylindrical structure, we could classify the waves appearing in our simulations. For the vertical case (θ = 0°), where the magnetic field and flux tubes are oriented along the driver polarization direction, we identified deformed fast sausage surface waves and slow sausage body waves. There might have additionally been a plane-like wave excited, which is difficult to distinguish from the fast sausage surface mode. For the inclined flux tubes and magnetic field, where the driver polarization now has a component perpendicular to the tubes, a fast kink surface wave is excited in conjunction with either a (deformed) fast sausage surface wave or a plane-like wave. Moreover, we also find slow sausage body waves.

In addition, we investigated the mode conversion from the initially acoustic waves to magnetic waves. We compared the outcome of a simple formula for a mode conversion coefficient by Schunker & Cally (2006) with the ratio of the magnetic energy flux to the sum of the magnetic and acoustic energy flux. We find that both methods give similar results with a maximum absolute error of 0.1 for inclination angles from θ = 0° to 10° and a maximum relative error of 40% for angles from θ = 10° to 30°. The deviation of the simple formula from the other method is remarkably small, given that we anticipated a large influence of the cross-field wave speed structuring. This validates the frequent use of the simple formula. We note, however, that the influence of a cutoff region was not tested in the present work.

According to our simulations, vertical gravity-acoustic waves from the photosphere are converted to waves with partial magnetic properties in areas with flux tubes (especially in between the tubes), if the magnetic field lines and the flux tubes are inclined from the vertical. In that case the initially vertically propagating plane waves changed direction to propagate along the magnetic field above the transition region and were transformed into kink and sausage modes. In the case with vertical magnetic field and flux tubes, we only observed sausage waves without any significant magnetic wave properties.

There are some important limitations of this work we would like to mention. First, because of our model containing straight flux tubes, the magnetic field does not change along the loops to satisfy ∇ ⋅ B = 0. However, this leads to an unusually high magnetic field strength of 300 G in the corona. In addition, there are regions in our model with β >  1 in the corona, which is a much higher value than expected for that height (see, e.g., Gary 2001). Realistic flux tubes are expected to strongly expand between photosphere and lower corona. Such a model would not only allow the magnetic field to decrease with height, but the expansion would also affect waves guided along the flux tubes. We assume the biggest change would be that the wave fronts are refracted along the field lines and would broaden, which would lead to damping of the waves. This was also mentioned in the results of Mumford et al. (2015). In fact, we assume the simple geometry of our flux tubes to be the main reason why the damping of the waves in our simulations is far smaller than in the observations of Grant et al. (2015). Despite these drawbacks we decided to first study straight flux tubes as they simplify the analysis of the excited waves. A similar study with expanding flux tubes is currently in progress.

A second limitation of this study are the high frequency and low amplitude of our driver. While we do not think that a higher amplitude would change the core of our results, we assume that a lower frequency of the driver with a realistic period of about five minutes leads to less waves being transmitted into the corona. This is because of the acoustic cutoff region, which prohibits the propagation of acoustic waves below the cutoff frequency. Instead, many waves would be reflected from that region. Since the wavelength in the current study is already much larger than the flux tube radius, we do not expect other big changes by decreasing the frequency. We plan to use a driver period of approximately five minutes in future work.

Movies

Movie of Fig. 7 Access here

Movie of Fig. 11 Access here

Acknowledgments

We would like to thank Paul Cally for helpful comments and suggestions. Furthermore, we would like to thank the referee for helping us to improve the manuscript. This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No. 724326).

References

All Figures

	Fig. 1. Magnetohydrostatic model of four loops for θ = 0° (top row) and θ = 15° (bottom row). Left column: magnitude of the total magnetic field. Middle column: pressure (logarithmic scaling). Right column: temperature. All length scales are in Mm.
In the text

Fig. 2. Profiles and profile ranges of the atmospheric quanities for the vertical case as a function of z. The lines indicate the minimum and maximum values at each height, respectively. The striped region indicates the transition region. Left: pressure according to the VAL-C model (orange dashed lines) and with added exponential term (red solid lines). Right: temperature (green solid lines) and density (blue dashed line, constant in x and y).

In the text

	Fig. 3. Contour of β = 1 for θ = 0° (left) and θ = 15° (right). At the bottom of the domain (photosphere) β ≫ 1.
In the text

	Fig. 4. Orthogonal decomposition vectors parallel (${\widehat{\mathit{e}}}_{\Vert }$), normal (${\widehat{\mathit{e}}}_{\perp }$), and azimuthal (${\widehat{\mathit{e}}}_{\mathrm{a}}$) to the magnetic iso-surfaces (cylindrical shapes) for a vertical (left) and inclined (right) loop.
In the text

	Fig. 5. Phase speed diagram of sausage modes and kink modes for conditions similar to the model atmosphere at a height of 2 Mm. The horizontal gray lines indicate various characteristic speeds.
In the text

	Fig. 6. Ratio of longitudinal displacement and perpendicular displacement of plasma for sausage and kink modes for conditions similar to the model atmosphere at a height of 2 Mm. Left: inside the loop at r = 0.69 R. Right: at the surface of the loop (r = R).
In the text

Fig. 7. Components of the velocity perturbation in a horizontal cut at 2 Mm for θ = 0° at a time of 106 s after the start of the simulation. The velocities are given in m s−1 and the spatial scales are in Mm. The black lines show the β = 1 contour, with β ≫ 1 inside and β <  1 outside the loops. The values of the first two pixels next to the margin are set to zero for the normal and azimuthal component, as the corresponding vectors were badly defined in that region. The temporal evolution is available as an online movie.

In the text

	Fig. 8. Vectors of the horizontal velocity perturbations at 2 Mm at time t = 106 s. The color scale shows the longitudinal velocity component in m s−1 and the black contours show the β = 1 border. The spatial scales are given in Mm.
In the text

	Fig. 9. Location of the vertical lines of Fig. 10. The black lines show the β = 1 contour at a height of 2 Mm, where there is high β inside and low β outside the “circles”. The red and blue crosses denote the location of the vertical line in the loop interior (Fig. 10 left) and the loop exterior (Fig. 10 right), respectively.
In the text

Fig. 10. Time-distance plot for the longitudinal (top), normal (middle), and azimuthal (bottom) velocity component for a vertical line in the loop interior (left) and loop exterior (right). The dotted red lines denote the transition region and the solid red line shows the β = 1 height with high β beneath the line. Plots without that line have high β for all heights. The dashed black lines show the sound speed, dotted black lines show the Alfvén speed, dash dotted black lines show the kink speed, and solid black lines show the tube speed. The extreme values in the middle right plot are saturated to make the other structures visible as well.

In the text

Fig. 11. Components of the velocity perturbation in a horizontal cut at 2 Mm for θ = 15° at a time of 112 s after the start of the simulation. The velocities are given in m s−1 and the spatial scales are in Mm. The black pseudo-circles show the β = 1 contour, with β ≫ 1 inside and β <  1 outside the loops. The values of the first two pixels next to the margin are set to zero for the normal and azimuthal component, as the corresponding vectors were badly defined in that region. The temporal evolution is available as an online movie.

In the text

	Fig. 12. Velocity perturbation perpendicular to the magnetic field in the loop center as a function of height at time t = 118 s for three different magnetic field inclinations.
In the text

Fig. 13. Time-distance plot for the parallel (top), normal (middle), and azimuthal (bottom) velocity component for a line in the loop interior (left) and loop exterior (right) in a simulation with a magnetic field that is inclined 15°. The dotted red lines indicate the transition region and the solid red line shows the β = 1 height with high β beneath the line. The plots without that line have high β for all heights. The dashed lines show the sound speed, dotted lines show the Alfvén speed, dash dotted lines show the kink speed, and solid lines show the tube speed. The black lines show these speeds for wave propagation along the magnetic field, while the brown lines show speeds for vertical wave propagation. The extreme values (only in reflection zones) in the plots for the parallel and azimuthal component are saturated to make the other structures visible as well.

In the text

Fig. 14. Top: conversion coefficient describing mode conversion from acoustic to magnetic waves as a function of inclination angles. The full red lines are mode conversion according to Schunker & Cally (2006) (Eq. (22)) in three different points within the loop exterior at the equipartition layer. The dotted black line is calculated according to Eq. (25) at z = 2 Mm and averaged over the whole horizontal plane. The dashed black line is the same, but only considering fluxes in regions with β <  1 (see text). Middle: absolute error of the red curves (Schunker & Cally 2006) compared to the dashed curve (Eq. (25) for β <  1). Bottom: same, but relative error.

In the text