© ESO 2018

1 Introduction

Observations of transition region spectral lines reveal the presence of redshifts for lines formed in the temperature range from ~10 000 K to temperatures of about 2.5 × 10⁵ K. Above that temperature one often finds that the lines are blueshifted. A large amount of literature has been devoted to understanding the observed average Doppler shifts of emission lines formed in the solar transition region and corona (Doschek et al. 1976; Mariska et al. 1978; Gebbie et al. 1981; Feldman et al. 1982; Dere et al. 1984, 1989; Klimchuk 1987; Athay & Dere 1989; Rottman et al. 1990; Brekke 1993; Brekke et al. 1997; Chae et al. 1998; Peter & Judge 1999; Peter 1999; Dadashi et al. 2011). These observations have revealed that mean line profiles of transition region lines in the chromospheric network are persistently redshifted by up to 15 km s⁻¹ in lines formed between T = 0.1 to 0.25 MK. These redshifts imply the presence of plasma flows or wave motions in the quiet Sun with amplitudes that are substantial fractions of the sound speed.

A large number of models have been proposed to explain the redshifts in the transition region, including downward propagating compressive waves generated by, for example, nanoflares in the corona (Hansteen 1993), downflows of previously heated spicular material (Pneuman & Kopp 1977; Athay & Holzer 1982; Athay 1984), rapid episodic heating at low heights of the upper chromospheric plasma to coronal temperatures (Hansteen et al. 2010), and plasma draining from reconnecting loops systems (Zacharias et al. 2011), but no definitive concensus has emerged. With the development of realistic three-dimensional radiation magnetohydrodynamics (3D rMHD) models (Gudiksen & Nordlund 2005; Hansteen et al. 2007; Bingert & Peter 2011; Gudiksen et al. 2011) over the past two decades, detailed comparison between model and observation has become feasible.

Peter et al. (2006) found that three-dimensional (3D) numerical models that span the photosphere to corona produce redshifts in transition region lines. In these models, coronal heating is due to the Joule dissipation of currents that are produced as the magnetic field is stressed and braided by photospheric motions. Hansteen et al. (2010) extended this work and showed that transition region redshifts are naturally produced in episodically heated 3D models in which the average volumetric heating scale height lies between the chromospheric pressure scale height of 200 km and the coronal scale height of 50 Mm.

However, even with full access to all relevant variables on an Eulerian grid it has proven difficult to follow the evolution of mass in the outer atmosphere, especially in the transition region where timescales are short, temperature and density gradients are very high and the magnetic topology is very complicated. Previously, this problem was tackled by injecting tracer fluids for certain temperature intervals into Bifrost simulations (Guerreiro et al. 2013) in order to follow mass motions, but naturally, this approach only covers certain regions and not the entire atmosphere, and it also proved difficult to follow the dynamics of mass backward in time. The evolution of visible structures in the solar atmosphere, such as fibrils in the chromosphere or magnetic bright points in the photosphere, has been used to follow magnetic field lines and waves over time, in simulations as well as observations (Leenaarts et al. 2015). However, Leenaarts et al. (2015) concluded that using the swaying motion of fibrils as a tracer of chromospheric transverse oscillations must be done carefully. Magnetic field lines can be traced easily in the simulations by integrating the magnetic field vector from a given grid point for a given distance. However, it is not straight forward to follow magnetic flux as it evolves in time. We have found that the implementation of passive tracer particles, so-called corks, into the Bifrost simulations offers a unique and consistent approach to follow magnetic field lines over time and address the evolution of mass flows between the various layers of the solar atmosphere.

The structure of the paper is as follows. In Sect. 2, the rMHD simulation and implementation of the passive tracer particles are described. In Sect. 3, the setup and cork analysis are explained. In Sect. 4, the results are presented. We finish with a discussion and outlook section (Sects. 5 and 6).

2 Simulations

The simulations were performed with the Bifrost stellar atmosphere code (Gudiksen et al. 2011) in a box extending from the upper convection zone to the low corona. The box size is 24 × 24 Mm² in horizontal direction and ~17 Mm in the vertical reaching 2.4 Mm below the solar surface and 14.4 Mm above. The simulation grid has a size of 504 × 504 × 496 grid points. This corresponds to a resolution of 48 × 48 km² in the horizontal direction. In the vertical direction, a nonequidistant grid spacing was applied and the resolution ranges between 19 km in the photosphere and transition region to ~100 km in the corona, where the temperature and density gradients are less steep.

The Bifrost code includes an equation of state based on solar abundances in local thermal equilibrium (LTE), optically thick radiative transfer, and losses from the photosphere and chromosphere including scattering. In the middle and upper chromosphere and in the transition region, a recipe for effectively thin radiation is used (Carlsson & Leenaarts 2012). In the corona, thermal conduction along the magnetic field plays an important role. To avoid extremely short time steps, this is treated by operator splitting, where the resulting implicit operator for conduction is solved using a multigrid method. For more details, see Gudiksen et al. (2011).

The magnetic field configuration was chosen in a way that it leads to a small network-like configuration. It was created by specifying the magnetic field at the bottom boundary using a potential field extrapolation to compute the field in the entire computational domain. The magnetic field was inserted into a relaxed quasi-steady state hydrodynamical simulation, which is allowed to evolve freely. The simulation was run for 3000 s of solar time using LTE ionization, before nonequilibrium hydrogen ionization was switched on for 830 s. Afterward, starting from t = 3850 s, the simulation was run for ~1500 s of solar time without nonequilibrium hydrogen ionization including passive tracer particles. More details are provided in Carlsson et al. (2016).

Bifrost is an explicit code with diffusive terms in the equations in order to ensure stability and to ensure that structures smaller thanthe grid size cannot develop. The diffusive operator employed is split in a small global diffusive term and a location specific hyper diffusion term (see Eq. (9) in Gudiksen et al. 2011). Spatial derivatives and the interpolation of variables are performed using high order polynomials. The equations are stepped forward in time using the explicit third-order predictor-corrector procedure described by Hyman (1979).

The simulations are set up in a way to study processes in the solar chromosphere and transition region with a magnetic field configuration characterized as “enhanced network”. The magnetic field at the bottom boundary consists of two patches of opposite polarity separated by 8 Mm (see Fig. 1, left panel) with an overall balanced flux. The average unsigned magnetic field strength in the photosphere is 48 G (5 mT). The magnetic field distribution does not change significantly during the simulation timespan. Both the top and the bottom boundaries are transparent. At the bottom boundary, the magnetic field is passively advected with no extra field fed into the computational domain. The effective temperature of the simulation is not set directly, but only indirectly by specifying the incoming entropy flux at the lower boundary. The average temperature in this model (Fig. 1, right panel) is maintained by fluid motions in the convection zone, radiative transfer in the photosphere, the balance between acoustic shocks and radiative losses in the lower chromosphere, and the balance of Joule and viscous heating, thermal conduction, and radiative losses in the upper chromosphere, transition region, and corona.

The energy equation in the Bifrost code is expressed in terms of the internal energy density per unit volume e as $$\frac{\partial e}{\partial t}=-\nabla \cdot (eu)-p\nabla \cdot u+Q_{\text{rad}}+Q_{\text{other}},$$(1)

where u is the velocity field, p is the gas pressure, Q_rad is the heating due to radiation, and Q_other is the heating due to heat conduction, viscous and Ohmic dissipation. The various terms contributing to the heating by radiation are $$Q_{\text{rad}}=Q_{\text{phot}}-L_{\text{chrom}}+Q_{\text{Ly}\alpha }+Q_{\text{EUV}},$$(2)

where Q_phot is the radiative heating from the photosphere described in Gudiksen et al. (2011), L_chrom describes losses from the chromosphere due to strong lines, Q_Lyα is heating from the Lyα line of hydrogen, and Q_EUV is the heating from extreme ultraviolet (EUV) photons, corresponding to the thin radiative losses from the transition region and corona (negative Q_EUV) absorbed in the chromosphere (positive Q_EUV). The three latter contributions are described in Carlsson & Leenaarts (2012). The Joule heating and viscous heating terms are defined as $$Q_{\text{joule}}=j^2/\sigma$$(3)

and $$Q_{\text{visc}}=2\eta {\xi }_{i,j}{\xi }_{i,j}.$$(4)

j is the current density, σ is the electric conductivity, η is the viscosity, and ${\xi }_{i,j}=\frac{1}{2}\left(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}-\frac{2}{3}{\delta }_{i,j}\frac{\partial u_m}{\partial x_m}\right)$ is the symmetric part of the viscous stress tensor. The thermal conduction term implemented in Bifrost takes the form $$Q_{\text{spitz}}=-\nabla F_c,$$(5)

with F_c = −T^5∕2∇_∥T (Spitzer 1956), where the gradient of T is taken only along the magnetic field and κ₀ is the thermal conduction coefficient. The conduction across the field is significantly smaller under the conditions present in the solar atmosphere and is smaller than the numerical diffusion, so it is ignored. The first two terms on the right-hand side of Eq. (1) describe the contributions due to the advection of energy and the compression or expansion of the plasma.

Fig. 1 Left panel: vertical magnetic field strength at height z = 0 Mm for timestep t = 5000 s. The color bar range is [–1.5 kG, +1.5 kG]. A set of magnetic field lines connecting the main magnetic polarities is overplotted in yellow in the left panel. Field line apices and cork positions of corks with T = 105 K are indicated by blue and red crosses. These field lines are further discussed in Sect. 4.4. Right panel: temperature distribution of the simulation grid points as a function of height for timestep t = 5000 s of the simulation. The average temperature at a given height is indicated by the gray line.

3 Method and analysis

As a new analysis tool, we have implemented passive tracer particles, so-called corks, into the Bifrost simulations. The tracking of the corks is based on the advection of the corks by the surrounding plasma flows. For the tracking of each cork in the simulation, the plasma flow at the respective cork position is considered for every single timestep. The position r_i of a cork labeled i evolves as $\frac{\partial r_i}{\partial t}=v(r_i)$, with v(r_i) being the plasma velocity at the location of the cork. Initially, corks are placed at every grid point of the simulation cube, which corresponds to a total of 504 × 504 × 496 = 125 991 936 individual corks that are followed over time. The snapshot cadence for writing out the basic MHD parameters during the simulation is 10 s; we chose the same cadence for writing out the cork positions. This doubles the amount of output data, however, little extra time is needed for the calculation and writing of the additional corks data. The number of corks in the box is assumed to remain constant throughout the simulation. Thus, corks that leave the simulation box on the top or bottom of the box are not lost, but they remain fixed at the same position, which is outside the domain considered for the analysis presented here. Since periodic boundary conditions are applied in the x- and y-direction, corks that leave the box on either side re-enter the box on the opposite side of the simulation domain.

The main difference between our cork method and the tracer fluid employed by Guerreiro et al. (2013) is that the latter is designed to address issues limited to certain temperature regions, whereas the corks cover the entire solar atmosphere without any constraints. The tracer fluid in the work of Guerreiro et al. (2013) is placed in a certain volume of the solar atmosphere defined by a specific temperature interval. Inside this specified temperature interval, the plasma density is set to the plasma density of each cell; outside the specified volume, the density is set to ρ × 10⁻⁵. It is therefore much more difficult to follow material in sufficient detail to be able to ask questions such as, where does material currently at 10⁵ K come from? In addition, the use of corks allows one to study the time evolution of specific magnetic field lines in much greater detail than when using tracer fluids. The applications are manifold, as this method allows the tracing of individual plasma structures as well as magnetic field lines backward and forward in time in the most accurate way.

Corks are injected at t = 3850 s, which is the snapshot that has been used as a reference snapshots in many earlier studies, (e.g., Carlsson et al. 2016; Leenaarts et al. 2015, and references therein). In particular, our time series is thus suitable for a comparison with the publicly available simulation of an enhanced network region (Carlsson et al. 2016).

In principle, many different initial setups for the test particles are possible. In this study, the focus is on mass flows throughout the entire solar atmosphere, therefore corks are placed at every single point of the simulation grid. Each cork has a unique ID, which allows us to investigate individual cork tracks and the collective behavior of the corks. In addition, various methods can be applied to follow the corks during the simulation. In this work, we simply follow the corks over time, as they are being advected by the plasma flows at the respective positions in the simulation cube.

The large number of corks that were tracked during the simulation allows us to perform statistical studies of the collective cork behavior, e.g., plasma flows throughout the entire solar atmosphere, as well as the evolution of various structures, e.g., coronal loops, including their physical parameters such as temperatures, densities, and velocities. Extracting the trajectories of the corks is straight forward in this setup.

In the following, corks with z > 0 Mm are referred to as atmospheric corks. Initially, 84% of the corks are atmospheric corks; 25 min later, only 79% of the corks are located above the surface. Approximately 18% of the atmospheric corks reach the upper boundary (z = 14.36 Mm) within 25 min, while ~6% reach below the surface. After 25 min solar time, 15% of all corks have reached the top boundary, and 4% of all corks have reached the bottom of the box (z = −2.44 Mm).

4 Results

In the following, we present results of a numerical experiment, in which we followed corks for approximately 25 min of solar time. The corks are initially placed at every grid point of the simulation box. At the time that corks are injected into the simulation (at t = 3850 s), the coronal part of the atmosphere is in a phase of vigorous heating. Thereafter, heating is balanced by energy losses through thermal conduction and radiation, and the average temperature structure remains stable during the second half of the simulation (see Fig. 1, right panel). Compared to the chromosphere and corona, the transition region fills only a very small volume and is characterized by high velocities. Because of these two facts corks do not spend much time in the transition region; hence the population of corks is low in this temperature layer. In Fig. 2, three examples of cork trajectories are shown, in which corks typical of their starting location have been followed for 800 s following (t = 4500 s). The evolution of the corks in height, temperature, sum of the Joule heating term and viscous heating term, and vertical velocity is shown on the right.

As an example of a cork that is originally at chromospheric temperatures and then becomes heated to upper transition region temperatures, we consider the cork labeled (a) in Fig. 2. This cork moves upward in the simulation box, from approximately 5 to 14 Mm, as its temperature rises from log T = 3.79 to log T = 5.79. We note that the cork traverses the lower transition region temperature range rapidly, in less than 50 s, following a strong upflow event at t = 4550 s. The rise in temperature can be associated with a series of small-scale heating events (Joule plus viscousheating, see Eqs. (3) and (4)) taking place over a timespan of ~300 s starting roughly at t = 4500 s. In addition, the cork experiences strong heating through heat conduction (Eq. (5)) as its temperature increases. This heating balances the cooling due to adiabatic expansion, while the cork is rising in the solar atmosphere. Radiative losses (Eq. (2)), on the other hand, are small and can be neglected. After having reached upper transition region temperatures, there are fewer episodic and burstily occurring Joule heating events and thus the heating rate of the cork decreases. This is to be expected in a model in which the heating is highly episodic and concentrated in certain spatial domains, e.g., in current sheets (Gudiksen & Nordlund 2005; Hansteen et al. 2015).

We find that material at lower transition region temperatures typically is either heated (as above), or cools within a short timespan out of the lower transition region temperature range. To illustrate the latter cooling case, we consider a cork that moves at a height of ~2–3 Mm during the first 500 s, while its temperature varies between 300 000 K–600 000 K (cork labeled (b)). At t = 4930 s, the cork temperature quickly drops down to temperatures below 20 000 K, while the cork moves downward in the solar atmosphere into denser plasma regions. During this period, radiative cooling (Eq. (2)) is the most efficient energy term. Preceding this cooling phase, a heating event located on the same loop (see Sect. 4.4) causes an increase in heat conduction (Eq. (5)) along the loop and a short-term temperature increase between t = 4900−4930 s. After the cooling phase, the cork remains at low temperature for the rest of the period.

The third and final example shows a cork, which resides in the upper transition region/low corona for the entire period (cork labeled (c)). During this time, the cork temperature varies between 600 000 K and 800 000 K, as it undergoes a series of episodic heating and cooling events, caused by a combination of Joule and viscous heating (Eqs. (3) and (4)) as well as adiabatic compression and expansion (see Eq. (1)), as the cork is moving up and down in the solar atmosphere. During these heating events, the heating rate rises by approximately one order of magnitude compared to the average heating rate. The cork remains between 8 and 13 Mm during the entire period and is a typical representative of the coronal and upper transition region corks in this simulation.

In the following subsections, we present the spatial distribution of the corks (Sect. 4.1) and their temperature distribution (Sect. 4.2) in more detail. The analysis of mass flows between different temperature layers and an explanation of the observed transition region redshifts is presented in Sect. 4.3. A scenario explaining the origin of these redshift for a particular set of magnetic field lines is outlined in Sect. 4.4.

Fig. 2 Left panel: cork trajectories of three representative examples described in Sect. 4.1. The temporal evolution for 800 s is shown. As time evolves, the colors change from blue to red. The respective projections in the xy-, xz- and yz-plane are indicated in different shades of gray for the different corks, which have been labeled (a), (b), and (c). Right panel: shown from top to bottom for the three examples on the left are the temporal evolution of cork height, cork temperature, vertical velocity, and heating and cooling terms per particle, i.e., the sum of Joule heating and viscous heating (black), adiabatic compression/expansion (red), Spitzer heat conduction (blue) and radiative losses (green). In the velocity panel, negative velocities are pointing upward, positive velocities are pointing downward. The first column shows the evolution of cork (a), which is heated from chromospheric to upper transition region temperatures starting at ~ t = 4600 s; the second column shows the evolution of cork (b), which quickly cools from transition region to chromospheric temperatures at ~t = 4900 s; and the third column shows the evolution of cork (c), which remains in the upper transition region throughout the entire simulation.

4.1 Overall distribution of the corks

Figure 3 shows a snapshot of the (x, y, z)-positions of the corks 15 min after their injection. The corks were initially placed in the xz-plane at y = 12 Mm. A movie of the entiretime series can be found in the online material. Every single dot indicates the position of one cork. Dark regions indicate regions of high cork density; white regions show regions where no corks are present. Very quickly, the corks spread throughout the simulation box and start outlining different morphological structures, including the photospheric granulation pattern close to the z = 0 Mm plane, large loop structures connecting the two main magnetic polarities located at approximately x = 7 Mm and x = 15 Mm, and short, low-lying loop-like structures. After a few minutes, much fewer corks are found above the magnetic field concentrations compared to more quiet magnetic regions (see Figs. 3 and 8, left panel). These voids are due to flow divergence and strong plasma flows above the strong magnetic field regions. In the following, we focus mostly on the vertical motion of the corks, but one should note that the corks are moving in the horizontal directions as well, as indicated by the projections of the trajectories in Fig. 2.

As the simulation progresses, corks rapidly distribute in height throughout the box. Depending on the timestep, it takes between 3 min and 6 min until the corks with an initial height of 5 Mm have traveled all the way up to the top of the box, about 5–8 min until corks with initial height of 3 Mm reach the top and 15–21 min for corks with an initial height of 1.5 Mm to reach that high. On the other hand, none of the corks with an initial height of 5 Mm ever reach below the surface within the time sequence of 25 min we analyzed. It takes between 20 min or longer for the corks with an initial height of 3 Mm to reach the surface and 7–19 min for the corks with an initial height of 1.5 Mm. Thus, the downward moving corks take much longer to reach the surface than the upward moving corks take to reach the top.

Fig. 3 Spatial distribution of corks, which have initially been placed in the y = 12 Mm plane of the simulation box, 15 min into the simulation. Every single dot indicates the position of one cork. Dark regions indicate regions, where the cork density is high; white regions represent voids of corks, i.e., regions where no corks are present. As time evolves, the corks start moving in both horizontal and vertical direction, thus slowly moving out of the y = 12 Mm plane (see Sect. 4.1). A movie showing the evolution over 25 min is available online.

4.2 Temperature evolution of transition region corks

The evolution of mass flows is most complicated for mass at transition region temperatures, hence to clarify the mass cycle for this gas we concentrate explicitly on the details of this temperature range. Perhaps the most striking development when looking at the evolution of the temperature distribution of the corks is the large drop in average temperature over large parts of the transition region temperature range within a short time. In Fig. 4, probability density functions (PDFs) of cork temperatures are presented for t = 4940 s (left panel) and t = 5060 s (right panel). The average temperature is shown by the red line, and values of the reference timestep t = 5000 s are indicated by a turquoise line. Around T = 100 000 K, the drop in average temperature within 2 min amounts to an order of magnitude.

Following a similar approach as Guerreiro et al. (2013), we investigated the temperature evolution of corks at transition region temperatures in more detail by following all corks in the temperature range log (T∕[K]) = [4.9, 5.1] over time for several minutes. Because our cork-based approach allows for both backward and forward tracing, we can thus also draw conclusions on the situation prior to an observed state. Figure 5 shows the temperature, density, and height evolution of corks in the temperature range log(T∕[K]) = [4.9, 5.1] at t = 5000 s for ±30 s (we have considered several other timespans with similar results). We find that after 1 min these corks have cooled down from an average temperatureof ~150 000 to ~ 50 000 K. At the same time, the average mass density of the corks increases by roughly half an order of magnitude and the average height decreases by ~ 0.5 Mm indicating that, on average, the corks are sinking lower into the atmosphere. The presence of global oscillations in the simulation complicates this analysis for certain timesteps, but on average, the same trend is observed for the entire time series.

After 2 min of solar time, only 4% of the corks are still in the same temperature range (log T = [4.9, 5.1]), whereas 86% have cooled to lower temperatures and 64% all the way down to below 20 000 K (see Table 1). At the same time, only 10% of the corks have been heated, i.e., around 4%, to temperatures above 300 000 K. Only very few corks, around 0.02%, reach temperatures above log(T∕[K]) = 6.0. Similar numbers are found for various other timesteps, suggesting a fast cooling process for the transition region corks once they reach the T = 100 000 K temperature range. Our results are in agreement with those of Guerreiro et al. (2013) who injected a tracer fluid filling the volume between 20 000 K and 300 000 K into the Bifrost simulations. After 2 min of solar time, they found that 95% of the material has left this region, about 70% of which has cooled to chromospheric temperatures, while 25% of the material was heated to higher upper transition region or coronal temperatures. The remaining 5% of the material remained within the initial temperature range.

On the other hand, in contrast to what was speculated by Guerreiro et al. (2013), we find no evidence for rapid heating to coronal temperatures shortly before the cooling phase sets in. Two minutes prior to t = 5000 s, the vast majority of the corks (i.e., 85%) already exist at temperatures above 100 000 K, while 4% are within the same temperature interval, whereas only 11% have lower temperature. Indeed, most of the corks already have temperatures higher than 100 000 K as early as 10 min before they transit the T = 100 000 K interval.

Fig. 4 Probability density functions (PDFs) of cork temperature. The PDFs are shown with respect to t = 5000 s for t = 4940 s (left panel) and t = 5060 s (right panel; see Sect. 4.2). The distribution of corks at t = 5000 s is shown by a turquoise line in each panel. The average cork temperature at each timestep is indicated by the red line. The blue lines show the minimum and maximum temperature of all corks for the respective height.

Fig. 5 One-minute evolution of temperature, density, and height distributions of corks in the temperature range log (T∕[K]) = [4.9, 5.1] at t = 5000 s. The evolution between t – 30 s and t + 30 s is shown by the different colors. The respective mean values are indicated by dashed lines.

4.3 Massflows in the solar atmosphere

The mass assigned to each cork is determined from the particle mass density times the volume of the cell that the cork is in at t = 3850 s. Throughout the entire simulation, the mass represented by each cork is assumed to be constant. Because of the continuous flows between the lower and upper layers of the solar atmosphere, the number of corks in each temperatureinterval varies with time, and the mass within each temperature layer is changing constantly. In order to disentangle the flows between different atmospheric layers, the corks were divided into various log T intervals ranging from log(T∕[K]) = 3.3 to 6.3 at intervals of 0.2 dex in logT. Owing to the density stratification of the solar atmosphere, corks with high masses are overwhelmingly found at low temperatures (see Fig. 6, upper left panel).

In the following, the mass evolution of corks in different temperature bins is investigated. Results are presented inFig. 6 and described in more detail in the following. Only the contributions to each selected bin are considered and all panels have been normalized to 100. As a representative timestep, t = 5000 s is chosen, and the mass evolution of the corks is investigated for ±300 s. The arrangement of the panels is of ascending order in temperature.

Most plasma at low transition region temperatures is shown to be episodically cycling between chromospheric and transition region temperatures. A small fraction of the mass (~3%) is found to have its origin at upper transition region temperatures, cooling to low transition region temperatures on a timescale of a few minutes; similarly a small fraction (~2%) is reheated to upper transition region temperatures within five minutes.

As we move to higher temperature channels, more temperature channels contribute in terms of mass. The contribution from lower temperature channels dominates, however, the fraction of hot mass that is cooling and cool mass that is being heated to upper transition region temperatures increases constantly. In the low and middle transition region (Fig. 6, second row), mass exchange takes place on relatively short timescales and the contributions of the various temperature channels are rather symmetric.

In the middle transition region temperature intervals log T = [4.9, 5.1], [5.1, 5.3], and [5.3, 5.5] (Fig. 6, third row), the asymmetry of the lines indicates that the fraction of hot mass that cools is higher than the fraction of cool mass that gets heated. This imbalance grows with higher temperature. On average, mass is cooling as it transits these bins. The temperature interval log T = [5.3, 5.5] is the temperature range, where a significant amount of material is cooling. On timescales of 10 min, a large percentage of material is found to cool all the way from log T = [5.7, 5.9] to log T = [3.9, l4.1], after passing through the logT = [5.3, 5.5] temperature band.

As we move up higher in the solar atmosphere, the distributions of cork masses become symmetric again, and the heating and cooling processes happen on longer timescales of several minutes. Most mass at upper transition region temperatures, i.e., logT = [5.7, 5.9] is already found in the same temperature interval several minutes earlier, and remains at the same temperature for several minutes after. The other half of the material has originally been heated from the temperature channel just below and some 10% of the mass has cooled from coronal to upper transition region temperatures within 5 min. A small amount of hot plasma (~2%) cools to temperatures of logT = [3.9, 4.3] within 5 min.

In the corona, about half of the mass in the temperature range log T = [5.9, 6.1] is heated from logT = [5.7, 5.9] on timescales of a few minutes. The other half of the material is found to remain at coronal temperatures for several minutes before and after. Slightly less than half of the mass cools to the log T = [5.7, 5.9] interval afterward and about half of the mass remains at coronal temperatures, however, a fraction of hot material (~8%) is found to cool all the way down to chromospheric and lower transition region temperatures (log T = [3.9, 4.3]).

In the following, density-squared-weighted vertical velocities are considered for the various temperature layers discussed above. These velocities are of interest, since they are proportional to the intensity of optically thin emission lines and thus to the observed Doppler shifts (ρ² ∝ I ∝ v_D). To reduce the effects of global oscillations present in the simulation, we analyze temporal averages of the vertical velocities for the timesteps t > 4500 s, when the simulation has settled down. The chromosphere is permeated by (mainly) upwardly propagating acoustic/slow mode waves and shocks. Therefore, we rather concentrate on average flows in the transition region and corona, i.e., on plasma with temperatures above log T= 3.9.

In the right panel of Fig. 7, the density-squared-weighted vertical velocity (black line) is shown together with the average vertical velocity (light gray line) as a function of temperature. On the left, the respective mass fractions of upward (blue) and downward (red) moving plasma are shown (solid lines) together with the respective fractions of surface area (dashed lines). This analysis is based on the actual MHD variables to avoid uncertainties due to uneven cork distributions. However, a comparison between corks and MHD variables shows that similar results are obtained for both. On average, more mass is upflowing than downflowing in the upper chromosphere, i.e., in the temperature interval log T = [3.9, 4.1]. This large percentage of mass is due to spatially localized upflows of high density plasma. Even though the surface area containing downflows is larger than surface area containing upflows, and the average vertical velocity in this temperaturerange is slightly pointing downward, these high density upflows result in observed blueshifts in the model.

Throughout the low and middle transition region, i.e., log T = [4.1, 5.5], the fractionof downflowing mass outweighs the fraction of upflowing mass by ~20%, and the surface area with downflows is larger than the surface area with upflows in the temperature range log T = [4.1, 5.3]. The averagevertical velocity is pointing downward, and since the average density of the plasma that is downflowing is higher than the average density of the plasma that is upflowing, strong average redshifts are observed. A similar situation (higher percentage of downflowing than upflowing mass) is found in the corona (log T = [5.9, 6.3]). On the other hand, in the upper transition region (logT = [5.5, 5.9]), the situation is reversed. The fraction of upflowing mass exceeds the fraction of downflowing mass, and the surface area covered by upflows is larger than the surface area covered by downflows. However, the average density is lower in regions with upflows than in regions with downflows. Therefore, the strong average upflows observed at upper transition region temperatures do not carry as much weight, and thus the average blueshift (∝ ρ²) observed in this model is notably smaller than the average upflow velocity.

In summary, the average density-squared weighted velocity, which reflects the trend of the average Doppler shifts with line formation temperatures, is pointing downward in the lower and middle transition region, because the fraction of downflowing mass outweighs the fraction of upflowing mass. This explains the observed transition region redshifts in the model. The situation is reversed in the uppertransition region, where we find that more material is upflowing than downflowing. Furthermore, in addition to having a larger mass fraction, downflows at transition region temperatures are found to last longer than upflows. For example, at t = 5000 s, 25% of the corks at logT = [4.9, 5.1] show downflows for 1 min or longer, while only 2% show upflows for 1 min or longer. On average, ~20% of the corks show downflows for 1 min or longer, while only 7% show upflows for 1 min or more. However, we note that material that is upflowing or downflowing for a minute or longer, no longer contributes to the log T = [4.9, 5.1] interval, since it is likely to be hotter or cooler after some 30 s.

Fig. 6 Upper left panel: temporal evolution of cork masses in a given temperature bin for a 10-min timespan around t = 5000 s. Only atmospheric corks are considered. The temperature bins are divided into temperature intervals of log T = 0.2 dex, and each bin is represented by a different color (see legend in panel labeled [4.9, 5.1]). The fraction with respect to the total mass of all atmospheric corks is shown (see Sect. 4.3). Rows 1 (middle and right panel)–4: temporal evolution (t  ±   300 s) of the mass of corks within a given temperature bin. Only corks that contribute to the selected bin are considered, and all panels are normalized to 100. The same legend applies for all panels.

Fig. 7 Left panel: mass fraction (solid lines) and surface area (dashed lines) of upflowing (blue) and downflowing (red) plasma in various temperature intervals in the range log T = [3.5, 6.1]. Temporal averages based on the timesteps between t = 4500 s and t = 5300 s are presented. Right panel: average vertical velocity (light gray line) and density-squared-weighted vertical velocity (black line) as a function of temperature. The latter represents the trend of the average Doppler shift with line formation temperature. The error bars show the mean error.

4.4 Flows along magnetic field lines

So far, we only considered vertical (mass) flows. Since plasma motions and thermal conduction are constrained to flow along the magnetic field in a low-β plasma, we focus on plasma flows along magnetic field lines in this section.

We again concentrate on typical corks, which are initially at mid-transition region temperatures, as the dynamics of these corks has proven to be the most complicated. Magnetic field lines have been traced through all 27 246 corks in the temperature range log(T∕[K]) = [4.99, 5.01] at t = 5000 s. A very small fraction, 0.2%, of those field lines do not connect back to the photosphere and have only one footpoint rooted in the photosphere. These open field lines are located in strong network regions.

Roughly 70% of the corks at T = 100 000 K show downflows, ~30% show upflows, and the average vertical velocity is found to be ⟨u_z⟩ = 6.7 km s⁻¹. Their spatial distribution is shown in the left panel of Fig. 8; the color coding represents the vertical velocity of each cork. The cork velocity map bears close resemblance to the C iv Doppler shift map, in regions where corks are present. A synthesized Doppler shift map of the C iv (1548 Å) line formed at T ~ 100 000 K is shown in Fig. 8 (right panel). The average Doppler velocity is found to be 8.9 km s⁻¹ at t = 5000 s. In the strong network regions, where strong upflows dominate, very few corks with T = 100 000 K are found.

The cork heights are distributed between 1.3 Mm and 7.8 Mm with a mean of 4.1 Mm. The loop apices are distributed between 2.5 Mm and 13.2 Mm with a mean of 4.8 Mm. Only about 1% of the corks lie on field lines with apex heights greater than 8 Mm. Both cork heights and apex heights are slightly lower for the downflowing corks as compared to those that are upflowing. Most field lines are roughly aligned in x-direction, which is dueto the underlying magnetic field configuration (see Fig. 1, left panel).

In the following, we focus on a bundle of magnetic field lines connected to corks with T = 100 000 K, which are showing relatively strong downflows of 15–18 km s⁻¹. In Fig. 1 (left panel) these field lines are indicated in yellow. The positions of the respective corks and field line apices are indicated by blue and red crosses. The structure we are considering is a rather cool, low-lying loop with an apex height of 4.3 Mm and an apex temperature of ~500 000 K. Figure 9 shows an example of such a magnetic field line. The top row shows projections of the magnetic field line on the xz-, yz-, and xy-planes for seven consecutive timesteps, i.e., t = 5000 ± 30 s. The evolution of various field line parameters, including temperature, pressure, pressure gradients, velocity components, and various heating and cooling terms, is shown below. The magnetic field line does not change much in shape or position over time, and the cork itself is found to be sliding down along the leg of the field line. While moving downward, the cork is cooling from approximately T = 200 000 K to T = 50 000 K.

The central part of the loop has been strongly heated. Since the loop is fairly hot, this heat rapidly spreads throughout the parts of the loop that are hot enough to support significant thermal conduction, leading to a high pressure plug in the upper part of the loop and low pressure dips in regions that have not been heated yet. The details of the two different legs depend on the previous loop history, i.e., the initial velocity and density structure, which easily leads to asymmetrical behavior in the two loop legs.

For a more detailed analysis of the (direction of the) flows acting upon the corks, a decomposition of the velocity vector was performed into velocity components parallel (u_∥) and perpendicular to the magnetic field vector (i.e., u_N in the plane of the field line and u_P perpendicular to the plane of the field line) at each position along the magnetic field line. These components are given by u_∥ = u ⋅ b, u_N = u ⋅ N, and u_P = u ⋅ P, where $T=\frac{dr}{ds}=b$, $N=\frac{dT}{ds}/\left|\frac{dT}{ds}\right|$, and P = T × N are the unit vectors for a curve r(s) in 3D space as a function of arc length s. In this coordinate system, the vertical velocity u_z can be expressed as the sum of u_∥,z + u_P,z + u_N,z = u_∥e_∥e_z + u_Pe_Pe_z + u_Ne_Ne_z. From a comparison of the parallel and perpendicular velocity components, we can obtain information on the kind of waves and disturbances that are present in the solar atmosphere.

At the position of the cork, the parallel velocity component is found to be stronger than the perpendicular velocity components. For t = 5000 s, we find u_∥,z = 12.5 km s⁻¹, u_N,z = 2.7 km s⁻¹, u_P,z = 0.74 km s⁻¹, which add up to a total vertical velocity of u_z = 16.0 km s⁻¹. In this case, the position of maximum vertical velocity in the loop is very close to (slightly below) the cork position.

A closer look at the pressure (gradient) along the field line reveals a higher pressure above the cork than below. The evolution of the plasma pressure shows that the cork is pushed down by a downward propagating pressure front from a high to a low pressure region. This process indicates a correlation between the pressure gradient and the velocity along the loop, where high downflows (redshifts) occur in regions with a large downward directed pressure gradient. Such a correlation is observed for the entire magnetic field line bundle shown in Fig. 1 (see Fig. 10). Because of their high density, such low-lying loops are associated with strong radiative losses (~ ρ²). This contributes to the cooling of the corks, as they are pushed down along the field line by a pressure wave from above. The cooling from higher to lower temperatures is typical for corks at transition region temperatures, as described in Sect. 4.2. It is not only observed for all corks of the field line bundle shown in Fig. 1, but for many other examples that we investigated. This process is found frequently both in this and other simulations (e.g., Hansteen et al. 2010; Guerreiro et al. 2013).

Fig. 8 Left: xy-positions of corks with logT = [4.99, 5.01] at t = 5000 s, color-coded with vertical cork velocity. Red indicates downflows; blue shows upflows. Right: C iv (1548 Å) Doppler shift map for t = 5000 s. Blue indicates blueshifts (upflows), red shows redshifts (downflows). The color bar on the right shows the velocity scales; the same range has been chosen for both panels.

Fig. 9 Evolution of field line parameters extracted for a low-lying field line including a downward moving transition region cork. Top row: projections of the field line on the xz-, yz-, and xy-planes. The field line parameters are shown as a function of the coordinate s along the field line, where s = 0 Mm corresponds to the position of the right footpoint of the loop in the panels of the top row. The evolution of the field line parameters is shown for t = 5000  ± 30 s, and the field line evolves from dark blue colors to red colors. The cork positions are indicated by asterisks and the loop apices by crosses.

Fig. 10 Correlation between maximum pressure gradient along the loop and u∥,z for a set of downward moving corks on low-lying field lines. The set of magnetic field lines is shown together with the vertical magnetic field strength at height z = 0 Mm in Fig. 1 (left panel).

5 Discussion and conclusions

We have presented the first quantitative study on passive tracer particles, so-called corks, which have been implemented in the Bifrost stellar atmosphere code. In particular, this novel approach sheds new light into mass motions and the formation of redshifts and blueshifts in the solar transition region. Our results show that on timescales of approximately 1 min, the majority of corks in the low and middle transition region are cooling and moving downward, while corks in the upper transition region are barely changing temperature and are mostly moving upward. The dominance of redshifts in the model can be explained as a result of a higher percentage of downflowing mass compared to upflowing mass, even though the average upflow velocity can be larger than the average downflow velocity in the lower and middle transition region. Blueshifts in the upper transition region can be explained by a combination of both, strong upflows and a higher percentage of upflowing mass compared to downflowing mass in this temperature region.

It has been shown that at any given time more than half of the corks in the transition region cool rapidly, i.e., within a minute, from temperatures above 100 000 K to temperatures below 20 000 K. Our results are in agreement with the findings of Guerreiro et al. (2013) who used tracer fluids at 20 000 K and 300 000 K in the Bifrost simulations to study mass flows. Whereas Guerreiro et al. (2013) speculated that there is rapid heating of the transition region material shortly before it undergoes cooling, our cork based approach has shown that most of the material already has resided at temperatures above 100 000 K for several minutes before undergoing cooling.

In this simulation, the cooling process is found to last longer than the heating process and downflows are found to last longer than upflows. The most commonly found situation in the transition region around 100 000 K is that corks are cooling and traveling downward. For this set of corks, a strong correlation between the pressure gradient along magnetic field lines and the velocity is found, indicating that a pressure-wave related mechanism is responsible for the observed downflows, i.e., material is being pushed downward at high velocity (~0.3 c_s, where c_s is the speed of sound).

In order to prevent the solar corona from draining, the prevailing downflows in the low and middle transition region must be balanced by upflowing material. This kind of material has been rather difficult to observe. Even less understood is how this material is heated to upper transition region and coronal temperatures on its way upward. As described in Sect. 4.3, our findings indicate that the heating of material from lower to upper transition region and coronal temperatures is a slow process, which takes several minutes. Approximately 4% of the mass in the temperature interval T = 20 000−50 000 K reaches temperatures above 500 000 K within 5 min; 0.5% reaches coronal temperatures above 800 000 K within the same time. Over a large range of heights, the heating of corks from low to upper transition region temperatures is comparable to the average heating rate in the simulation at the respective height. After having reached upper transition region temperatures, corks tend to remain in this temperature region for several minutes, thus constituting a continuous mass reservoir for the corona.

Previous studies by Hansteen et al. (2010) and Guerreiro et al. (2013) suggested a common formation mechanism for both observed redshifts and blueshifts based on rapid heating scenarios leading to a high pressure plug of material at upper transition region temperatures that relaxes toward equilibrium by expelling material. Pressure wave driven downflows along the field lines as observed in this work complement the picture: the larger the downward directed pressure gradient, the higher the observed downflow along the loop (see Fig. 10). The scenario described above (Sect. 4.4) is a process found frequently both in this simulation and in others. It is a natural consequence of strong heating in a certain region leading to overpressure, forcing material down one or both loop legs.

The cork trajectories were used to follow magnetic field lines over time in a highly accurate way. The full exploitation of this approach will be pursued in a follow-up paper, which will focus on the various kinds of loops present in the simulation. In this study, the method was applied to all corks with temperatures of 100 000 K at t = 5000 s. The field lines, which were traced for the transition region corks, have an average apex height of some 5 Mm, indicating that almost all transition region material is found in low-lying and relatively cool loops. Treating upflowing and downflowing corks separately did not reveal any differences with respect to the loop height. Using a decomposition of the velocity vector into parallel and perpendicular components along the magnetic field line, we were able to show that for downflowing corks the velocity component parallel to the field is larger than the perpendicular velocity component.

The validity of this approach is limited by the magnetic Reynolds number R_m, i.e., the ratio of the advective and diffusive term of the induction equation, which determines how well the assumption of frozen-in magnetic flux holds. Bifrost employs a space- and time-variable nonisotropic magnetic diffusion coefficient η, and we cannot therefore give a single value for R_m. Leenaarts et al. (2015) showed that the frozen-in condition holds for a large portion of the simulation domain, however, there are instances in which the field lines evolve independently of the plasma, indicating areas of strong diffusion or reconnection. In this work, we have ignored such field lines.

6 Outlook

In this study, we have focused on corks that travel into and out of the solar transition region. The corks were implemented into a so-called enhanced network simulation, but many different magnetic field configurations, including flux emergence models, are possible. Their analysis is currently underway.

We have given a description and first insights into a new method for tracing magnetic flux, which is based on the cork trajectories. The full exploitation of this approach and the analysis of different sets of field lines will be subject of another paper. Future workwill also include a more detailed discussion on the various kinds of waves based on the analysis of corks in the Bifrost models.

Movie

Movie of Fig. 3 Access here

Acknowledgements

The authors thank the anonymous referee for his or her thorough review and comments, which contributed to improving the manuscript. The research leading to these results has received funding from the European Research Council under the European Union’s Seventh Framework Programme (FP7/2007-2013) / ERC Grant agreement No. 291058. This research was supported by the Research Council of Norway through the grant “Solar Atmospheric Modelling”. This research was supported in part with computational resources at NTNU provided by NOTUR, http://www.sigma2.no. Some computations were performed on resources provided by the Swedish National Infrastructure for Computing (SNIC) at the High Performance Computing Center North at Umeå University and the PDC Centre for High Performance Computing (PDC-HPC) at the Royal Institute of Technology in Stockholm. CHIANTI is a collaborative project involving George Mason University, the University of Michigan (USA) and the University of Cambridge (UK).

References

All Tables

All Figures

Fig. 1 Left panel: vertical magnetic field strength at height z = 0 Mm for timestep t = 5000 s. The color bar range is [–1.5 kG, +1.5 kG]. A set of magnetic field lines connecting the main magnetic polarities is overplotted in yellow in the left panel. Field line apices and cork positions of corks with T = 105 K are indicated by blue and red crosses. These field lines are further discussed in Sect. 4.4. Right panel: temperature distribution of the simulation grid points as a function of height for timestep t = 5000 s of the simulation. The average temperature at a given height is indicated by the gray line.

In the text

Fig. 2 Left panel: cork trajectories of three representative examples described in Sect. 4.1. The temporal evolution for 800 s is shown. As time evolves, the colors change from blue to red. The respective projections in the xy-, xz- and yz-plane are indicated in different shades of gray for the different corks, which have been labeled (a), (b), and (c). Right panel: shown from top to bottom for the three examples on the left are the temporal evolution of cork height, cork temperature, vertical velocity, and heating and cooling terms per particle, i.e., the sum of Joule heating and viscous heating (black), adiabatic compression/expansion (red), Spitzer heat conduction (blue) and radiative losses (green). In the velocity panel, negative velocities are pointing upward, positive velocities are pointing downward. The first column shows the evolution of cork (a), which is heated from chromospheric to upper transition region temperatures starting at ~ t = 4600 s; the second column shows the evolution of cork (b), which quickly cools from transition region to chromospheric temperatures at ~t = 4900 s; and the third column shows the evolution of cork (c), which remains in the upper transition region throughout the entire simulation.

In the text

Fig. 3 Spatial distribution of corks, which have initially been placed in the y = 12 Mm plane of the simulation box, 15 min into the simulation. Every single dot indicates the position of one cork. Dark regions indicate regions, where the cork density is high; white regions represent voids of corks, i.e., regions where no corks are present. As time evolves, the corks start moving in both horizontal and vertical direction, thus slowly moving out of the y = 12 Mm plane (see Sect. 4.1). A movie showing the evolution over 25 min is available online.

In the text

Fig. 4 Probability density functions (PDFs) of cork temperature. The PDFs are shown with respect to t = 5000 s for t = 4940 s (left panel) and t = 5060 s (right panel; see Sect. 4.2). The distribution of corks at t = 5000 s is shown by a turquoise line in each panel. The average cork temperature at each timestep is indicated by the red line. The blue lines show the minimum and maximum temperature of all corks for the respective height.

In the text

	Fig. 5 One-minute evolution of temperature, density, and height distributions of corks in the temperature range log (T∕[K]) = [4.9, 5.1] at t = 5000 s. The evolution between t – 30 s and t + 30 s is shown by the different colors. The respective mean values are indicated by dashed lines.
In the text

Fig. 6 Upper left panel: temporal evolution of cork masses in a given temperature bin for a 10-min timespan around t = 5000 s. Only atmospheric corks are considered. The temperature bins are divided into temperature intervals of log T = 0.2 dex, and each bin is represented by a different color (see legend in panel labeled [4.9, 5.1]). The fraction with respect to the total mass of all atmospheric corks is shown (see Sect. 4.3). Rows 1 (middle and right panel)–4: temporal evolution (t  ±   300 s) of the mass of corks within a given temperature bin. Only corks that contribute to the selected bin are considered, and all panels are normalized to 100. The same legend applies for all panels.

In the text

Fig. 7 Left panel: mass fraction (solid lines) and surface area (dashed lines) of upflowing (blue) and downflowing (red) plasma in various temperature intervals in the range log T = [3.5, 6.1]. Temporal averages based on the timesteps between t = 4500 s and t = 5300 s are presented. Right panel: average vertical velocity (light gray line) and density-squared-weighted vertical velocity (black line) as a function of temperature. The latter represents the trend of the average Doppler shift with line formation temperature. The error bars show the mean error.

In the text

	Fig. 8 Left: xy-positions of corks with logT = [4.99, 5.01] at t = 5000 s, color-coded with vertical cork velocity. Red indicates downflows; blue shows upflows. Right: C iv (1548 Å) Doppler shift map for t = 5000 s. Blue indicates blueshifts (upflows), red shows redshifts (downflows). The color bar on the right shows the velocity scales; the same range has been chosen for both panels.
In the text

Fig. 9 Evolution of field line parameters extracted for a low-lying field line including a downward moving transition region cork. Top row: projections of the field line on the xz-, yz-, and xy-planes. The field line parameters are shown as a function of the coordinate s along the field line, where s = 0 Mm corresponds to the position of the right footpoint of the loop in the panels of the top row. The evolution of the field line parameters is shown for t = 5000  ± 30 s, and the field line evolves from dark blue colors to red colors. The cork positions are indicated by asterisks and the loop apices by crosses.

In the text

	Fig. 10 Correlation between maximum pressure gradient along the loop and u∥,z for a set of downward moving corks on low-lying field lines. The set of magnetic field lines is shown together with the vertical magnetic field strength at height z = 0 Mm in Fig. 1 (left panel).
In the text