Interpreting magnetic helicity flux in solar flux emergence

Abstract

Magnetic helicity flux gives information about the topology of a magnetic field passing through a boundary. In solar physics applications, this boundary is the photosphere and magnetic helicity flux has become an important quantity in analysing magnetic fields emerging into the solar atmosphere. In this work we investigate the evolution of magnetic helicity flux in magnetohydrodynamic (MHD) simulations of solar flux emergence. We consider emerging magnetic fields with different topologies and investigate how the magnetic helicity flux patterns correspond to the dynamics of emergence. To investigate how the helicity input is connected to the emergence process, we consider two forms of the helicity flux. The first is the standard form giving topological information weighted by magnetic flux. The second form represents the net winding and can be interpreted as the standard helicity flux less the magnetic flux. Both quantities provide important and distinct information about the structure of the emerging field and these quantities differ significantly for mixed sign helicity fields. A novel aspect of this study is that we account for the varying morphology of the photosphere due to the motion of the dense plasma lifted into the chromosphere. Our results will prove useful for the interpretation of magnetic helicity flux maps in solar observations.

1 Introduction

The year 2018 marks the 60th anniversary of Woltjer’s observation (Woltjer 1958) that magnetic helicity in a tangent magnetic field is an invariant of ideal magnetohydrodynamics (MHD). That is, for a simply connected domain $\unicode[STIX]{x1D6FA}$ ,

(1.1) $$\begin{eqnarray}\frac{\text{d}}{\text{d}t}\int _{\unicode[STIX]{x1D6FA}}\boldsymbol{A}\boldsymbol{\cdot }\boldsymbol{B}\,\text{d}V=0,\quad \boldsymbol{B}\boldsymbol{\cdot }\boldsymbol{n}=0\quad \text{ on }\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA},\end{eqnarray}$$

where $\boldsymbol{B}$ is the magnetic field, $\boldsymbol{A}$ is the vector potential of the magnetic field and $\boldsymbol{n}$ is the surface normal. This form of helicity is also gauge invariant in the sense that the transformation $\boldsymbol{A}\rightarrow \boldsymbol{A}+\unicode[STIX]{x1D735}\unicode[STIX]{x1D712}$ , for some scalar function $\unicode[STIX]{x1D712}$ , does not change the value of helicity. A few years later, Moffatt (Moffatt 1969) discovered the topological connection of magnetic helicity to Gauss’ linking number. Thus, the subject of magnetic topology was born and has been influential in many areas of MHD ever since (e.g. Taylor 1974; Frisch et al. 1975; Berger 1984; Schindler, Hesse & Birn 1988; Pevtsov, Canfield & Metcalf 1995; DeVore 2000; Russell et al. 2015). Some important examples of applications include the prediction of relaxed magnetic states in reverse pinch experiments (Taylor 1974) and the analysis of emerging solar active regions (e.g. Pevtsov et al. 1995; Démoulin et al. 2002; Leka, Fan & Barnes 2005; LaBonte, Georgoulis & Rust 2007). The first example we have given is an instance of magnetic topology applied to magnetic relaxation. Taylor (1974) conjectured that (total) magnetic helicity is approximately conserved in a system with very small resistivity. Minimizing the magnetic energy whilst preserving the magnetic helicity, the relaxed state is predicted to be a linear force-free field. This simple and elegant theory has had good experimental corroboration and is one of the successes of the application of magnetic helicity.

The second application we gave, magnetic flux emergence, is the subject of this paper. ‘Bundles’ of magnetic field rise up through the solar convection zone until they reach the Sun’s surface, the photosphere. Here, the magnetic bundles emerge into the solar atmosphere where they can produce interesting phenomena such as flares and coronal mass ejections (e.g. Hood, Archontis & MacTaggart 2012; Cheung & Isobe 2014). Since we cannot observe magnetic fields in the solar atmosphere directly, we rely on indirect information about their structure. This is where magnetic helicity input can play an important role.

When observers study the magnetic fields of solar active regions, the photosphere acts as a lower boundary for the magnetic field. This fact introduces a complication in the original definition of magnetic helicity. It means that active region magnetic fields have a non-trivial normal boundary component at the photosphere and, hence, magnetic field lines which leave the domain. In such circumstances the helicity is no longer gauge invariant. We can bypass the above complication, however, by considering relative magnetic helicity (Berger & Field 1984) which provides a topological invariant for fields connected to a boundary. Magnetic helicity can, therefore, be used in the analysis of solar active regions, in particular as an indicator of eruptive activity (e.g. Yeates & Hornig 2016; Guo et al. 2017; Pariat et al. 2017).

Although some researchers study relative helicity in solar atmosphere (e.g. Valori, Démoulin & Pariat 2012; Yeates & Hornig 2016; Guo et al. 2017; Pariat et al. 2017), this quantity requires knowledge of the magnetic field’s full structure which is not observationally available. Observational studies instead measure the relative magnetic helicity flux through the photosphere (e.g. Chae 2001; Kusano et al. 2002; Pevtsov, Maleev & Longcope 2003; Yokoyama et al. 2003; Pariat et al. 2006; LaBonte et al. 2007; Schuck 2008; Vemareddy 2015). Photospheric maps of helicity flux during flux emergence can reveal interesting and complex patterns. Often, however, these results are interpreted based on a ‘standard model’ of flux emergence where the emerging magnetic field is a twisted flux tube just below the photosphere (e.g. Archontis & Török 2008; MacTaggart & Hood 2009; Moreno-Insertis & Galsgaard 2013). Twist is an important part of helicity that, in flux emergence studies, leads to eruptive behaviour in the atmosphere. Observations of active region helicity input, however, imply that even bipolar regions, associated with single sign helicity flux rope/sheared arcade formation, can often input both signs of helicity (e.g. Leka & Skumanich 1999; Pevtsov et al. 2008; Vemareddy 2015; Vemareddy & Démoulin 2017; Bi et al. 2018). Modelling and interpreting the emergence of mixed helicity regions should be a priority in order to provide an understanding of their behaviour.

The purpose of this article is to show how (relative) magnetic helicity flux can be interpreted in solar flux emergence with a variety of emerging magnetic topologies. We first examine the ‘topological meaning’ of helicity flux in terms of winding (field line entanglement). We then study the magnetic helicity and winding fluxes in flux emergence simulations with different initial profiles of the magnetic field. The paper ends with a summary and a discussion.

2 Magnetic helicity flux

Consider an emerging magnetic field subject to ideal motion and two points $\boldsymbol{a}_{1}(t)$ and $\boldsymbol{a}_{2}(t)$ which represent the intersection of two field lines and the photospheric plane $P$ . The motion of these points in time follows the motion of the points of intersection of these field lines in $P$ as the field emerges. The net rotation of these points within $P$ as a function of time is

(2.1) $$\begin{eqnarray}c(\boldsymbol{a}_{1},\boldsymbol{a}_{2},t)=\frac{1}{2\unicode[STIX]{x03C0}}\int _{t_{0}}^{t_{1}}\frac{\text{d}\unicode[STIX]{x1D703}_{12}(\boldsymbol{a}_{1},\boldsymbol{a}_{2},t)}{\text{d}t}\,\text{d}t,\end{eqnarray}$$

where $\unicode[STIX]{x1D703}_{12}(\boldsymbol{a}_{1},\boldsymbol{a}_{2},t)$ is the angle of the line segment linking $\boldsymbol{a}_{1}$ and $\boldsymbol{a}_{2}$ . Equation (2.1) represents the entanglement of a pair of field lines due to emergence or in plane motions of the field lines, hence $\text{d}c/\text{d}t$ represents the rate of input of entanglement through the photosphere and into solar atmosphere due to these field lines.

In order to determine the time derivative of the angle $\unicode[STIX]{x1D703}_{12}$ , first note that

(2.2) $$\begin{eqnarray}\unicode[STIX]{x1D703}_{12}=\text{arctan}\displaystyle \left(\frac{(\boldsymbol{a}_{2}-\boldsymbol{a}_{1})\boldsymbol{\cdot }\hat{\boldsymbol{y}}}{(\boldsymbol{a}_{2}-\boldsymbol{a}_{1})\boldsymbol{\cdot }\hat{\boldsymbol{x}}}\right).\end{eqnarray}$$

Differentiating gives,

(2.3) $$\begin{eqnarray}\frac{\text{d}\unicode[STIX]{x1D703}_{12}}{\text{d}t}=\frac{1}{(\boldsymbol{a}_{2}-\boldsymbol{a}_{1})\boldsymbol{\cdot }(\boldsymbol{a}_{2}-\boldsymbol{a}_{1})}\hat{\boldsymbol{z}}\boldsymbol{\cdot }(\boldsymbol{a}_{2}-\boldsymbol{a}_{1})\times \left(\frac{\text{d}\boldsymbol{a}_{2}}{\text{d}t}-\frac{\text{d}\boldsymbol{a}_{1}}{\text{d}t}\right).\end{eqnarray}$$

The (relative) magnetic helicity $H$ represents the total winding $c(\boldsymbol{a}_{1},\boldsymbol{a}_{2},t)$ over all pairs of paths weighted by the photospheric flux over a given period of time (Berger 1988; Démoulin & Berger 2003). Thus, the rate of input of magnetic helicity into the solar atmosphere, through the photosphere $P$ , is

(2.4) $$\begin{eqnarray}\frac{\text{d}H}{\text{d}t}=-\frac{1}{2\unicode[STIX]{x03C0}}\int _{P}\int _{P}B_{z}(\boldsymbol{a}_{1})B_{z}(\boldsymbol{a}_{2})\frac{\text{d}\unicode[STIX]{x1D703}_{12}(\boldsymbol{a}_{1},\boldsymbol{a}_{2},t)}{\text{d}t}\,\text{d}x_{1}\,\text{d}y_{1}\,\text{d}x_{2}\,\text{d}y_{2}.\end{eqnarray}$$

In order to use (2.4), we need to know how to track points $\boldsymbol{a}_{1}$ and $\boldsymbol{a}_{2}$ . This motion is given by both in-plane motions and the projection due to motion out of the plane (flux emergence or submergence). If $\boldsymbol{u}$ is the velocity field and $\boldsymbol{u}_{\Vert }$ is the projection of $\boldsymbol{u}$ onto $P$ , it can be shown that

(2.5) $$\begin{eqnarray}\frac{\text{d}\boldsymbol{a}}{\text{d}t}=\boldsymbol{u}_{\Vert }-\frac{u_{z}}{B_{z}}\boldsymbol{B}_{\Vert }=\boldsymbol{u}_{\Vert }-u_{z}\left(\frac{\text{d}\boldsymbol{a}}{\text{d}z}\right)_{\Vert },\end{eqnarray}$$

(e.g. Berger 1988; Démoulin & Berger 2003). The term $\boldsymbol{u}_{\Vert }$ accounts for the motions of the field lines within the plane. The term $u_{z}(\text{d}\boldsymbol{a}/\text{d}z)_{\Vert }$ accounts for the apparent motion within the plane due to the emergence (or submergence if $u_{z}<0$ ) of field lines. The total derivative $\text{d}\boldsymbol{a}/\text{d}t$ assumes that the field lines are advected with the plasma (ideal motion).

The helicity flux through $P$ represents the transmission of topological information into the atmosphere and can change in time. In magnetic relaxation, although magnetic helicity controls the final equilibrium, the path to that equilibrium can involve very complicated dynamics. Similarly in flux emergence, although magnetic helicity affects the overall evolution of the emerging field, the magnetic helicity flux is associated closely with the dynamics of emergence – the interplay of both plasma and magnetic field. For this reason, it is important to interpret carefully what maps of magnetic helicity flux mean.

A related quantity is the field line helicity input (e.g. Berger 1988; Vemareddy 2015; Vemareddy & Démoulin 2017) which we define here as

(2.6) $$\begin{eqnarray}\frac{\text{d}{\mathcal{H}}}{\text{d}t}(\boldsymbol{a}_{0})=-\frac{1}{2\unicode[STIX]{x03C0}}B_{z}(\boldsymbol{a}_{0})\int _{P}B_{z}(\boldsymbol{a})\frac{\text{d}\unicode[STIX]{x1D703}(\boldsymbol{a}_{0},\boldsymbol{a})}{\text{d}t}\,\text{d}x\,\text{d}y.\end{eqnarray}$$

Equation (2.6) describes the contribution to the helicity flux $\text{d}H/\text{d}t$ from a single point $\boldsymbol{a}_{0}\in P$ (thus the contribution to helicity flux from a single field line). Although the integrand of (2.6) has a singularity at $\boldsymbol{a}_{0}$ , the integral still converges as the pole is of a lower order than the domain of integration. In the solar physics literature, $\text{d}{\mathcal{H}}(\boldsymbol{a}_{0})/\text{d}t$ is often labelled ‘ $G_{\unicode[STIX]{x1D703}}$ ’ and is used as a proxy for helicity flux density (e.g. Pariat, Démoulin & Berger 2005).

The contribution from a neighbourhood around the pole (call it $a_{0}^{\unicode[STIX]{x1D716}}$ ) can be shown to be decomposed into the local twisting of the field around this point (which relates to the local electric current) and a contribution called the writhe, due to the geometry of the field line passing through $\boldsymbol{a}_{0}$ (Călugăreanu 1959, 1961; White 1969; Berger & Prior 2006). Considering the contribution of the integral (2.6) from $a_{0}^{\unicode[STIX]{x1D716}}$ , we can relate this to self-helicity as it measures the winding of the field in the locality of the field line passing through $\boldsymbol{a}_{0}$ . The contribution due to the rest of the field, on the set $P-a_{0}^{\unicode[STIX]{x1D716}}$ , with this local neighbourhood is related to mutual helicity since the contribution represents the winding of the rest of the field with the small ‘flux tube’ in $a_{0}^{\unicode[STIX]{x1D716}}$ . This idea of decomposing the helicity into mutual and self-components has been used in numerous studies. Often it is the case that the field will have finite volume flux ropes, as in Pariat et al. (2006) and Guo et al. (2017), be the result of clear large scale shearing motion, as in Démoulin et al. (2002), or have regions of similar squashing factor, such as in Guo et al. (2017). In these cases it is possible to extend the notion of self-helicity to these finite bundles of field lines and gain further insight into the field’s morphology by tracking these self/mutual interactions. In this study, we will consider fields that exhibit highly complex local current structure and topology. Hence, a decomposition of the helicity into self and mutual components is not a simple task or necessarily well defined. Therefore, in this study, we do not pursue the self/mutual approach to describing helicity.

As a final note, before presenting our analysis, we add that Pariat et al. (2006) showed that it is possible to get spurious input values from (2.6) due to field structures with no helicity. This is not the case for the emergence/submergence events analysed here but is a possibility that should always be considered when looking at the distributions of this quantity (such spurious contributions vanish for the net input (2.4)).

2.1 Separating field strength and topological information

To disentangle the field strength (flux) and topological information associated with the helicity input, we can also consider the net winding input

(2.7) $$\begin{eqnarray}\frac{\text{d}L}{\text{d}t}=-\frac{1}{2\unicode[STIX]{x03C0}}\int _{P}\int _{P}\unicode[STIX]{x1D70E}_{z}(\boldsymbol{a}_{1})\unicode[STIX]{x1D70E}_{z}(\boldsymbol{a}_{2})\frac{\text{d}\unicode[STIX]{x1D703}_{12}(\boldsymbol{a}_{1},\boldsymbol{a}_{2},t)}{\text{d}t}\,\text{d}x_{1}\,\text{d}y_{1}\,\text{d}x_{2}\,\text{d}y_{2},\end{eqnarray}$$

where

(2.8) $$\begin{eqnarray}\unicode[STIX]{x1D70E}_{z}(\boldsymbol{x})=\left\{\begin{array}{@{}cc@{}}1 & \text{ if }B_{z}(\boldsymbol{x})>0,\\ -1 & \text{ if }B_{z}(\boldsymbol{x})<0,\\ 0 & \text{ if }B_{z}(\boldsymbol{x})=0,\end{array}\right.\end{eqnarray}$$

(e.g. Prior & Yeates 2014). The rate $\text{d}L/\text{d}t$ measures the input of field line entanglement through the photosphere and ignores the extra flux information provided by the helicity. A similar suggestion made in the literature is to divide the helicity input by the square of the flux (e.g. Yamamoto et al. 2005; LaBonte et al. 2007; Vemareddy & Démoulin 2017). We propose that $\text{d}L/\text{d}t$ is more meaningful in this context as it is based on a quantity, the net winding, which is also an ideal invariant (Prior & Yeates 2014). Indeed, it was shown recently to be more topologically meaningful than helicity as it can be used, under certain circumstances, to classify field topologies and hence precisely measure changing field connectivity (Prior & Yeates 2018).

2.2 Evolving photosphere

The standard practice in both theoretical and observational studies of helicity is to treat the lower (photospheric) boundary as fixed. In almost all studies, the photosphere is a plane (e.g. Pariat et al. 2005, 2006; Jeong & Chae 2007; Sturrock et al. 2015; Vemareddy 2015; Vemareddy & Démoulin 2017), although calculations in spherical geometry have also been performed (Berger 1985; MacTaggart et al. 2016; Moratis et al. 2018). In our calculations, which we will describe shortly, our fixed photospheric plane $P$ will represent the fixed height $z=0$ in a Cartesian domain.

In addition to this approach, we also consider the effect of an evolving photosphere. Including this effect means that the lower boundary is no longer fixed but is a surface that deforms in space and time. In our calculations, we will define the photospheric boundary as the surface where the plasma density $\unicode[STIX]{x1D70C}$ has the non-dimensional value of 1. In order to evaluate how a changing photosphere will affect the helicity and winding fluxes, we use the following procedure. First, we calculate the varying surface on which $\unicode[STIX]{x1D70C}=1$ , i.e. the set

(2.9) $$\begin{eqnarray}\left\{P_{v}(\boldsymbol{a},t)\,|\,\boldsymbol{a}\in P,\,\unicode[STIX]{x1D70C}(P_{v})=1\right\}.\end{eqnarray}$$

$P$ will always represent the horizontal plane at $z=0$ in our model (more details will be given later). At $t=0$ in the simulations, $P_{v}\equiv P$ . As indicated, $P_{v}$ is a function of the coordinates $(x,y)$ of $P$ .

We then evaluate $\text{d}\boldsymbol{a}/\text{d}t$ at $P_{v}(\boldsymbol{a},t)$ rather than at $P(\boldsymbol{a})$ . We can use this quantity to calculate an effective rotational derivative $\text{d}\unicode[STIX]{x1D703}_{12}/\text{d}t$ using (2.3). To account for the surface geometry, we replace the area element $\text{d}x_{1}\text{d}y_{1}$ with $J(\boldsymbol{a}_{1})\text{d}x_{1}\text{d}y_{1}$ , where $J$ the Jacobian of the surface $P_{v}$ at the coordinates of $\boldsymbol{a}_{1}$ (and similarly for $\boldsymbol{a}_{2}$ ). For example, the calculation (2.7) will become

(2.10) $$\begin{eqnarray}\displaystyle \frac{\text{d}L}{\text{d}t} & = & \displaystyle -\frac{1}{2\unicode[STIX]{x03C0}}\int _{P}\int _{P}\bigg[\unicode[STIX]{x1D70E}_{z}(P_{v}(\boldsymbol{a}_{1},t))\unicode[STIX]{x1D70E}_{z}(P_{v}(\boldsymbol{a}_{2},t))\frac{\text{d}\unicode[STIX]{x1D703}_{12}(P_{v}(\boldsymbol{a}_{1},t),P_{v}(\boldsymbol{a}_{2},t))}{\text{d}t}\nonumber\\ \displaystyle & & \displaystyle \,\times J(P_{v}(\boldsymbol{a}_{1},t))J(P_{v}(\boldsymbol{a}_{2},t))\bigg]\text{d}x_{1}\,\text{d}y_{1}\,\text{d}x_{2}\,\text{d}y_{2}.\end{eqnarray}$$

Equation (2.10) does not represent the exact winding flux through a spatially non-uniform surface. We could, for example, have made use of the differential form definition of relative helicity given in Chapter 3 of Arnold & Khesin (1999) in order to calculate the actual flux. However, we argue that what the above procedure represents is more akin to the effective projection of information onto the plane $P$ which occurs in the creation of magnetogram and vector magnetogram data (e.g. Démoulin & Pariat 2009; Scherrer et al. 2012). These magnetogram data are obtained from line-of-sight information, Zeeman splitting and normal to line-of-sight linear polarization information. These data are effectively projected onto a plane $P$ . Since we are trying to gauge the consequences of an effect which is already uncertain (the interpretation of optical information) we feel the above approach is sensible first step. To the best of our knowledge, this is the first theoretical helicity study which has attempted to take account of a moving photosphere.

3 Simulation set-up and initial conditions

In this study, we will consider small active regions of photospheric area ${\sim}75~\text{Mm}^{2}$ . Studying regions of this size is common in the flux emergence literature (Hood et al. 2012; Cheung & Isobe 2014) and is justified on the grounds of achieving suitable spatial and temporal resolution. Of course, care must be taken when comparing the results of such studies to observations of regions of different sizes. Although qualitative behaviour may be found across spatial and temporal scales, quantitative information will be different.

As our focus is on the dynamics of the magnetic field, we will consider an idealized description of the solar atmosphere. The bulk properties of the plasma and magnetic field dynamics are described by compressible MHD. The three-dimensional (3-D) resistive and compressible MHD equations are solved using a Lagrangian remap scheme (Arber et al. 2001). In dimensionless form, the MHD equations are

(3.1) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{D}\unicode[STIX]{x1D70C}}{\text{D}t}=-\unicode[STIX]{x1D70C}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}, & \displaystyle\end{eqnarray}$$

(3.2) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{D}\boldsymbol{u}}{\text{D}t}=-\frac{1}{\unicode[STIX]{x1D70C}}\unicode[STIX]{x1D735}p+\frac{1}{\unicode[STIX]{x1D70C}}(\unicode[STIX]{x1D735}\times \boldsymbol{B})\times \boldsymbol{B}+\frac{1}{\unicode[STIX]{x1D70C}}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D748}+\boldsymbol{g}, & \displaystyle\end{eqnarray}$$

(3.3) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{D}\boldsymbol{B}}{\text{D}t}=(\boldsymbol{B}\boldsymbol{\cdot }\unicode[STIX]{x1D735})\boldsymbol{u}-(\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u})\boldsymbol{B}+\unicode[STIX]{x1D702}\unicode[STIX]{x1D6FB}^{2}\boldsymbol{B}, & \displaystyle\end{eqnarray}$$

(3.4) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{D}\unicode[STIX]{x1D700}}{\text{D}t}=-\frac{p}{\unicode[STIX]{x1D70C}}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}+\frac{1}{\unicode[STIX]{x1D70C}}\unicode[STIX]{x1D702}|\boldsymbol{j}|^{2}+\frac{1}{\unicode[STIX]{x1D70C}}Q_{\text{visc}}, & \displaystyle\end{eqnarray}$$

(3.5) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{B}=0, & \displaystyle\end{eqnarray}$$

with specific energy density

(3.6) $$\begin{eqnarray}\unicode[STIX]{x1D700}=\frac{p}{(\unicode[STIX]{x1D6FE}-1)\unicode[STIX]{x1D70C}}.\end{eqnarray}$$

The basic variables are the density $\unicode[STIX]{x1D70C}$ , the pressure $p$ , the magnetic induction $\boldsymbol{B}$ (referred to as the magnetic field) and the velocity $\boldsymbol{u}$ . $\boldsymbol{j}$ is the electric current density, $\boldsymbol{g}$ is the surface gravitational acceleration (uniform in the $z$ -direction) and $\unicode[STIX]{x1D6FE}=5/3$ is the ratio of specific heats. The dimensionless temperature $T$ can be found from

(3.7) $$\begin{eqnarray}T=(\unicode[STIX]{x1D6FE}-1)\unicode[STIX]{x1D700}.\end{eqnarray}$$

We make the variables dimensionless against photospheric values, a standard procedure in flux emergence studies (e.g. Hood et al. 2012, and references therein). We have pressure $p_{\mathit{ph}}=1.4\times 10^{4}$ Pa; density $\unicode[STIX]{x1D70C}_{\mathit{ph}}=2\times 10^{-4}$ kg m $^{-3}$ ; scale height $H_{\mathit{ph}}=170$ km; surface gravitational acceleration $g_{\mathit{ph}}=2.7\times 10^{2}$ m s $^{-2}$ ; speed $u_{\text{ph}}=6.8$ km s $^{-1}$ ; time $t_{\mathit{ph}}=25$ s; magnetic field strength $B_{\mathit{ph}}=1.3\times 10^{3}$ G and temperature $T_{\mathit{ph}}=5.6\times 10^{3}$ K. In the non-dimensionalization of the temperature we use a gas constant ${\mathcal{R}}=8.3\times 10^{3}$ m $^{2}$  s $^{-2}$  K $^{-1}$ and a mean molecular weight $\tilde{\unicode[STIX]{x1D707}}=1$ . $\unicode[STIX]{x1D702}$ is the resistivity and we take its value to be $10^{-3}$ . This value is close to the lowest physical resistivity that can be chosen before numerical resistivity dominates (see Arber, Haynes & Leake 2007; Leake, Linton & Török 2013). Note that in the simulations we present, reconnection plays a minor role. Thus, to a good approximation, the motion in the simulations can be considered to be ideal. The fluid viscosity tensor and the viscous contribution to the energy equation are respectively

(3.8) $$\begin{eqnarray}\unicode[STIX]{x1D748}=2\unicode[STIX]{x1D707}\left[\boldsymbol{D}-{\textstyle \frac{1}{3}}(\text{tr}\boldsymbol{D})\unicode[STIX]{x1D644}\right]\end{eqnarray}$$

and

(3.9) $$\begin{eqnarray}Q_{\text{visc}}=\unicode[STIX]{x1D748}:\unicode[STIX]{x1D735}\boldsymbol{u},\end{eqnarray}$$

where

(3.10) $$\begin{eqnarray}\boldsymbol{D}={\textstyle \frac{1}{2}}\left(\unicode[STIX]{x1D735}\boldsymbol{u}+\unicode[STIX]{x1D735}\boldsymbol{u}^{\text{T}}\right)\end{eqnarray}$$

is the symmetric part of the rate of strain tensor and $\unicode[STIX]{x1D644}$ is the identity tensor. We take $\unicode[STIX]{x1D707}=10^{-5}$ and use this form of viscosity primarily to aid stability. The code accurately resolves shocks by using a combination of shock viscosity (Wilkins 1980) and Van Leer flux limiters (Van Leer 1979), which add heating terms to the energy equation. Values will be expressed in non-dimensional form unless explicitly stated otherwise.

The equations are solved in a Cartesian computational box of (non-dimensional) sizes $[-45,45]\times [-45,45]\times [-30,65]$ in the $x$ , $y$ and $z$ directions respectively. The boundary conditions are closed on the top and base of the box and periodic on the sides. Damping layers are included at the side and top boundaries to reduce the reflection/transmission of waves. The computational mesh contains 486 $\times$ 486 $\times$ 729 points.

3.1 Initial background atmosphere and velocity perturbation

The idealized initial equilibrium atmosphere is given by prescribing the temperature profile

(3.11) $$\begin{eqnarray}T(z)=\left\{\begin{array}{@{}cc@{}}1-{\displaystyle \frac{\unicode[STIX]{x1D6FE}-1}{\unicode[STIX]{x1D6FE}}}z, & z<0,\\ 1, & 0\leqslant z\leqslant 10,\\ 150^{[(z-10)/10]}, & 10<z<20,\\ 150, & z\geqslant 20.\end{array}\right.\end{eqnarray}$$

Starting from the top of (3.11), the sections represent the solar interior, the photosphere/chromosphere, the transition region and the corona. The above model temperature profile is standard in many works on flux emergence (e.g. Murray et al. 2006; Fan 2009; Sturrock et al. 2015). The above temperature profile is also used in Prior & MacTaggart (2016).

The other state variables, pressure and density, are found by solving the hydrostatic equation in conjunction with the ideal equation of state

(3.12) $$\begin{eqnarray}\frac{\text{d}p}{\text{d}z}=-\unicode[STIX]{x1D70C}g,\quad p=\unicode[STIX]{x1D70C}T.\end{eqnarray}$$

To study emergence, we must place a particular form for the magnetic field in the solar interior and apply a perturbation to allow it to emerge. In the following simulations, we will apply an initial velocity perturbation of the form

(3.13) $$\begin{eqnarray}\boldsymbol{u}\boldsymbol{\cdot }\hat{\boldsymbol{z}}=u_{0}\exp \left(-\frac{x^{2}}{x_{0}^{2}}\right)\exp \left(-\frac{y^{2}}{y_{0}^{2}}\right)\exp \left(-\frac{(z+\unicode[STIX]{x1D716}_{0}-R)^{2}}{z_{0}^{2}}\right)\sin \left(\frac{t}{t_{0}}\unicode[STIX]{x03C0}\right),\end{eqnarray}$$

where we set the constants $u_{0}=0.05$ , $x_{0}=5$ , $y_{0}=3$ , $z_{0}=5$ , $\unicode[STIX]{x1D716}_{0}=2.5$ and $t_{0}=6$ ; $R$ is the major radius of the tube and is described below. After $t=6$ the perturbation is switched off.

We will consider two magnetic field models. The first model is a twisted toroidal tube that has been used in many other studies (e.g. MacTaggart & Hood 2009; MacTaggart & Haynes 2014; MacTaggart et al. 2015; Sturrock et al. 2015). The second model is a ‘mixed helicity’ field which has a complex topology but zero total helicity and will serve as a model for mixed helicity emergence. We now discuss the construction of such fields.

3.2 How to construct mixed helicity fields

The following represents a specific case of the general mathematical form of magnetic flux ropes with general field line topology, introduced in Prior & Yeates (2016) and first applied to flux emergence in Prior & MacTaggart (2016). We assume the flux rope has a toroidal shape with an axis curve $\boldsymbol{r}(s)$ parameterized by its arclength $s$ as

(3.14) $$\begin{eqnarray}\boldsymbol{r}(s)=(-R\cos (s/R),0,R\sin (s/R)+z_{0}),\quad s\in [0,\unicode[STIX]{x03C0}R],\end{eqnarray}$$

where $R$ is the major radius of the torus and $z_{0}$ is the height of the rope’s anchoring footpoints at the base of the computational domain. This expression can be used to define a moving orthonormal frame $\{\boldsymbol{T},\boldsymbol{d}_{1},\boldsymbol{d}_{2}\}$ for $\boldsymbol{r}(s)$ , which can, in turn, be used to define a tubular coordinate system through the map $\boldsymbol{f}(s,x_{1},x_{2})$ , where

(3.15) $$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{f}(s,x_{1},x_{2})=\boldsymbol{r}(s)+x_{1}\boldsymbol{d}_{1}+x_{2}\boldsymbol{d}_{2}, & \displaystyle\end{eqnarray}$$

(3.16) $$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{d}_{1}=\left(\cos (s/R),0,-\sin (s/R)\right),\quad \boldsymbol{d}_{2}=(0,1,0). & \displaystyle\end{eqnarray}$$

From this mapping, we can find a set of basis vector fields $\{\frac{\unicode[STIX]{x2202}\boldsymbol{f}}{\unicode[STIX]{x2202}s},\frac{\unicode[STIX]{x2202}\boldsymbol{f}}{\unicode[STIX]{x2202}x_{1}},\frac{\unicode[STIX]{x2202}\boldsymbol{f}}{\unicode[STIX]{x2202}x_{2}}\}$ which can be used to define general vector fields in this toroidal coordinate system as

(3.17) $$\begin{eqnarray}\boldsymbol{B}=B_{s}\frac{\unicode[STIX]{x2202}\boldsymbol{f}}{\unicode[STIX]{x2202}s}+B_{1}\frac{\unicode[STIX]{x2202}\boldsymbol{f}}{\unicode[STIX]{x2202}x_{1}}+B_{2}\frac{\unicode[STIX]{x2202}\boldsymbol{f}}{\unicode[STIX]{x2202}x_{2}}.\end{eqnarray}$$

For such fields to be divergence free the following condition must hold

(3.18) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}JB_{s}}{\unicode[STIX]{x2202}s}+\frac{\unicode[STIX]{x2202}JB_{1}}{\unicode[STIX]{x2202}x_{1}}+\frac{\unicode[STIX]{x2202}JB_{2}}{\unicode[STIX]{x2202}x_{2}}=0,\end{eqnarray}$$

where $J=(R-x_{1})/R$ is the Jacobian of the map $\boldsymbol{f}$ . It follows that a magnetic field defined in a cylinder can be transferred to a torus simply by dividing the components by $J$ . In particular, we take the mixed helicity (braided) field used in numerous studies (e.g. Wilmot-Smith, Pontin & Hornig 2010; Wilmot-Smith et al. 2011; Pontin & Hornig 2015), which is composed of a series of overlapping counter-twists. The magnetic field has no helicity (the average total entanglement) but subsets of its field lines are braided. In our toroidal coordinate system this field is composed of exponential twists $\boldsymbol{B}_{t}(b_{0},k,a,l,x_{1c},x_{2c},s_{c})$ given by

(3.19) $$\begin{eqnarray}\displaystyle \hspace{-12.0pt} & \displaystyle \boldsymbol{B}_{t}(b_{0},k,a,l,x_{1c},x_{2c},s_{c})=\frac{2b_{0}k}{aJ}\text{exp}\left(-\frac{(x_{1}-x_{1c})^{2}+(x_{2}-x_{2c})^{2}}{a^{2}}-\frac{(s-s_{c})^{2}}{l^{2}}\right)\boldsymbol{R},\qquad & \displaystyle\end{eqnarray}$$

(3.20) $$\begin{eqnarray}\displaystyle \hspace{-10.79993pt} & \displaystyle \boldsymbol{R}=-(x_{2}-x_{2c})\frac{\unicode[STIX]{x2202}\boldsymbol{f}}{\unicode[STIX]{x2202}x_{1}}+(x_{1}-x_{1c})\frac{\unicode[STIX]{x2202}\boldsymbol{f}}{\unicode[STIX]{x2202}x_{2}}, & \displaystyle\end{eqnarray}$$

where the parameter $b_{0}$ determines the strength of the field, $a$ the horizontal width of the twist zones, $l$ their vertical extent and $k$ the handedness of the twist ( $k=1$ is right handed). The centre of rotation is $(x_{1c},x_{2c},s_{c})$ . The mixed helicity field is then defined as a superposition of $n$ pairs of positive and negative twists and an axial background field

(3.21) $$\begin{eqnarray}\displaystyle \boldsymbol{B}_{br}(b_{0},a,l,d,R,n) & = & \displaystyle \mathop{\sum }_{i=1}^{n}\left[\boldsymbol{B}_{t}(b_{0},1,a,l,0,-d,s_{d}i)\right.\nonumber\\ \displaystyle & & \displaystyle +\,\left.\boldsymbol{B}_{t}(b_{0},-1,a,l,0,d,s_{d}(i+1))\right]+\frac{b_{0}}{J}\frac{\unicode[STIX]{x2202}\boldsymbol{f}}{\unicode[STIX]{x2202}s},\end{eqnarray}$$

where $s_{d}=\unicode[STIX]{x03C0}R/(2n+1)$ and $d$ is the axial offset of the two opposing twists (see figure 1). The final component is the axial field (see figure 1 b). Finally, to express (3.21) in the ambient Cartesian coordinate system $(x,y,z)$ , we use the following maps between the tubular coordinates $(s,x_{1},x_{2})$ and the Cartesian coordinates

(3.22) $$\begin{eqnarray}s=\arctan (z,x)+\frac{\unicode[STIX]{x03C0}}{2},\quad x_{1}=(x+R\cos s)\cos s-(z-R\sin s)\sin s,\quad x_{2}=y,\end{eqnarray}$$

where the branch cut for the arctan function is at $\unicode[STIX]{x03C0}$ . For a magnetic field confined to a tube of minor radius $r$ , the field $\boldsymbol{B}_{b}(x,y,z)$ takes the following form

(3.23) $$\begin{eqnarray}\boldsymbol{B}_{b}(x,y,z)=\left\{\begin{array}{@{}cc@{}}\boldsymbol{B}_{br}(b_{0},a,l,d,R,n) & \text{ if }x_{1}^{2}(x,y,z)+x_{2}^{2}(x,y,z)\leqslant r,\\ 0 & \text{ if }x_{1}^{2}(x,y,z)+x_{2}^{2}(x,y,z)>r.\end{array}\right.\end{eqnarray}$$

For completeness, we can define a uniformly twisted field $\boldsymbol{B}_{t}$ as

(3.24) $$\begin{eqnarray}\boldsymbol{B}_{tw}(x,y,z)=\left\{\begin{array}{@{}cc@{}}{\displaystyle \frac{\unicode[STIX]{x1D719}}{J}}\boldsymbol{R} & \text{ if }x_{1}^{2}(x,y,z)+x_{2}^{2}(x,y,z)\leqslant r,\\ 0 & \text{ if }x_{1}^{2}(x,y,z)+x_{2}^{2}(x,y,z)>r.\end{array}\right.\end{eqnarray}$$

with $\unicode[STIX]{x1D719}$ determining the rate of rotation of the field and $\boldsymbol{R}$ given by (3.20). The actual twisted model that we will use in the simulations is that in (3.24) but weighted with an exponential term (see MacTaggart & Hood 2009). In these studies, we set the tube parameters to be $z_{0}=-30$ , $R=17.5$ and $r=2.5$ , so that the flux rope is anchored at the bottom boundary and its maximum initial height is $-10$ . The mixed helicity parameters are $a=\sqrt{0.04}$ , $l=0.04\unicode[STIX]{x03C0}R$ , $d=2.5/3$ , $b_{0}=5$ and we consider $n=2$ .

Figure 1. An example of the mixed helicity field we consider. (a) A field with $n=2$ , i.e. two right-handed twist elements (red) and two left handed (blue). (b) Shows the background axial field which is added to (a) in equation (3.21).

Finally, we remark that the more general field specification in Prior & Yeates (2016) allows for arbitrary tube shapes $\boldsymbol{r}(s)$ (as well as varying tube radius). The frame vectors $\boldsymbol{d}_{1}(s)$ and $\boldsymbol{d}_{2}(s)$ are defined by parallel transport, which requires an arclength integration. Thus, in general, the relationship between the tube coordinates $(s,x_{1},x_{2})$ and $(x,y,z)$ will require numerical integration.

4 Quantities analysed

We now list the quantities we will make use of in our analysis of the flux emergence simulations.

4.1 Field line helicity input rate

For completeness, we restate (2.6),

(4.1) $$\begin{eqnarray}\frac{\text{d}{\mathcal{H}}}{\text{d}t}(\boldsymbol{a}_{0})=-\frac{1}{2\unicode[STIX]{x03C0}}B_{z}(\boldsymbol{a}_{0})\int _{P}B_{z}(\boldsymbol{a})\frac{\text{d}\unicode[STIX]{x1D703}(\boldsymbol{a}_{0},\boldsymbol{a})}{\text{d}t}\,\text{d}x\,\text{d}y,\end{eqnarray}$$

which represents the contribution to the magnetic helicity input rate $\text{d}H/\text{d}t$ from a single point $\boldsymbol{a}_{0}\in P$ . For the moving photosphere calculation we evaluate the field at the points of the surface $P_{v}(\boldsymbol{a},t)$ and account for the surface Jacobian, i.e. $\text{d}x\,\text{d}y\rightarrow J\text{d}x\,\text{d}y$ . We label this quantity $\text{d}{\mathcal{H}}_{v}/\text{d}t$ .

4.1.1 Net helicity input

The net helicity flux is given by

(4.2) $$\begin{eqnarray}\frac{\text{d}H}{\text{d}t}=\int _{P}\frac{\text{d}{\mathcal{H}}}{\text{d}t}(\boldsymbol{a}_{0})\,\text{d}x\,\text{d}y,\end{eqnarray}$$

where the integration is taken over all $\boldsymbol{a}_{0}\in P$ . The total helicity input over a period $[t_{0},t]$ through $P$ is then given by

(4.3) $$\begin{eqnarray}H(t)=\int _{t_{0}}^{t}\frac{\text{d}H}{\text{d}t^{\prime }}\,\text{d}t^{\prime }.\end{eqnarray}$$

In this study we choose $t_{0}$ to be the time at which the emerging magnetic field first reaches $z=0$ . The equivalent moving photosphere calculations will be labelled $\text{d}H_{v}/\text{d}t$ and $H_{v}(t)$ . We point out that whilst the calculation concerns a moving surface, the quantity $\text{d}H/\text{d}s$ is calculated on a fixed domain $P$ by projection (see 2.10). This applies to all calculations which consider the moving domain.

4.2 Winding input rate

The field line winding input is

(4.4) $$\begin{eqnarray}\frac{\text{d}{\mathcal{L}}}{\text{d}t}(\boldsymbol{a}_{0})=-\frac{1}{2\unicode[STIX]{x03C0}}\unicode[STIX]{x1D70E}_{z}(\boldsymbol{a}_{0})\int _{P}\unicode[STIX]{x1D70E}_{z}(\boldsymbol{a})\frac{\text{d}\unicode[STIX]{x1D703}(\boldsymbol{a}_{0},\boldsymbol{a})}{\text{d}t}\,\text{d}x\,\text{d}y,\end{eqnarray}$$

and the integral of this over $P$ gives $\text{d}L/\text{d}t$ . The net winding input over a time period $[t_{0},t]$ is

(4.5) $$\begin{eqnarray}L(t)=\int _{t_{0}}^{t}\frac{\text{d}L}{\text{d}t^{\prime }}\,\text{d}t^{\prime }.\end{eqnarray}$$

The moving photosphere calculations will be labelled with a $v$ subscript as for the helicity inputs, with analogous changes to the calculations.

4.3 Weak field corrections

In our simulations, there is no magnetic field outside the emerging regions. Therefore, there is a thin current sheet around the emerging field where the field strength rapidly decreases to zero. In this region the, field line topology is incoherent and not likely to be resolved numerically (typical field strengths in this layer are measured to be ${<}0.01\,\%$ of the peak field strength). Such regions make a negligible contribution to the helicity, due to the field strength weighting, but can have a significant effect on the winding calculations. To remove such potentially spurious information we will, in some calculations, modify the definition of the indicator function $\unicode[STIX]{x1D70E}_{z}$ to be

(4.6) $$\begin{eqnarray}\unicode[STIX]{x1D70E}_{z}(\boldsymbol{x})=\left\{\begin{array}{@{}cc@{}}1 & \text{if }B_{z}(\boldsymbol{x})>0\text{ and }|\boldsymbol{B}|>\unicode[STIX]{x1D716},\\ -1 & \text{if }B_{z}(\boldsymbol{x})<0\text{ and }|\boldsymbol{B}|>\unicode[STIX]{x1D716},\\ 0 & \text{if }B_{z}(\boldsymbol{x})=0\text{ or }|\boldsymbol{B}|\leqslant \unicode[STIX]{x1D716}.\end{array}\right.\end{eqnarray}$$

Where we use this definition in the article, we will be explicit and specify the value of $\unicode[STIX]{x1D716}$ . Otherwise it is to be assumed that $\unicode[STIX]{x1D70E}_{z}$ has its usual definition (2.8).

4.4 Winding and helicity ratios

Both the helicity and winding densities can be positive and negative corresponding to right- and left-handed field entanglement. For mixed helicity flux ropes there can be significant entanglement, locally, but little overall average. This is true of the mixed helicity fields $\boldsymbol{B}_{b}$ which have zero total helicity. To quantify whether the helicity input for a given simulation is biased to one sign or is mixed, we calculate the following ratios

(4.7) $$\begin{eqnarray}\displaystyle H_{t}^{r}=\frac{\text{d}H}{\text{d}t}\Big/\frac{\text{d}H^{abs}}{\text{d}t},\quad \frac{\text{d}H^{abs}}{\text{d}t}=\frac{1}{2\unicode[STIX]{x03C0}}\int _{P}\int _{P}\left|B_{z}(\boldsymbol{a}_{1})B_{z}(\boldsymbol{a}_{2})\frac{\text{d}\unicode[STIX]{x1D703}_{12}(\boldsymbol{a}_{1},\boldsymbol{a}_{2},t)}{\text{d}t}\right|\text{d}x_{1}\,\text{d}y_{1}\,\text{d}x_{2}\,\text{d}y_{2}, & & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

(4.8) $$\begin{eqnarray}L_{t}^{r}=\frac{\text{d}L}{\text{d}t}\Big/\frac{\text{d}L^{abs}}{\text{d}t},\quad \frac{\text{d}L^{abs}}{\text{d}t}=\frac{1}{2\unicode[STIX]{x03C0}}\int _{P}\int _{P}\left|\unicode[STIX]{x1D70E}_{z}(\boldsymbol{a}_{1})\unicode[STIX]{x1D70E}_{z}(\boldsymbol{a}_{2})\frac{\text{d}\unicode[STIX]{x1D703}_{12}(\boldsymbol{a}_{1},\boldsymbol{a}_{2},t)}{\text{d}t}\right|\text{d}x_{1}\,\text{d}y_{1}\,\text{d}x_{2}\,\text{d}y_{2}.\end{eqnarray}$$

4.5 Weighted velocity flux ${\mathcal{V}}_{z}$

The plasma flow across the photosphere will transpire to be an important quantity in what follows. We calculate the net rate of plasma flow through the photosphere weighted by the absolute plasma flow rate,

(4.9) $$\begin{eqnarray}{\mathcal{V}}_{z}=\frac{\int _{P}u_{z}\,\text{d}x\,\text{d}y}{\int _{P}|u_{z}|\text{d}x\,\text{d}y}.\end{eqnarray}$$

5 Simulation analysis

5.1 Twisted flux emergence

First, we analyse the helicity input of a twisted field with a non-dimensional twist of $\unicode[STIX]{x1D719}=-0.4$ and a negative (left-handed) chirality (the other parameters are as stated in § 3.2). This value of twist is a common choice for flux emergence simulations (e.g. Hood et al. 2012, and references therein) and represents a level of twist slightly below what would be required for the onset of the kink instability.

5.1.1 Physical characteristics of emergence

Figure 2. Illustrative distributions characterizing the emergence of the twisted flux rope into the Sun’s atmosphere. (a,d,g) are contour plots of the current density, showing the buoyancy instability-triggered rise and expansion into the corona. Also indicated is the plane $z=0$ . (b,e,h) represent subsets of the field lines at these times. (c,f,i) are the corresponding magnetograms with a clear bipole structure.

Figure 3. Photospheric helicity and winding input for the emergence of the twisted field. (a) The helicity input rate $\text{d}H/\text{d}t$ . It is negative on average and has a consistent oscillation. The vertical lines indicate the times $t=67$ –72 at which the field’s helicity distribution is analysed in what follows. (b) The net helicity input $H(t)$ which is always negative and increases in magnitude over time. (c) The winding input rate $\text{d}L/\text{d}t$ and (d) the total winding input $L(t)$ .

A description of the emergence of twisted flux tubes has been described in detail in many other studies (e.g. MacTaggart & Hood 2009; Hood et al. 2012) so we only highlight aspects critical to the subsequent helicity input analysis. Illustrative figures are shown in figure 2 at times $t=35,55$ and $71$ . The magnetic field reaches the photosphere where it remains until the plasma $\unicode[STIX]{x1D6FD}$ drops to order unity and the field becomes subject to the magnetic buoyancy instability, subsequently expanding into the corona (Hood et al. 2012). The initial toroidal flux tube leads to a bipole emergence which is clearly visible in the magnetograms.

Figure 4. Time series of the ratios (a) $H_{t}^{r}$ and (b) $L_{t}^{r}$ for the twisted field emergence.

5.1.2 Helicity and winding input time series

The temporal evolutions of both $\text{d}H/\text{d}t$ and $H$ are shown in figures 3(a) and 3(b) respectively. The tracking of the helicity input begins when the magnetic field reaches the photospheric boundary which, in this simulation, is almost coincident with the onset of the magnetic buoyancy instability ( $t=18$ ). With the initial expansion of the field into the corona there is an input of negative helicity (figure 3 a), as is to be expected for this negative helicity flux rope (MacTaggart & Hood 2009). Later, the helicity input settles into an oscillatory pattern with the input rate, occasionally, becoming net positive. The net helicity input $H$ , shown in figure 3(b), exhibits, on average, an almost linear input increase of negative helicity (in line with the results of Sturrock et al. (2015) accounting for the opposing sign of twist with a smaller oscillation about this trend). The time series $\text{d}L/\text{d}t$ and $L$ , shown in figure 3(c,d), are qualitatively similar to their helicity counterparts so we do not focus on them in what immediately follows. In contrast to these results, the oscillations in the $\text{d}H/\text{d}t$ time series calculated in Sturrock et al. (2015) are relatively small (though similarly coherent) and do not cause the change in sign we find here. We will expand more on this comparison later. Finally, the ratios $H_{t}^{r}$ and $L_{t}^{r}$ , whose time series are shown in figure 4, indicate, as expected, a preference towards one sign of helicity input, i.e. the quantity is often close to 1 in magnitude so that most of the winding/helicity density is coherent in sign (the sign, of course, oscillates from positive to negative). An important question to address is what is causing the oscillations and, further, why are they so significant in comparison to previous studies?

Figure 5. Helicity input rate density distributions $\text{d}{\mathcal{H}}/\text{d}t$ of the emerged region during the period $t\in [67,72]$ over which the sign of $\text{d}H/\text{d}t$ varies form positive to negative (net). The plots shown in (a–f) correspond to times $t=67$ – $72$ indicated by the set of vertical lines on figure 3(a). (a) The density is dominantly positive around the PIL. In (b–f), the helicity input at the PIL changes from dominantly positive to negative. The helicity input in the two flux poles becomes dominantly negative over the cycle.

Figure 6. Magnetic and velocity field distributions in the photosphere at $t=67$ . (a) The transverse field $\boldsymbol{B}_{\Vert }$ superimposed on a scalar plot of its magnitude (stronger red colouring implies a stronger field). The twist at the positive pole is left-handed and the twist at the negative pole is right-handed. (b) The transverse velocity field $\boldsymbol{u}_{\Vert }$ superimposed on a scalar plot of its magnitude (stronger blue colouring implies a stronger field). Both vortices have right-handed rotation.

Figure 7. Distributions of the components of the helicity producing field $\boldsymbol{w}$ , at times $t=67$ and $72$ . Panels (a,b) are the distributions of $\sqrt{w_{x}(\boldsymbol{a}_{0})}$ and $\sqrt{w_{y}(\boldsymbol{a}_{0})}$ respectively at $t=67$ . Panels (c,d) are the same distributions but at $t=72$ . For both $w_{x}(\boldsymbol{a}_{0})$ and $w_{y}(\boldsymbol{a}_{0})$ , the sign of the distribution either side of the PIL reverses.

We now focus on the time period $t\in [67,72]$ , indicated in figure 3(a) by a set of vertical lines. There is a variation between positive and negative input rates over this period. There is no particular reason why we choose to present the results of this period over the rest of the input cycle (except the rise stage). Indeed we checked that the following analysis would have led to similar conclusions for the behaviour of the system over the period $t\in [40,76]$ .

Distributions of $\text{d}{\mathcal{H}}/\text{d}t$ , during the considered time interval, are shown in figure 5. At $t=67$ (a) and $t=72$ (f) there are respectively strong positive and then negative field line helicity input densities around the polarity inversion line (PIL). For the intermediate times (b–d), the field line helicity density in this region shows a gradual variation from positive to negative. Patches of strong negative density at the two poles of flux distribution also develop over the period. There are fluid vortices, in the in-plane velocity field, which are centred on these poles. It is shown in figure 6 that at the magnetic field’s positive pole, the vortex opposes the direction of magnetic twist, whilst at the negative pole the signs of the vortical and magnetic twists agree. We might speculate that some kind of relative balance between the twisted field input and fluid rotation might lead to the oscillation in helicity input around the inversion line. However, we now show this oscillation occurs instead due to the cyclic submergence and re-emergence of the flux, that is, the movement of part of the magnetic field below and above the photospheric plane. In what follows, submergence does not necessarily indicate that the field passes deep beneath the photosphere. As will be demonstrated, any part of the field that makes contact with the photospheric boundary, and perhaps passing only slightly beneath it, will register a signal in the magnetograms and in the helicity and winding inputs. Therefore, in this work, submergence is any previously emerged field moving down from the atmosphere or otherwise making contact with the photosphere from above.

5.2 Helicity sign change around the PIL

Figure 8. Plots of the winding input rate over the period $t=67$ to $71$ for (a) the unrestricted calculation $\text{d}L/\text{d}t$ and (b) the capped version $\text{d}L^{c}/\text{d}t$ for which the threshold is $c=2$ .

The distributions of the $x$ and $y$ -components of $\boldsymbol{w}=\text{d}\boldsymbol{a}/\text{d}t$ at $t=67$ and $t=72$ are shown in figure 7. The most significant helicity producing velocities arise in the neighbourhood of the PIL (the plots are of $\sqrt{w}_{x}$ and $\sqrt{w}_{y}$ are shown for clarity, thus relatively exaggerating the magnitudes in weaker regions). There is a clear switch in sign of both the $w_{x}$ and $w_{y}$ densities on both sides of the PIL over this cycle. To test if it is the temporal change in the sign of the field line velocity field $\boldsymbol{w}$ close to the PIL which determines the change in sign of the helicity input, we consider a modified vector field $\boldsymbol{w}^{c}$ as follows

(5.1) $$\begin{eqnarray}\boldsymbol{w}^{c}=\left\{\begin{array}{@{}cc@{}}\boldsymbol{w} & \text{ if }\left|\boldsymbol{w}\right|\leqslant c,\\ 0 & \text{ if }\left|\boldsymbol{w}\right|>c,\end{array}\right.\end{eqnarray}$$

where $c$ is some chosen threshold. We choose $c=2$ here, which represents a cutoff speed much smaller than typical field line speeds at the PIL ( ${\approx}10$ ). As long as the cutoff is small enough, the following results are robust. The main behaviour of the following results was also obtained for $c=1$ and $c=5$ .

In essence, we cut out the higher helicity producing velocities which reside primarily in the region of the PIL and not at the centres of the magnetic footpoints. We also calculate the modified field line winding input rate $\text{d}L^{c}/\text{d}t$ using $\boldsymbol{w}^{c}$ . We see in figure 8(a) that the full winding input rate $\text{d}L/\text{d}t$ changes sign from positive to negative over the time period $t\in [67,72]$ . By contrast, for the restricted input $\text{d}L^{c}/\text{d}t$ (b), the sign of winding input is always negative over this period and two orders of magnitude smaller. Therefore, ignoring the winding input due to the field line motions $\boldsymbol{w}$ around the PIL leads to a time series which does not have a positive input rate over this period (the story is the same for the helicity input rate). It was confirmed that this is true throughout the simulation (for the time period when the oscillations in $\text{d}H/\text{d}t$ occur).

Figure 9. Distributions of the vertical velocity $u_{z}(\boldsymbol{a}_{0})$ at times (a) $t=67$ and (b) $72$ . (c) A plot of the net velocity flux ${\mathcal{V}}_{z}$ through $P$ . (d) Plots of $\text{d}H/\text{d}t$ and ${\mathcal{V}}_{z}$ scaled to have values between $0$ and $1$ .

The cause of the helicity sign change is found to be related to a switch in the sign of $u_{z}$ in the region around the PIL. To establish this, we first evaluate the relative magnitudes of the terms $\boldsymbol{u}_{\Vert }$ and $\boldsymbol{B}_{\Vert }u_{z}/B_{z}$ . The submergence/emergence term $\boldsymbol{B}_{\Vert }u_{z}/B_{z}$ is typically found to be an order of magnitude higher in the regions where $\boldsymbol{w}$ is (relatively) large. The $u_{z}$ distributions at the start ( $t=67$ ) and end ( $t=72$ ) of the cycle are shown in figure 9(a,b). Initially, there is a net negative velocity with most of the flow surrounding the PIL. Then at the end of the period the velocity is net positive with most of the flow being at the PIL. A time series plot of the net velocity flux ${\mathcal{V}}_{z}$ is shown in figure 9(c). A scaled comparison of the maxima and minima of $\text{d}H/\text{d}t$ and ${\mathcal{V}}_{z}$ is shown in figure 9(d) emphasizing that it is motion across the photosphere that is dominating the helicity in this phase of emergence.

5.2.1 Flux rope centre at the photosphere

Figure 10. Arrows indicating the magnetic field in the plane orthogonal to the PIL at $(0,0,0)$ at $t=67$ . This is superimposed on a plot of the current magnitude in the same plane. The core of the flux rope is clearly visible and is centred at the PIL.

The fact that the oscillations in helicity (and winding) are sufficient to lead to a sign of input which opposes the field’s chirality is a result of the fact that the bulk of the field’s initially twisted flux rope remains trapped at the photosphere, see figure 10. The plasma $\unicode[STIX]{x1D6FD}$ $(=2p/|\boldsymbol{B}|^{2})$ has a value of approximately 5 at the photosphere in this simulation. This implies that the magnetic field is not dynamically dominant (as in the low- $\unicode[STIX]{x1D6FD}$ corona) and so can moved by the surrounding plasma. A combination of upward motion, from emergence, and downward motion, from draining plasma (e.g. Hood et al. 2012) causes the flux rope centre (axis) to oscillate about the $z=0$ plane. In this topologically simple field (relative to mixed helicity model that we will study shortly) the bulk of the topological information is concentrated at the flux rope centre. Therefore, if the centre crosses the photospheric plane, the response in the helicity and winding rates is large. As mentioned before, it is transport across the photospheric plane and not (un)twisting motions on the plane that causes the largest changes in helicity and winding.

If the plasma $\unicode[STIX]{x1D6FD}$ were smaller, it would be expected that the oscillations would be less pronounced. In Sturrock et al. (2015), who perform a very similar simulation but with a flux rope of almost double the initial field strength to the one we consider, they still find oscillations in $\text{d}H/\text{d}t$ but much less pronounced than those found here (their oscillations do not change the sign of the helicity input rate). Two factors are important in contributing to this change in behaviour. The first is that the axes of flux tubes with higher field strengths can emerge further in the atmosphere compared to those with weaker field strengths (MacTaggart & Hood 2009). This makes the flux rope centre further away from the $z=0$ plane and so there is less flux crossing the plane in any partial submergence event. Secondly, the plasma $\unicode[STIX]{x1D6FD}$ is smaller and the magnetic field is less susceptible to surrounding plasma motions.

5.3 Tracking the moving photosphere

Figure 11. Photospheric helicity and winding inputs for the emergence of the twisted field which account for the tracking of the changing photosphere geometry. (a) The helicity input rate $\text{d}H_{v}/\text{d}t$ . (b) The net helicity input $H_{v}(t)$ . (c) The winding input rate $\text{d}L_{v}/\text{d}t$ . The vertical lines are at times $t=30,31$ when the field structure is analysed further in what follows. (d) The total winding input $L_{v}(t)$ .

Figure 12. Vertical velocity distributions at $t=28$ . (a) The velocity density $u_{z}(P_{v}(x,y))J$ at the moving photosphere line $\unicode[STIX]{x1D70C}=1$ . (b) The difference between the varying photosphere velocity density $u_{z}(P_{v}(x,y))J$ and the $z=0$ velocity density $u_{z}(x,y)$ .

The velocity flux oscillations indicated in figure 9 imply that the surface $\unicode[STIX]{x1D70C}=1$ will be varying in space and time. As discussed in the introduction, we can track this variation and calculate modified quantities which account for this motion. The time series of the adjusted quantities $\text{d}H_{v}/\text{d}t$ , $H_{v}$ , $\text{d}L_{v}/\text{d}t$ and $L_{v}$ are shown in figure 11. On a qualitative level, the helicity time series (panels (a,b) respectively) are effectively the same as the $z=0$ calculations shown in figure 3. The magnitudes, however, of the various peaks of the input rate $\text{d}H_{v}/\text{d}t$ are roughly four times larger than those of $\text{d}H/\text{d}t$ . This magnification is due to increased variations in velocities at the $\unicode[STIX]{x1D70C}=1$ surface. As an example, figure 12 displays the varying photosphere velocity distribution $Ju_{z}(P_{v}(x,y))$ in (a) and the difference $Ju_{z}(P_{v}(x,y))-u_{z}(x,y)$ between the moving surface and static surface velocity distributions in (b). In (a), the strongest velocities are, as before, at the PIL. In (b), major differences are at the main footpoints. We note that the differences can be greater in magnitude than a typical velocity $Ju_{z}(P_{v}(x,y))$ . Therefore, it is not simply a case of larger velocity magnitudes at the $\unicode[STIX]{x1D70C}=1$ surface compared to the $z=0$ plane, rather a difference in the distribution of flows on these surfaces.

Figure 13. A slice of the vector field in the plane orthogonal to the PIL at $t=30$ , just prior to the spike in the time series $\text{d}L_{v}/\text{d}t$ shown in figure 11(c). Also shown as a green curve is the intersection of the surface $P_{v}$ and this plane. The background distribution is the out of plane component of the current density.

Figure 14. Magnified representations of a slice of the vector field in the plane orthogonal to the PIL at $t=30$ and $t=31$ , prior to and at the spike in the time series $\text{d}L_{v}/\text{d}t$ shown in figure 11(d). Shown as a green curve is the intersection of the surface $P_{v}$ and this plane. The background distribution is the out of plane component of the current density. (a) shows ( $t=30$ ) a magnification of the region of figure 13 which contains the curve $\unicode[STIX]{x1D70C}=1$ . (b) shows ( $t=31$ ) the same distribution as in (a) but at the time of the spike. The $\unicode[STIX]{x1D70C}=1$ curve has dropped vertically by a value of approximately $0.5$ .

The winding time series, panels (c,d) of figure 9, show significant qualitative differences compared to the $z=0$ input calculations shown in panels (c,d) of figure 3. There is one significant period of negative input rate $\text{d}L_{v}/\text{d}t$ at $t=31$ . This is near the end of the period when the central part of the field (containing the original tube axis) emerges at the photosphere (see figure 13). The spike coincides with a relatively sharp (negative) change in the height of the $\unicode[STIX]{x1D70C}=1$ surface from $t=30$ (figure 14 a) to $t=31$ (figure 14 b). This change results in a significant change in the negative helicity input as the flux rope core moves further beyond the photospheric boundary. It was confirmed that this jump in the moving photospheric surface was atypically large from $t=30$ to $t=31$ by comparison to typical changes in its morphology over single time step of size $1$ . After this, the input rate $\text{d}L_{v}/\text{d}t$ falls to a relatively small rate. It is interesting to note that this occurs due to a relatively fast change in the geometry of the $\unicode[STIX]{x1D70C}=1$ surface, rather than a sudden change in field topology, indicating that this may be an important factor in interpreting observational inputs of field topology. The oscillations shown in the helicity input rate $\text{d}H_{v}/\text{d}t$ are still present but, as shown in figure 11(d), these have a small effect on the net helicity injected into the photosphere, by comparison to the sharp negative input around $t=30$ .

It is intriguing that the (moving photosphere) winding input more clearly represents the transition between the emergence of the field’s twisted core to the photosphere (the $\unicode[STIX]{x1D70C}=1$ surface) followed by the lack of changing topological input due to the core getting stuck at this surface. The significantly reduced (at least on a relative scale) size of the oscillations in comparison to the helicity input results from the lack of field strength weighting in the winding input. The inclusion of field strength in the helicity input rate magnifies the oscillatory signal due to the central core of the field being immediately surrounded by strong field – slight dips under the photospheric boundary produce high signals because of the strong field strengths. Another way to view this result is that the winding input treats the oscillation at the photosphere as not especially interesting at later times as it is not really indicating the input of new topological information into the atmosphere.

5.4 Mixed helicity emergence

We now consider a mixed helicity field $\boldsymbol{B}_{b}$ , specified by (3.23) with $n=2$ . The emergence of mixed helicity (braided) magnetic fields and their associated dynamics was first described in Prior & MacTaggart (2016).

5.4.1 Physical characteristics of the emergence

Figure 15. Illustrative distributions characterizing the emergence of a mixed helicity flux rope with $n=2$ twist pairs. (a,d,g) are magnetograms. (a) depicts ( $t=27$ ) a bipole structure associated with the initial emergence phase. In (d) ( $t=38$ ), this distribution has separated slightly and there also appear to be some weak thin horizontal structures appearing at the centre of the domain. In (g) ( $t=62$ ), these weak additional structures have developed into pairs of bipoles whose polarity oppose that of the larger initial bipole. (b,e,h) are current contours ( $|\boldsymbol{j}|=0.01$ ) which indicate the field’s expansion. Also indicated is the plane $z=0$ . (c,f,i) are representative field lines.

Various illustrative visualizations of the emerging field’s evolution are shown in figure 15. As for the twisted tube, the emerging field gets trapped below the photosphere before becoming subject to the magnetic buoyancy instability and rising into the coronal region. The magnetograms (a,d,g) show that, in addition to a larger bipole structure developing, there is also an additional formation of smaller bipole structures in between these larger poles with orientations in opposition to that of the large bipole. In (b,e,h), contours of the current density structure show the field’s expansion, which is restricted along the PIL of the main bipole. This behaviour is in contrast to the single structure found in the twisted case (figure 2) which shows a more uniform expansion. The field lines, shown in  (c,f,i) indicate a pair of intertwined field line arcades linked by a centrally dipped section (which partially submerges, as we will see later). There is no significant twisting in the field lines. It was established in Prior & MacTaggart (2016) that the dips in the current contours coincide with plasma draining and are the nonlinear manifestations of a Rayleigh–Taylor instability.

Figure 16. Photospheric helicity and winding inputs for the emergence of the mixed helicity field. (a) The helicity input rate $\text{d}H/\text{d}t$ . The input is both positive and negative. Also shown are lines at $t=25,42$ and $48$ . (b) The net helicity input $H(t)$ . (c) The winding rate $\text{d}L/\text{d}t$ . (d) The total winding input $L(t)$ , the sign of this input is generally opposite to that of the helicity.

5.4.2 Helicity and winding input: static photosphere

Plots of the time series of the quantities $\text{d}H/\text{d}t,\,H,\,\text{d}L/\text{d}t$ and $L$ are shown in figure 16. The helicity rate $\text{d}H/\text{d}t$ (a) is both positive and negative and shows oscillatory behaviour (less temporally regular than for the twisted field). The net helicity input $H(t)$ (b) switches from negative to net positive. The winding input rate $\text{d}L/\text{d}t$ (c) shows much more rapid variations. The net input $L$ (d) has a pattern qualitatively similar to the net helicity but with opposing sign. The ratios $H_{t}^{r}$ and $L_{t}^{r}$ are found to be typically 1–2 % (figure 17), indicating the topological input is not far off neutral, as to be expected.

The distributions $\text{d}{\mathcal{H}}/\text{d}t$ and $\text{d}{\mathcal{L}}/\text{d}t$ are shown in figure 18. The snapshots are chosen to display critical characteristics of these distributions which are present at various times of the emergence process.

There are always four distinct regions of substantial positive and negative helicity (and winding) and these occur in pairs. These regions are generally of similar size and magnitude, which explains why the net helicity (and winding) input is much smaller than its absolute total. In practice, it is difficult to detect the imbalances in the maps of $\text{d}{\mathcal{H}}/\text{d}t$ that lead to the time series variations shown in figure 16(a). By contrast, the winding distributions $\text{d}{\mathcal{L}}/\text{d}t$ capture more features related to emergence and submergence (we will return to this point later). We ignore the distribution values around the edge of the magnetic domain as the field is very weak there. The field line structure in the boundary current sheet is incoherent and, thus, any apparent winding input is not meaningful.

Figure 17. Time series of the ratios (a) $H_{t}^{r}$ and (b) $L_{t}^{r}$ for the $n=2$ mixed helicity field emergence.

Figure 18. Distributions of $\text{d}{\mathcal{H}}/\text{d}t$ and $\text{d}{\mathcal{L}}/\text{d}t$ which indicate the various characteristic stages of the distribution evolutions (indicated with vertical lines on figure 16). (a) $\text{d}{\mathcal{H}}/\text{d}t$ and (b) $\text{d}{\mathcal{L}}/\text{d}t$ at $t=25$ . In both cases there are two significant regions of positive input and two with negative input, a quadrupolar distribution. In (b) the straight PIL is visible along the $y$ axis. (c) $\text{d}{\mathcal{H}}/\text{d}t$ and (d) $\text{d}{\mathcal{L}}/\text{d}t$ at $t=42$ . In (c) the four regions from (a) remain but are surrounded by thin strips of helicity input of opposing sign. In (d), by contrast, the sign of the four regions have swapped compared to (b) and there are additional sub regions of positive and negative input centred on the PIL. (e) $\text{d}{\mathcal{H}}/\text{d}t$ and (f) $\text{d}{\mathcal{L}}/\text{d}t$ at $t=48$ . The helicity input (e) still has four dominant domains but they have swapped sign from the distributions in (a,c). In (f) the flux region seen in the latter magnetograms (figure 15) appears as a distorted PIL and there are strong concentrations of helicity either side of it.

The PIL is clear in all $\text{d}{\mathcal{L}}/\text{d}t$ distributions and its shape varies significantly as various new helicity structures appear (note that all additional structures appear at the PIL).

5.5 Helicity and winding input: moving photosphere

Figure 19. Helicity and winding inputs for the emergence of the mixed helicity field, allowing for the varying photospheric geometry. (a) The helicity input rate $\text{d}H_{v}/\text{d}t$ . (b) The net helicity input $H_{v}(t)$ . (c) The winding rate $\text{d}L_{v}/\text{d}t$ . The vertical lines mark times $t=45,46,47$ at which the field is analysed in detail in the text. (d) The total winding input $L_{v}(t)$ .

Time series of $\text{d}H_{v}/\text{d}t$ , $H_{v}$ , $\text{d}L_{v}/\text{d}t$ and $L_{v}$ , which account for the varying photospheric geometry, are shown in figure 19. In order to prevent the winding measure being affected by any weak (unresolved) field, we utilize the modified definition of the function $\unicode[STIX]{x1D70E}_{z}$ (4.6) with $\unicode[STIX]{x1D716}=0.0001$ for the $\text{d}L_{v}/\text{d}t$ and $L_{v}$ calculations. As was the case for the twisted field, the only significant change in the helicity time series compared to those shown in figure 16(a,b) is in regard to the magnitudes. This difference could, again, be traced to differences in the vertical velocity distribution. The time series of the winding quantities $\text{d}L_{v}/\text{d}t$ and $L_{v}$ , shown in figure 19(c,d), differ significantly from those shown in figure 16(c,d). We find that this is, in part, due to the varying photosphere $P_{v}$ and partly due to ignoring the contributions from significantly weak regions of the field. As discussed previously, the winding is more sensitive to large input from topologically complex field lines of very weak field strength. Hence, the calculation of the winding requires some care.

The sign of the net input $L_{v}$ is largely dictated by a significant spike in the input $\text{d}L_{v}/\text{d}t$ at $t=47$ . After this, there are two large oscillations which are more balanced. The net inputs $H_{v}$ and $L_{v}$ now agree in terms of the net sign of their input over the simulation (as opposed to the quantities $H$ and $L$ , emphasizing the sensitivity of $L$ or $L_{v}$ ).

5.6 Field submergence and deformation

As mentioned above, a clearer transition is revealed by the winding compared to the helicity. The key question to answer here is, what is causing this transition? Some indications can be found in the distributions of $\text{d}{\mathcal{L}}_{v}/\text{d}t$ shown in figure 20. These distributions are at the times indicated by the parallel lines in figure 19(c). Figure 20(a) displays $\text{d}{\mathcal{L}}_{v}/\text{d}t$ when $\text{d}L_{v}/\text{d}t$ is changing from positive to negative. A small region of negative $\text{d}{\mathcal{L}}_{v}/\text{d}t$ is visible at the centre (0,0). The region grows in (b), which corresponds to when $\text{d}L_{v}/\text{d}t$ takes its largest negative value. In (c), which corresponds to the large positive spike in $\text{d}L_{v}/\text{d}t$ and positive jump in $L_{v}$ , the distribution changes to a clear antisymmetrical pattern. There is no longer the isolated patch of negative $\text{d}{\mathcal{L}}_{v}/\text{d}t$ at the centre and the PIL is highly deformed.

Figure 20. Distributions of the winding input density $\text{d}{\mathcal{L}}_{v}/\text{d}t$ covering the period of the time series shown in figure 19(d) at which there is a strong spike in input. (a) At $t=45$ , before the spike. (b) At $t=46$ . (c) At $t=47$ at the spike’s peak.

Figure 21. A contour plot of density $\unicode[STIX]{x1D70C}=0.2$ , slightly above the $\unicode[STIX]{x1D70C}=1$ photospheric level at $t=62$ . There is a an $s$ -shaped peak in the density profile at the centre of the domain which is parallel to the $y$ -direction and has been established to correspond to a pooling of dense plasma in this location. This shape matches the morphology of line observed in the winding input densities $\text{d}L/\text{d}t$ , shown in figure 18(f), and $\text{d}L_{v}/\text{d}t$ , shown in figure 20(c).

Figure 22. Slices in the $y$ – $z$ plane, at the PIL at (a) $t=27$ and (b) $t=42$ , revealing different phases of emergence. The slices display distributions of the out of plane component of the current density. Also shown, as a green line, the $\unicode[STIX]{x1D70C}=1$ surface of the moving photosphere. (a) ( $t=27$ ) There is a spiral structure, indicative of a locally twisted field just off the centre of the photospheric domain. It is part way through emerging through the $\unicode[STIX]{x1D70C}=1$ photosphere line. (b) ( $t=42$ ) The previous twist has split in two due to the deformation of the current structure in the atmosphere.

Figure 23. Slices in the $y$ – $z$ plane at the PIL at (a) $t=46$ and (b) $t=47$ . The slices display distributions of the out of plane component of the current density. Also shown, as a green line, the $\unicode[STIX]{x1D70C}=1$ surface of the moving photosphere.

Since the emerging field has very little twist, it cannot easily support dense plasma in the atmosphere. This manifests itself in the atmosphere as the buckling of the central part of the emerging field, where dense plasma drains and restricts the magnetic field in the lower atmosphere. This process is visualized in figure 15 and the pooling of dense plasma which forms can be seen in figure 21. The patch of negative $\text{d}{\mathcal{L}}_{v}/\text{d}t$ , described above, indicates the beginning of motion leading to submergence, that is, the movement of helicity carrying field from the atmosphere down to the photospheric boundary. We now characterize various stages of the emergence up to and including this submergence event.

Figure 22(a) displays the emergence of the mixed helicity field at the (relatively) early time of $t=27$ . The slice at $x=0$ displays the $x$ -component of current density, the projection of magnetic field arrows on the plane and the position of the moving photosphere. What is shown is a twisted structure beginning to emerge into the atmosphere. Comparison with figure 16 shows that the helicity and winding inputs are still relatively weak. Figure 22(b) displays a similar slice but for the much later time of $t=42$ . Comparison with figure 16 shows that this corresponds to the helicity approaching a local (positive) maximum and the winding near its most negative value. A current structure exits in the atmosphere in figure 22(b) that was not present in (a). This current developed due to the buckling of the field due to dense plasma drainage, as described previously. This nonlinear deformation of the field has resulted in two ‘twist units’, as opposed to the single unit displayed in (a). Where the two twist units meet, at $y=0$ and just above the photosphere, there is a small patch of negative current. This patch corresponds to the patch of negative winding displayed in figure 20(a,b), which reveals this patch growing at the later times of $t=45$ and $t=46$ .

As mentioned above, figure 20(c) reveals that the patch of negative winding found in figure 20(a,b) disappears. Instead, the winding distribution reveals a highly deformed PIL and winding magnitudes double that of the previous time step. This behaviour is caused by submergence. Consider figure 23. This figure shows $y$ – $z$ slices (at $x=0$ ) of the $x$ -component of the current density and the photospheric boundary at times (a) $t=46$ and (b) $t=47$ , corresponding to the times of the winding distributions in figure 20(b,c). In figure 23(a) the structure of positive current, created by dipped field in the atmosphere, is just glancing the photospheric boundary. In figure 23(b), this structure has passed slightly beneath the photospheric boundary. Despite the appearance of figure 23(a,b) being only marginally different, the submergence of a very small part of the field can have a large impact on the winding distributions (as well as the magnetograms). As submergence continues, the changing current structure leads to a large spike in the winding input (figure 19 c,d).

It is interesting to note that the original emerging field (such as the twisted unit in figure 22 a) did not lead to a spike in the winding input. It was only after the emerged field above the photosphere changed and this new structure was submerged that an event was detected in the time series.

Further confirmation that submergence dominates these helicity and winding time series is found, as we did for the twisted case, by examining the velocity flux. Time series of the velocity flux and the $\text{d}H/\text{d}t$ are compared in figure 24(a). After $t=40$ , ${\mathcal{V}}_{z}<0$ for the period corresponding to the draining and submergence described above. We also note that there is a correlation between oscillations in the helicity input time series and the velocity flux (figure 24 b).

Figure 24. Plots of the velocity flux ${\mathcal{V}}_{z}$ and its correlation to the helicity input rate. (a) The temporal variation of ${\mathcal{V}}_{z}$ . (b) A comparison of the temporal variation of the quantities ${\mathcal{V}}_{z}$ and the helicity input rate $\text{d}H/\text{d}t$ (both scaled to lie in the range $[0,1]$ ).

6 Conclusions and discussion

This article presents simulations of the emergence of magnetic flux tubes with varying internal geometries and topologies. One simulation considers a twisted tube which has a single sign helicity (negative in this case). The other simulation focusses on a tube of mixed helicity. The two main quantities analysed are the (relative) helicity $H$ and the winding $L$ , both found by considering their rate of change through the photosphere. We have performed calculations through a planar surface, representing the initial position of the photosphere in the simulations, as well as a surface of constant density which represents a (potentially) changing photosphere. Both quantities ( $H$ and $L$ ) rely on the same baseline information, the changing entanglement of the points at which magnetic field lines pierce the photosphere, i.e. the evolving field topology affected by both the in-plane fluid motion and the vertical motion of the field passing through the photosphere. The major difference between $H$ and $L$ is the weighting of magnetic flux for the helicity input. The following observations emerge from the analysis:

1. (i) Large oscillations in the helicity input are consistently related to the reversal of plasma flow across the photosphere, occurring in both the single sign and mixed helicity emergence simulations. These reversals appear in the spatial distributions of the helicity and winding rates as the temporally varying sign of clear structures. These structures can be further correlated to the appearance of flux structures in the magnetogram distributions, leading to several features in the magnetograms being linked to submergence rather than emergence.

2. (ii) The winding time series tends to be far more sensitive to the emergence (submergence) of field topology and can detect individual events of helicity/winding carrying structures intersecting the photospheric boundary. The helicity time series tends to give more gradual variations in input. It is sensitive to the emergence (submergence) of topologically complex field when this is linked to the transport of significant magnetic flux crossing the photospheric boundary.

3. (iii) The ratios of the input of helicity and winding to their unsigned equivalents is far larger for the twisted emergence than for the mixed helicity emergence. This result was to be expected based on the initial conditions: the twisted tube has net helicity, the mixed helicity field does not.

4. (iv) It is important to take the changing geometry of the photosphere into account, particularly for the calculation of the winding. It is also important to remove the effect of weak or unresolved field, which can affect the winding calculation both quantitatively and qualitatively.

Based on our analysis of the simulations, we make the following suggestions which may aid or enhance the analysis of observational flux emergence data:

1. (i) We recommend that both the helicity and winding input distributions and time series should be calculated. This task should be straight forward, in one sense, as the only main difference in the calculation of the winding compared to the helicity is to substitute an indicator function measuring the sign of $B_{z}$ instead of $B_{z}$ itself. Since the winding is more sensitive to noise, however, suitable preprocessing of the observational data will be required. If possible, the varying geometry of the photosphere should be accounted for as it can affect both the magnitude and morphology of the input time series.

2. (ii) We recommend that the behaviour of the helicity and winding rates and distributions should be correlated to the vertical velocity flux, where possible. This would aid in determining if a new structure in a magnetogram is the result of emergence or submergence, a task which would be difficult if reliable information is restricted only to photospheric magnetograms. On that note, we acknowledge that the $\text{d}H/\text{d}t$ and ${\mathcal{V}}_{z}$ time series for the twisted field simulation show a shift between anti-correlation and correlation (figure 9 d), while for the mixed helicity simulation these quantities appear to be more coherently correlated (figure 24 b). Although we have not investigated the relationship between the time series in detail, that is beyond recognizing when magnetic field is emerging or submerging, it is an area worthy of future study and may reveal important characteristics of the emergence of different types of magnetic field.

3. (iii) We recommend a combined analysis of the winding input time series and distributional information. These can be used to detect events in which magnetic field structures with complex topology pass through the photosphere (from above or below). This approach provides more information than magnetogram distributions, which indicate new emergence/submergence events but no information regarding their topological content. However, care must be taken to discount events which have significantly weak magnetic field and may be spurious. Such events may have little significance for the overall evolution of the emerging magnetic field. This point is made clear in the difference between the winding series for the mixed helicity field. In the non-moving surface case the sign of the total helicity and winding inputs differ (figure 16 b,d). It is entirely plausible that this can occur for a field which has significant small-scale but well defined topology (something with a high winding measure) that occurs in a region of weak field ( $B_{z}\ll 1$ ). Then the helicity assigns a very small weight to this topology (winding) due to the multiplication by a very small $B_{z}$ . However, as discussed above, in simulations of emergence into an, initially, magnetically empty region above the photosphere, there is a current sheet surrounding the emerging flux rope. In practice this leads to a small region on the boundary of the emerged field where the field strength rapidly decreases to zero. The field line structure in this boundary layer is not resolved and, therefore, not to be trusted. It was found that by removing these contributions, the helicity and winding inputs have the same overall sign (see figure 19 b,d). A significant part of this change was found to be due to the cutoff in field strength applied to the winding calculations. The size of $\unicode[STIX]{x1D716}$ (used in the cutoff) that needs to be prescribed will depend on the implementation of how the magnetic field is approximated numerically. In particular, it will depend on the numerical scheme (finite difference, finite element, etc.) and the grid resolution. Without rigorous theoretical results, a ‘trial and error’ approach to $\unicode[STIX]{x1D716}$ selection is recommended. If the value is too low, the winding results will be significantly different for slightly higher values of $\unicode[STIX]{x1D716}$ . If the value is too high, field that can be resolved numerically may be removed from the calculation.

4. (iv) We recommend that both the helicity and winding ratios, $H_{t}^{r}$ and $L_{t}^{r}$ , should be calculated. Again, these calculations should be straight forward as they use the same basic information as the helicity and winding inputs. The evidence from our study suggests that a consistently low value of these ratios is indicative of an emerging mixed helicity structure where the internal twisting in the field is not dominated by one sign. Our results suggest that a ratio below $5\,\%$ is indicative of a mixed helicity field. That being said, we have only considered a strongly twisted case. For emerging fields with lower twists, the ratio may also be small and this is something to be considered in future work.

Apart from the points mentioned above, there are several directions in which to carry this work forward. So far, we have considered emerging regions before the formation of any eruptions, i.e. the motion we have studied is effectively ideal. How non-ideal effects (e.g. reconnection) change the helicity and winding signatures, including their interpretation, remains to be studied. It is known that for high magnetic Reynolds number plasmas, dissipative losses to the (relative) helicity in the volume are small (Berger 1984; Pariat et al. 2015; Sturrock et al. 2015). This is true typically of dissipative helicity loss at the photospheric surface (Pariat et al. 2015; Sturrock et al. 2015). However, this does not necessarily imply that local reconnection rates are low (which depend on the electric current), just that the average losses are small. Thus, since there are small regions of relatively high current at the photosphere in these simulations, there may be some regions where local losses are significant. This, in turn, may lead to significant local changes in field topology.

Since our focus in this study was not on the field’s evolution above the photosphere, we did not make direct calculations of the relative helicity in the volume above the photosphere. This calculation has been performed in various flux emergence studies in which the changing structure of the field above the photosphere was of significant focus (e.g. Pariat et al. 2015, 2017; Sturrock et al. 2015). As an example comparison, the helicity input in the twisted case of this manuscript has a linear total input with an oscillatory component. This behaviour is in line with the results of Sturrock et al. (2015) (allowing for the fact that the oscillations in our study are larger due to the field’s twisted core getting stuck at the photosphere), who perform the relative helicity calculation in the volume and confirm, as expected, that the total helicity input is approximately equal to the relative helicity in the volume. Theoretically, as long as the emerging field does not interact with the boundaries of the computational domain and the evolution is ideal, the helicity flux and direct calculations should give identical results. Errors in different numerical implementations and resistive effects can play a role in making these results diverge, so checking both calculations is, generally, a good idea. This time, however, we feel that the above comparison with the results of Sturrock et al. (2015) justifies our decision to postpone such calculation checks for a thorough study of the topology of the field above the photosphere.

Further to this, previous analyses of helicity dissipation have only considered the effects of scalar diffusion. In flux emergence applications, a diffusion tensor, with components parallel and perpendicular to the magnetic field, represents a more accurate diffusion model for the region between the photosphere and the corona. Further still, diffusion perpendicular to the magnetic field can be orders of magnitude larger than that parallel to the field, in this region (e.g. Arber et al. 2007). Therefore, more work needs to be done to investigate both the global and local dissipative effects on helicity and winding.

Another important aspect for further study, revealed by our results, is the sensitivity of helicity and winding to changes at the photosphere. The present model should be updated to include a convectively unstable solar interior.