Photospheric Motions and Their Effects on the Corona:
a Numerical Approach

The simulations are carried out with the 3-D MHD code extensively described in Gudiksen & Nordlund (2005b), which solves the fully compressible MHD equations:

Here $\rho$, u, $P$, J, B, g, $\mu$, E and $e$ denote the mass density, velocity vector, gas pressure, electric current density, magnetic flux density, gravitational acceleration, the vacuum permeability, electric field vector and the internal energy per unit volume respectively; The terms ${\underline{\underline{\tau}}}$ and $\eta$ refer to the viscous stress tensor and the magnetic diffusivity; $Q_{visc}$ and $Q_{Joule}$ represent the viscous and Joule heating.

The set of MHD equations is computed on a staggered mesh of $100^{3}$ grid points corresponding to a box of $20\times 20\times 20~\rm{Mm^{3}}$ (distance between two grid points of $0.2~\rm{Mm}$), employing 6th-order derivative operators, 5th-order interpolation operators to align physical variables in space, and a 3rd-order Runge-Kutta time-stepping scheme.

Simulations on this work use a classical heat conduction along the magnetic field, Spitzer conductivity (Spitzer, 1956) and also a cooling function representing radiative losses in the optically thin coronal plasma. Parker (1972) argued that the vigorous formation of tangential discontinuities (current sheets), magnetic dissipation is practically independent of the conductivity. Moreover, it was showed numerically by Galsgaard & Nordlund (1996) and also Hendrix et al. (1996) that the energy created through dissipation is not strongly dependent on small scales. One should expect that, by resolving the smallest scales in a simulation of DC heating, the dissipated energy could have a different value. Instead, high resolution does not give dramatic changes in results as the stress applied on the magnetic field is not increased but rather shifted to a fine scale. The important point on the ohmic/nanoflare heating hypothesis is the level of stress experienced by the magnetic field which does not depend on the resolution scale, so consistent results can be retrieved from simulations without the necessity of using very high resolution.

The MHD equations are solved in a 3D x-y-z coordinate system with periodic boundary conditions in the horizontal direction. The top of the box has zero-gradient of temperature at the boundary so no heat exchanges are allowed. Also the vertical gradient of the horizontal velocity is null at the upper boundary. Both the density and potential field are extrapolated in the ghost zones (regions outside the computational domain). Vertical velocity is set to zero on both top and bottom boundaries and the driver establishes the horizontal velocity at the bottom. The lower boundary is stressed by a time dependent velocity field, constructed from a Voronoi tesselation, reproducing the granulation pattern. The velocity field reproduces the geometrical pattern as well a the amplitude power spectrum of the velocity and the vorticity, leaving no free parameters.

The photospheric motions are produced through a Voronoi tessellation driver (Okabe et al., 2000) which reproduce the geometry and flow pattern of the granular surface of the Sun. A particular model of a Voronoi tesselation investigated by Schrijver et al. (1997), the so-called “multiplicatively weighted Voronoi” tessellation, was found to reproduce the granular pattern of the Sun very well. It segments the two-dimensional space into Voronoi regions $V(\textbf{p}_{i})$ defined by

with $w_{i}$ representing the different weights (greater than zero) associated with the generator points $\textbf{x}_{i}$, and $p$ denoting the index of the power law function which controls the tessellation. By using these weighted distances between the generated points, the generation of very large or very small granules can be avoided. By choosing the weights $w_{i}$ carefully, a uniform granular structure is produced with only acceptable sizes for the granules, ranging in size between 1 and 2 Mm, in the case of solar photospheric granules, to $\sim$ 30 Mm or larger for supergranules (Leighton, 1963).

The velocity profile used by the model is taken from Simon & Weiss (1989), who fitted observational data from the SOUP instrument on Spacelab 2, and then adjusted to better fit the solar velocity spectrum by

In the previous equation $v_{i}$ is a time-dependent quantity, proportional to $w_{i}(t)$ which is defined by

with $W_{i}$ being the maximum weight and $q$ an even integer related with the growth time of the granule. The velocity at the center of granules evolve in time along a Gaussian profile $v_{i}\propto e^{-(t/t_{0})^{2}}$.

To quantify the effect of the photospheric granular network on the corona, we set up three sets of simulations. In each set we vary one of three parameters describing the granules: The granular radius $r$, the granular lifetime $\tau$ and the velocity amplitude of the plasma $v_{0}$. Each set contains three simulations with one simulation being identical for all the sets so the total number of simulations is 7. The details of the simulations and sets are provided in Table 1.

Simulation I is common to all the three sets and corresponds to the closest reproduction of solar conditions, by choosing $r=2~\rm{Mm}$, $\tau=5~\rm{min}$ and $v_{0}=3~\rm{km~s^{-1}}$. It is important here to point out that these parameters are not entirely independent and setting these values does not produce granules that have exactly these characteristics. The actual values for the granular radius, granular lifetime and average velocity in the granular field will be discussed in Sec. 4.1.

In set 1 both the lifetime and velocity amplitude are kept constant, while the granular radius is varied from 2 to 5 Mm. In the second set we simulate granules with lifetimes of 5, 8 and 12 minutes, maintaining the size of the granules and velocity amplitude constant and equal to 2 Mm and $3~~\rm{km~s^{-1}}$ respectively. The third set includes simulations with three different flow velocities, 3, 6 and 10 $\rm{km~s^{-1}}$ and granule sizes and lifetimes set to 2 Mm and 5 minutes, respectively.

Our simulation domain is initially a current-free space containing vertical magnetic field lines extending from the base which coincides with the photosphere ($z=0$) to the top of the box, 20 Mm up in the corona ($z=H=20~\rm{Mm}$). We then let the code run for several hundreds of minutes. The temporal resolution, i.e. the interval between two consecutive snapshots produced by the model is roughly half a second, thus providing us with an almost continuous evolution of the desired variables. We should note that the time span of each simulation is dependent on how fast the electric current density reaches a steady-state. After this point (designated here as the turn-over point) we extend the simulation, in order to ensure that the steady-state was reached.

The details of the driving velocity field in the photosphere is important. As seen, all the three parameters controlling the velocity field have an impact on the $J_{A}^{2}$ and the turn over point.

In a simplified picture, the granular velocity field will advect magnetic footpoints radially in each granule, until they arrive at the intergranular boundary. Depending on the location of the surrounding granules, the magnetic footpoints will then either stay at the location in the intergranular lane or be advected almost perpendicular to their previous velocity until they reach an intersection between more than two granules, where the velocity is generally low. The development of new granules will make that process repeat itself. The direction of movement will be independent of the previous history of the magnetic footpoints, so they will experience a random walk, where the step length will be roughly the granular radius, unless the velocity of the flow inside the granules is so low that they will never have time to reach one granular radius before the granule disappears. In that case the random walk step length will be the product of the velocity inside the granule and the lifetime of the granule. A convincing argument that the magnetic footpoints do experience a random walk, is found by following magnetic footpoints in the granular velocity field, assuming they are passively advected (see Fig. 5 for an example of the path travelled by one footpoint). We have followed 1024 footpoints at the surface (forming a $32\times 32$ grid at the bottom of the box) over time for each simulation. We then plot the displacement from the original position ($R$) as a function of time. Figure 6 shows that total displacement of the footpoints in the granular field of simulation I, and a fit of the form

In the previous expression $C$ and $p$ are constants and a value of $p=0.5$ would correspond to theoretical random walk. The fitting curves have values of $p$ as indicated on Table 3 which are fairly close to the theoretical value. We will assume that the motion experienced by a magnetic footpoint is a random walk, and we can then fit the data with

The random walk is described by two parameters: the number of steps $N$ and the step length $\ell$. The two parameters cannot be found independently, so they are not easy to find as we have only snapshots in time of the simulation. They can be found indirectly if we can find the step length $\ell$ and the step time $\bar{\tau}$ (the time required for one step). If these can be found the number of steps can also be determined since $N=t/\bar{\tau}$ where $t$ is the elapsed time. Because we have decoupled the velocity inside the granules, their lifetimes and their size, we have to be careful with the estimation of the parameters. The step length is assumed to be the smallest of these two values: the product of velocity times lifetime or the average radius of the granules at the surface, i.e., the distance from the center of each granule to the intergranular lanes/ vertices. By making this assumption we assume that the footpoint starts its movement always at the center of a granular cell which is a strong simplification and it will be discussed later. A footpoint can initiate the motion anywhere inside the granule migrating then towards the intergranular lanes, probably lying there for a while until a new granule emerges and propels it again. So, the granular radius should be understood as an over-estimated step length and the maximum value that this parameter can take.

The average radius of the granules $r_{av}$ is calculated using the same approach taken in Wedemeyer-Böhm & Rouppe van der Voort (2009) for solar observations, in this case from a power spectrum of the velocity field on the $z=0$ plane (see Fig. 7 for an example of the spectrum as a function of wavenumber for one snapshot of simulation I), where the value on the x-axis corresponding to the peak of the curve provides the diameter of the granules. From Fig. 7 it is immediately seen that the power spectrum has a broad peak, around $D\simeq 2.8~\rm{Mm}$, which gives $r_{av}\simeq 1.4~\rm{Mm}$ for this simulation.

The average velocity $\bar{v}$ inside the granules is then calculated from the velocity profile (Eq. 8 and Fig. 8) using the average radius from the powerspectrum of the velocity field:

The average lifetime $\bar{\tau}$ of the granules is found by looking at the distribution of the number of granules according to their development. Figure 9 shows an example (simulation I) which shows that most of the granules are living within 3.5 minutes around the instant of maximum development ($t-t_{0i}$), which will correspond to an average lifetime of 7 minutes.

Initially $\ell$ is set equal to the product $\bar{v}\bar{\tau}$. In case the product $\bar{v}\bar{\tau}$ is greater than $r_{av}$, $r_{av}$ is taken as the step length.

Table 4 organizes all these quantities discussed above for our simulations.

Figure 10 (right) shows the correlation between the granular size and the time of the turn-over point of the simulations. A plausible explanation of this behavior consistent with Parker’s model and the random motion of the footpoints relies on the critical angle and on the fact that footpoints must be moved a certain distance from their original location in order to reach that angle. The properties of the random motion means that the total distance traversed by the footpoints is more sensitive to the step length than the number of steps. All simulations of set 1 (where the size of the granules is changed keeping the other initial parameters constant) have approximately the same step length but one has to be careful as the size of granules affects the results. We have been assuming the center of the granule as the original location of footpoints which is not necessarily true of course. Eventually a footpoint may start to be dragged very close to the intergranular lane so it cannot travel the distance equivalent to the step length. In this first set of simulations $\ell\simeq 1.10~\rm{Mm}$ which means that for a footpoint to be able to move this distance it has to be located at least one step length from the intergranular lane. If the footpoint starts to move from a point farther away from the intergranular lane, its motion will be limited to $\ell$ by the lifetime of the granule. There will be then an area of the granule around the center where footpoints can start the motion with high probability of traveling along $\ell$. This area inside granules is obviously larger and larger as we increase the size of the granules meaning that, statistically, the real step length will closely match the value found for larger granules and will progressively reduces as we decrease the granular size. Our values of step length correspond in this way to the maximum value the step length can take. Taking the size of the granules into account it is reasonable to think then that our value of the step length will diverge more and more from the real value as we decrease the size of the granules. This can justify the inverse correlation observed between the granular size and the time for the turn-over point: smaller granules will have a smaller and smaller value of step length and then the equilibrium between current generated and released in the corona will happen at a later time when compared with larger granules.

Figure 11 (right) shows the variation of the instant of turn-over point as a function of granular lifetime. Here the correlation is linearly increasing. According to the random walk theory there should be no dependence on the granular lifetime. Nevertheless we should not exclude the influence of the adjacent granules and their lifetime. For shorter lifetimes the topology of velocity field changes faster, and the magnetic field lines trapped in intergranular lanes will be moved sooner due to the faster changing velocity field topology. A second effect is that the electric current density is lower as we will see in Sec. 4.3.

The time of turn-over point is also influenced by the velocity field and the correlation is shown in Figure 12 (right). It is obvious that as we increase the parameter $v_{0}$ in the simulations, the average step velocity ($\bar{v}$) also increases approximately by the same amount, which will induce an increase of the step length. Associated with this fact, the electric current produced also dramatically increases and the time needed to balance this process by dissipating the larger amount of current will be longer than for low velocities.

At this point there are enough elements to give an estimate (in a rough way as we are doing some strong assumptions) for the critical angle introduced by Parker in his model. This critical angle should correspond to the average deviation from the vertical of the magnetic footpoints when the dissipation starts to cancel the transverse component of the magnetic field at the same rate it is created by the drag of footpoints in the photosphere. In this study that moment will coincide approximately with the time of turn-over point, when after the fast initial growth of electric current density we observe a stabilization of this quantity caused by the balance between creation and dissipation of electric current.

Using the values of the time of turn-over point, time step and step length and applying them on Eq. 15, one can calculate the average displacement of the magnetic footpoints from its original location

where $t_{turn}$ is the time of turn-over point and $\bar{\tau}$ the time step, so $t_{turn}/\bar{\tau}=N_{turn}$.

The estimated deviation angle from the vertical will then be

with $H$ being the height of the box, $H=20~\rm{Mm}$. Results are indicated on Table 5.

Due to all the assumptions taken along this work, it is hard to find correlations between the angle and the initial parameters, or at least to put a high level of trust in them. Nevertheless one can observe that there are no huge variations of the estimated angle which varies from 5.5 to approximately 12.5, thus leading to an average value around 10 degrees. This is not very far from the 14^∘ postulated by Parker and could mean that our assumptions are justifiable. There is a particular discrepancy in set 1 of experiments where the estimated angle decreases from 8.76^∘ in simulation I to 5.50^∘ in simulation III. According to 18, the critical angle is proportional to the displacement at the moment of turn-over point, which is proportional to the step length from 17. We should expect a smaller angle for smaller granular size as the turn-over point is reached faster leading to a shorter displacement and the calculated step length is about the same in all these 3 simulations. However it was already explained at the beginning of this section that the true step length should be smaller for smaller granules, so this fact would uniform the critical angle, making it approximately the same for the three simulations.

The energy dissipated by currents in the atmosphere originates as Poynting flux injected into the magnetic field by the photospheric motions. The amount of dissipated energy in the atmosphere must therefore be proportional with the injected energy flux when a statistical equilibrium is reached. If this is correct there should exist simple relations between the granular parameters and the average electric current in the atmosphere. Parker (1983a) estimated the Poynting flux injected into the atmosphere:

where $v$ is the velocity which the magnetic footpoint is moved with (measured in km$/$s) during the time interval $\Delta t$, $B$ is magnetic field strength, $B_{\perp}$ is the field strength of the transverse component of the magnetic field generated by the movement of the footpoint and $H$ is the height of box. As there is only a dependency on the flow velocity, the relations should here be rather simple. It is though important to notice that the velocity here is the effective velocity, meaning the velocity that carries a magnetic footpoint further away from it initial position in a straight line. It is therefore not that simple to estimate the effect of the variables defined earlier, $r_{av}$, $\bar{\tau}$, and $\bar{v}$.

Figure 10 (left) shows the dependence of the average current density in the atmosphere as a function of the average granular size and the good linear fit between these two quantities: the current density generated goes up as the average size of the granules increases. The step length calculated for set 1 is roughly the same but because of the reason already pointed out in the previous discussion about the variation of the instant of turn-over point with granular size, one should expect a smaller value of the step length than the calculated value for smaller granules. So in a photosphere with larger granules the motions of the plasma will cause a higher displacement in a straight line of the footpoints of the magnetic field lines, twisting the magnetic field more than in the case of smaller granules and thus producing more electric current. With the module of the Poynting vector proportional to the transverse component of the magnetic field, $B_{\perp}$, a larger granule, produces higher $B_{\perp}$ by creating a stronger twist of the magnetic field and therefore produces an increase of the energy flux.

Figure 11 (left) evinces a correlation between the average granular lifetime and $J_{A}^{2}$. Magnetic footpoints caught inside a granule that lives longer follow a straight path radially inside the granule for a longer time, and therefore, the effective velocity is higher. This happens as long as the granules are large enough for the magnetic footpoints not to reach the edges of the granules within the granular lifetime, but in this set 2 of simulations the granular size is kept constant and then the influence of the size of the granules is the same for the three experiments of this set. On the other hand, inside short-lived granules footpoints will move straight in shorter time steps decreasing the effective velocity and consequently the $J_{A}^{2}$.

In Figure 12 (left) we show an identical scenario for the variation of the $J_{A}^{2}$ as functions of the average velocity field at the photosphere. From equation 19 it seems that there should be a correlation between the total current squared and the velocity squared, but again the velocity in equation 19 is the straight line velocity while the velocity experienced by the magnetic footpoints are not the effective straight line velocity. The correlation between the $J_{A}^{2}$ and the velocity can then be explained by the random walk nature of the velocity field.

One can think, based on expression 12 that a similar correlation should be obtained between the average current density or the turn-over point and the product $\bar{v}\bar{\tau}$. Following the previous discussed relation of average current density and the angle $\theta$, the latter should increase with $\bar{v}\bar{\tau}$, and the same should apply for the turn-over point. Figure 13 shows this behavior, although it is not possible to linearly fit all the data points. This non-linear increase of both the average current density and turn-over point with $\bar{v}\bar{\tau}$ implies that these variables are not connected in a simple way, and that the relation between the amount of electric current arising in the corona and the deviation angle $\theta(t)$ to the vertical should be more complex. This is probably again an effect of the difference between the $\bar{v}$ and the velocity used in equation 19.