2.1 Relative magnetic helicity

Importantly, we do not require $\unicode[STIX]{x2202}V$ to be a magnetic surface, but permit non-zero $B_{n}:=\boldsymbol{B}\boldsymbol{\cdot }\hat{\boldsymbol{n}}$ on $\unicode[STIX]{x2202}V$ . This is the typical situation in astrophysical applications. As is well known, the magnetic helicity,

(2.1) $$\begin{eqnarray}h(V)=\int _{V}\boldsymbol{A}\boldsymbol{\cdot }\boldsymbol{B}\,\text{d}^{3}x,\end{eqnarray}$$

is not gauge invariant in this situation, since under a gauge transformation from $\boldsymbol{A}$ to $\boldsymbol{A}^{\prime }=\boldsymbol{A}+\unicode[STIX]{x1D735}\unicode[STIX]{x1D712}$ , the helicity becomes

(2.2) $$\begin{eqnarray}h^{\prime }(V)=\int _{V}\boldsymbol{A}\boldsymbol{\cdot }\boldsymbol{B}\,\text{d}^{3}x+\oint _{\unicode[STIX]{x2202}V}\unicode[STIX]{x1D712}B_{n}\,\text{d}^{2}x,\end{eqnarray}$$

so that its value depends on the gauge function $\unicode[STIX]{x1D712}$ . A gauge-invariant ‘relative helicity’ may be defined as

(2.3) $$\begin{eqnarray}H_{R}=\int _{V}(\boldsymbol{A}+\boldsymbol{A}_{p})\boldsymbol{\cdot }(\boldsymbol{B}-\boldsymbol{B}_{p})\,\text{d}^{3}x,\end{eqnarray}$$

where $\boldsymbol{B}_{p}$ is some arbitrary reference field whose normal component $B_{pn}$ matches $B_{n}$ on $\unicode[STIX]{x2202}V$ (Berger & Field Reference Berger and Field1984; Finn & Antonsen Jr Reference Finn and Antonsen1985). Typically $\boldsymbol{B}_{p}$ is chosen to be the unique such field that is current free (potential) within  $V$ , although this is not essential. In (2.3), $\boldsymbol{A}$ is an arbitrary vector potential for $\boldsymbol{B}$ (so that $\boldsymbol{B}=\unicode[STIX]{x1D735}\times \boldsymbol{A}$ ) and $\boldsymbol{A}_{p}$ is an arbitrary vector potential for $\boldsymbol{B}_{p}$ (so that $\boldsymbol{B}_{p}=\unicode[STIX]{x1D735}\times \boldsymbol{A}_{p}$ ). It is easy to show from (2.3) that $H_{R}$ is invariant under gauge transformation of either $\boldsymbol{A}$ or  $\boldsymbol{A}_{p}$ , for a given  $\boldsymbol{B}_{p}$ .

2.2 Field-line helicity

The concept of field-line helicity arises if we consider the magnetic helicity in a subvolume $D_{i}\in V$ . This clearly makes sense for any material (comoving) subvolume $D_{i}$ provided that $B_{n}=0$ on $\unicode[STIX]{x2202}D_{i}$ , for then

(2.4) $$\begin{eqnarray}h(D_{i})=\int _{D_{i}}\boldsymbol{A}\boldsymbol{\cdot }\boldsymbol{B}\,\text{d}^{3}x\end{eqnarray}$$

would be both gauge invariant and conserved under ideal-MHD evolution. One can envisage sub-dividing $V$ into a finer and finer union of material subvolumes, in order to gain more and more detailed information about the topological structure. As mentioned in § 1, the finest possible decomposition is where each $D_{i}$ is a thin magnetic flux tube surrounding a single magnetic field line, because magnetic field lines are material lines in an ideal evolution.

In this paper, however, we will consider only field lines that are rooted at both ends in the boundary $\unicode[STIX]{x2202}V$ , so that their surrounding tubular domains have $B_{n}\neq 0$ on the ends, violating this boundary condition. Figure 1 shows such a field line, $L(\boldsymbol{x})$ , along with a surrounding flux tube $V_{\unicode[STIX]{x1D716}}(\boldsymbol{x})$ of radius $\unicode[STIX]{x1D700}$ (in some cross-section). We label $L(\boldsymbol{x})$ by one of its endpoints $\boldsymbol{x}\in \unicode[STIX]{x2202}V$ . Throughout this paper, $\boldsymbol{x}$ will exclusively denote a point on the boundary $\unicode[STIX]{x2202}V$ .

Because $B_{n}\neq 0$ on $\unicode[STIX]{x2202}V_{\unicode[STIX]{x1D716}}(\boldsymbol{x})$ , we expect $h(V_{\unicode[STIX]{x1D700}}(\boldsymbol{x}))$ to change due to either motion of the field line endpoints or gauge transformation of  $\boldsymbol{A}$ . The latter is a more significant problem, and will be addressed in detail in this paper. First, however, we take the limit to an infinitesimally thin tube. Since $h(V_{\unicode[STIX]{x1D700}}(\boldsymbol{x}))$ is a volume integral, we normalise by the flux $\unicode[STIX]{x1D6F7}_{\unicode[STIX]{x1D700}}(\boldsymbol{x})$ of the tube and define

(2.5) $$\begin{eqnarray}{\mathcal{A}}(\boldsymbol{x})=\lim _{\unicode[STIX]{x1D700}\rightarrow 0}\frac{h(V_{\unicode[STIX]{x1D700}}(\boldsymbol{x}))}{\unicode[STIX]{x1D6F7}_{\unicode[STIX]{x1D700}}(\boldsymbol{x})}.\end{eqnarray}$$

In the limit $\unicode[STIX]{x1D700}\rightarrow 0$ , the tube collapses to the line $L(\boldsymbol{x})$ , and ${\mathcal{A}}(\boldsymbol{x})$ tends to a well-defined limit – independent of the choice of cross-section – which Berger (Reference Berger1988) calls the field-line helicity. Note that, although ${\mathcal{A}}(\boldsymbol{x})$ is defined by integrating only over the tubular volume $V_{\unicode[STIX]{x1D700}}(\boldsymbol{x})$ , it nevertheless contains information about how this particular flux tube interacts with the magnetic field outside, due to the fact that $\boldsymbol{A}$ is defined globally. In this paper, we will assume that $\boldsymbol{B}$ has no ergodic field lines with infinite length, as we have not shown that a well-defined limit exists in that case.

Figure 1.

Definition sketch for (2.5), showing a thin flux tube $V_{\unicode[STIX]{x1D700}}(\boldsymbol{x})$ of radius $\unicode[STIX]{x1D700}$ and magnetic flux $\unicode[STIX]{x1D6F7}_{\unicode[STIX]{x1D700}}(\boldsymbol{x})$ around the field line $L(\boldsymbol{x})$ that is rooted at $\boldsymbol{x}\in \unicode[STIX]{x2202}V$ .

The definition (2.5) shows that ${\mathcal{A}}(\boldsymbol{x})$ has the dimensions of a magnetic flux. A more convenient formula may be obtained by writing the volume integral $h(V_{\unicode[STIX]{x1D700}}(\boldsymbol{x}))$ as an integral along $L(\boldsymbol{x})$ of cross-sectional surface integrals. In the limit $V_{\unicode[STIX]{x1D700}}(\boldsymbol{x})\rightarrow L(\boldsymbol{x})$ , the vector potential $\boldsymbol{A}$ is constant on each cross-section, so that

(2.6) $$\begin{eqnarray}{\mathcal{A}}(\boldsymbol{x})=\lim _{\unicode[STIX]{x1D700}\rightarrow 0}\frac{\displaystyle \int _{V_{\unicode[STIX]{x1D700}}(\boldsymbol{x})}\unicode[STIX]{x1D6F7}_{\unicode[STIX]{x1D700}}(\boldsymbol{x})\boldsymbol{A}\boldsymbol{\cdot }\text{d}\boldsymbol{l}}{\unicode[STIX]{x1D6F7}_{\unicode[STIX]{x1D700}}(\boldsymbol{x})}=\int _{L(\boldsymbol{x})}\boldsymbol{A}\boldsymbol{\cdot }\text{d}\boldsymbol{l}.\end{eqnarray}$$

This formula shows clearly that, were $L(\boldsymbol{x})$ to be a contractible closed loop, then ${\mathcal{A}}(\boldsymbol{x})$ would be gauge independent, and by Stokes’ theorem it would simply be the magnetic flux linked through that loop. Clearly ${\mathcal{A}}(\boldsymbol{x})$ would then be an ideal-MHD invariant, representing the magnetic flux linked with $L(\boldsymbol{x})$ . In our case, $L(\boldsymbol{x})$ is not a closed loop. Nevertheless, ${\mathcal{A}}(\boldsymbol{x})$ remains an ideal-MHD invariant provided that there are no boundary motions, and also that the gauge of $\boldsymbol{A}$ on the boundary is fixed. Indeed a gauge transformation from $\boldsymbol{A}$ to $\boldsymbol{A}^{\prime }=\boldsymbol{A}+\unicode[STIX]{x1D735}\unicode[STIX]{x1D712}$ will change the value of ${\mathcal{A}}(\boldsymbol{x})$ to

(2.7) $$\begin{eqnarray}{\mathcal{A}}^{\prime }(\boldsymbol{x})={\mathcal{A}}(\boldsymbol{x})+\unicode[STIX]{x1D712}(\boldsymbol{x}_{+})-\unicode[STIX]{x1D712}(\boldsymbol{x}),\end{eqnarray}$$

where $\boldsymbol{x}_{+}$ is the other end of the field line, as in figure 1. This is simply the limiting version of the formula (2.2). In fact, since $\unicode[STIX]{x2202}V$ is a single closed surface, it follows that ${\mathcal{A}}(\boldsymbol{x})$ , like $h(V)$ , depends only on $\hat{\boldsymbol{n}}\times \boldsymbol{A}$ . For if $\hat{\boldsymbol{n}}\times \boldsymbol{A}^{\prime }=\hat{\boldsymbol{n}}\times \boldsymbol{A}$ on $\unicode[STIX]{x2202}V$ , then $\hat{\boldsymbol{n}}\times \unicode[STIX]{x1D735}\unicode[STIX]{x1D712}=0$ on $\unicode[STIX]{x2202}V$ , so that $\unicode[STIX]{x1D712}$ is constant over all of $\unicode[STIX]{x2202}V$ and ${\mathcal{A}}^{\prime }(\boldsymbol{x})={\mathcal{A}}(\boldsymbol{x})$ . This would not be true for a domain where $\unicode[STIX]{x2202}V$ is not a single connected surface, for example a spherical shell. As we will see in this paper, the art of working with field-line helicity is to choose an informative gauge for $\hat{\boldsymbol{n}}\times \boldsymbol{A}$ .

Finally, we note the relation between ${\mathcal{A}}(\boldsymbol{x})$ and the overall helicity $h(V)$ . When the field lines $L(\boldsymbol{x})$ partition the whole volume (i.e. there are no closed or ergodic field lines), integrating (2.5) over all field lines, weighted by magnetic flux, will give $h(V)$ (Berger Reference Berger1988). In other words,

(2.8) $$\begin{eqnarray}\frac{1}{2}\oint _{\unicode[STIX]{x2202}V}{\mathcal{A}}|B_{n}|\,\text{d}^{2}x=h(V),\end{eqnarray}$$

where the factor half arises because each field line has two endpoints on $\unicode[STIX]{x2202}V$ .

3.1 Towards a definition

One option is to use ${\mathcal{A}}(\boldsymbol{x})$ directly, but restrict the gauge  $\boldsymbol{A}$ . We know that, whatever the gauge, ${\mathcal{A}}(\boldsymbol{x})$ is invariant under an ideal-MHD evolution for every field line provided that $\hat{\boldsymbol{n}}\times \boldsymbol{A}$ on $\unicode[STIX]{x2202}V$ remains fixed in time. But in an arbitrary gauge, condition (3.1) will not be satisfied, since $h(V)\neq H_{R}$ in general. Nevertheless, there are a restricted family of gauges of $\boldsymbol{A}$ where $h(V)=H_{R}$ . To see this, write (2.3) in the form

(3.2) $$\begin{eqnarray}H_{R}=\int _{V}\boldsymbol{A}\boldsymbol{\cdot }\boldsymbol{B}\,\text{d}^{3}x+\oint _{\unicode[STIX]{x2202}V}\boldsymbol{A}\times \boldsymbol{A}_{p}\boldsymbol{\cdot }\hat{\boldsymbol{n}}\,\text{d}^{2}x-\int _{V}\boldsymbol{A}_{p}\boldsymbol{\cdot }\boldsymbol{B}_{p}\,\text{d}^{3}x.\end{eqnarray}$$

It is well known that the boundary term vanishes if we restrict $\hat{\boldsymbol{n}}\times \boldsymbol{A}=\hat{\boldsymbol{n}}\times \boldsymbol{A}_{p}$ on $\unicode[STIX]{x2202}V$ (Barnes Reference Barnes1988; Berger Reference Berger1988). This is always possible since $B_{pn}=B_{n}$ on $\unicode[STIX]{x2202}V$ (§ 4.3 shows how to construct such an  $\boldsymbol{A}$ , given an arbitrary  $\boldsymbol{A}_{p}$ ). We then need only choose $\boldsymbol{A}_{p}$ such that $\int _{V}\boldsymbol{A}_{p}\boldsymbol{\cdot }\boldsymbol{B}_{p}\,\text{d}^{3}x=0$ , and we will have $H_{R}=h(V)$ and hence (3.1). In general, this is not possible, since $\boldsymbol{B}_{p}$ may have a closed-loop field line $L$ within $V$ whose gauge invariant $\oint _{L}\boldsymbol{A}_{p}\boldsymbol{\cdot }\text{d}\boldsymbol{l}$ does not vanish. However, in the case where $\boldsymbol{B}_{p}$ is a potential field, it cannot have any closed-loop field lines, and we can always find a gauge for $\boldsymbol{A}_{p}$ such that $\int _{V}\boldsymbol{A}_{p}\boldsymbol{\cdot }\boldsymbol{B}_{p}\,\text{d}^{3}x=0$ (e.g. Barnes Reference Barnes1988; Berger Reference Berger1988).

Even when $\boldsymbol{B}_{p}$ is a potential field, the conditions that $\int _{V}\boldsymbol{A}_{p}\boldsymbol{\cdot }\boldsymbol{B}_{p}\,\text{d}^{3}x=0$ and $\hat{\boldsymbol{n}}\times \boldsymbol{A}=\hat{\boldsymbol{n}}\times \boldsymbol{A}_{p}$ on $\unicode[STIX]{x2202}V$ do not uniquely define ${\mathcal{A}}(\boldsymbol{x})$ , because there is still freedom in the choice of  $\boldsymbol{A}_{p}$ . For example, if we have such a gauge $\boldsymbol{A}_{p}$ and change to $\boldsymbol{A}_{p}^{\prime }=\boldsymbol{A}_{p}+\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}$ , then

(3.3) $$\begin{eqnarray}\int _{V}\boldsymbol{A}_{p}^{\prime }\boldsymbol{\cdot }\boldsymbol{B}_{p}\,\text{d}^{3}x=\oint _{\unicode[STIX]{x2202}V}\unicode[STIX]{x1D719}B_{n}\,\text{d}^{2}x.\end{eqnarray}$$

Provided $B_{n}\neq 0$ , this integral may still vanish even if $\unicode[STIX]{x1D719}\neq 0$ on $\unicode[STIX]{x2202}V$ .

3.2 General definition

In practice, it may be inconvenient or (when $\boldsymbol{B}_{p}$ is non-potential) impossible to find a gauge $\boldsymbol{A}_{p}$ such that $\int _{V}\boldsymbol{A}_{p}\boldsymbol{\cdot }\boldsymbol{B}_{p}\,\text{d}^{3}x=0$ . Therefore, we propose instead an alternative definition for relative field-line helicity that does not impose this requirement on  $\boldsymbol{A}_{p}$ . Namely,

(3.4) $$\begin{eqnarray}{\mathcal{A}}_{R}(\boldsymbol{x}):={\mathcal{A}}(\boldsymbol{x})-{\mathcal{A}}_{p}(\boldsymbol{x}),\end{eqnarray}$$

where ${\mathcal{A}}(\boldsymbol{x})$ is the usual field-line helicity (2.6) and

(3.5) $$\begin{eqnarray}{\mathcal{A}}_{p}(\boldsymbol{x})=\int _{L_{p}(\boldsymbol{x})}\boldsymbol{A}_{p}\boldsymbol{\cdot }\text{d}\boldsymbol{l},\end{eqnarray}$$

which is the field-line helicity of the reference field $\boldsymbol{B}_{p}$ on its own field line rooted at the same point $\boldsymbol{x}\in \unicode[STIX]{x2202}V$ . Importantly, our definition assumes the gauge condition on $\boldsymbol{A}$ that

(3.6) $$\begin{eqnarray}\hat{\boldsymbol{n}}\times \boldsymbol{A}=\hat{\boldsymbol{n}}\times \boldsymbol{A}_{p}\quad \text{on}~\unicode[STIX]{x2202}V.\end{eqnarray}$$

It follows that ${\mathcal{A}}_{R}(\boldsymbol{x})$ will satisfy (3.1), whatever the gauge of  $\boldsymbol{A}_{p}$ . Note that, unlike ${\mathcal{A}}(\boldsymbol{x})$ , this ${\mathcal{A}}_{R}(\boldsymbol{x})$ will not, in general, have the same value at both ends of a given field line $L(\boldsymbol{x})$ , since the corresponding reference field line $L_{p}(\boldsymbol{x})$ at each end will be different.

A nice property of ${\mathcal{A}}_{R}(\boldsymbol{x})$ is that ${\mathcal{A}}_{R}(\boldsymbol{x})=0$ whenever $\boldsymbol{B}_{p}=\boldsymbol{B}$ all along $L(\boldsymbol{x})$ for some $\boldsymbol{x}\in \unicode[STIX]{x2202}V$ . For in that case, $L(\boldsymbol{x})=L_{p}(\boldsymbol{x})$ and we must have $\boldsymbol{A}=\boldsymbol{A}_{p}+\unicode[STIX]{x1D735}\unicode[STIX]{x1D712}$ along this line. Then

(3.7) $$\begin{eqnarray}{\mathcal{A}}_{R}(\boldsymbol{x})=\int _{L(\boldsymbol{x})}(\boldsymbol{A}-\boldsymbol{A}_{p})\boldsymbol{\cdot }\text{d}\boldsymbol{l}=\int _{L(\boldsymbol{x})}\unicode[STIX]{x1D735}\unicode[STIX]{x1D712}\boldsymbol{\cdot }\text{d}\boldsymbol{l}.\end{eqnarray}$$

However, since $\unicode[STIX]{x2202}V$ is a connected surface, and $\hat{\boldsymbol{n}}\times \unicode[STIX]{x1D735}\unicode[STIX]{x1D712}=0$ on this surface thanks to our gauge restriction, it follows that $\unicode[STIX]{x1D712}$ is constant on $\unicode[STIX]{x2202}V$ and hence ${\mathcal{A}}_{R}(\boldsymbol{x})=0$ . This property makes ${\mathcal{A}}_{R}$ useful for identifying non-potential regions within  $\boldsymbol{B}$ , and will not generally be shared by the basic field-line helicity ${\mathcal{A}}$ , even in a gauge where $\int _{V}\boldsymbol{A}_{p}\boldsymbol{\cdot }\boldsymbol{B}_{p}\,\text{d}^{3}x=0$ overall.

3.3 Gauge dependence and physical meaning

Unfortunately, our definition of ${\mathcal{A}}_{R}$ in (3.4), even with the restriction that $\hat{\boldsymbol{n}}\times \boldsymbol{A}=\hat{\boldsymbol{n}}\times \boldsymbol{A}_{p}$ on $\unicode[STIX]{x2202}V$ , is still not uniquely defined, and depends on the gauge of  $\boldsymbol{A}_{p}$ . To see this, consider a gauge transformation from $\boldsymbol{A}_{p}$ to $\boldsymbol{A}_{p}^{\prime }=\boldsymbol{A}_{p}+\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}$ . To preserve the boundary restriction, we must also change the gauge of $\boldsymbol{A}$ to $\boldsymbol{A}^{\prime }=\boldsymbol{A}+\unicode[STIX]{x1D735}\unicode[STIX]{x1D712}$ , where $\hat{\boldsymbol{n}}\times \unicode[STIX]{x1D735}(\unicode[STIX]{x1D719}-\unicode[STIX]{x1D712})=0$ on $\unicode[STIX]{x2202}V$ . It follows that the two gauge functions can differ on $\unicode[STIX]{x2202}V$ only by a global constant (since $\unicode[STIX]{x2202}V$ is connected). Using this fact, we find that ${\mathcal{A}}_{R}(\boldsymbol{x})$ goes to

(3.8) $$\begin{eqnarray}\displaystyle {{\mathcal{A}}_{R}}^{\prime }(\boldsymbol{x}) & = & \displaystyle {\mathcal{A}}_{R}(\boldsymbol{x})+\int _{L(\boldsymbol{x})}\unicode[STIX]{x1D735}\unicode[STIX]{x1D712}\boldsymbol{\cdot }\text{d}\boldsymbol{l}-\int _{L_{p}(\boldsymbol{x})}\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}\boldsymbol{\cdot }\text{d}\boldsymbol{l},\end{eqnarray}$$

(3.9) $$\begin{eqnarray}\displaystyle & = & \displaystyle {\mathcal{A}}_{R}(\boldsymbol{x})+\unicode[STIX]{x1D719}(\boldsymbol{x}_{+})-\unicode[STIX]{x1D719}(\boldsymbol{x}_{p+}).\end{eqnarray}$$

When $L(\boldsymbol{x})\neq L_{p}(\boldsymbol{x})$ , the other ends of the field lines (denoted $\boldsymbol{x}_{+}$ and $\boldsymbol{x}_{p+}$ ) will differ, so the value of ${\mathcal{A}}_{R}(\boldsymbol{x})$ will change.

To see the physical meaning of this gauge dependence, we can interpret ${\mathcal{A}}_{R}(\boldsymbol{x})$ as a magnetic flux (cf. Yeates & Hornig Reference Yeates and Hornig2016). This is illustrated in figure 2, where $L(\boldsymbol{x})$ is a field line of $\boldsymbol{B}$ and $L_{p}(\boldsymbol{x})$ is a field line of $\boldsymbol{B}_{p}$ rooted at the same point  $\boldsymbol{x}$ . The curve $\unicode[STIX]{x1D6FE}\subset \unicode[STIX]{x2202}V$ in figure 2(a) closes the loop and defines a surface in $V$ bounded by $L(\boldsymbol{x})$ , $\unicode[STIX]{x1D6FE}$ and $-L_{p}(\boldsymbol{x})$ . If $\unicode[STIX]{x1D6FE}$ is chosen so that $\int _{\unicode[STIX]{x1D6FE}}\boldsymbol{A}_{p}\boldsymbol{\cdot }\text{d}\boldsymbol{l}=0$ , then ${\mathcal{A}}_{R}(\boldsymbol{x})$ will be precisely the magnetic flux through this surface, by Stokes’ theorem. Provided the field-line endpoints and $\hat{\boldsymbol{n}}\times \boldsymbol{A}_{p}$ on $\unicode[STIX]{x2202}V$ and remain fixed in time, this flux will be an ideal-MHD invariant. In any gauge, it is possible to find such a curve $\unicode[STIX]{x1D6FE}$ along which the line integral vanishes, as proven by Yeates & Hornig (Reference Yeates and Hornig2016). In fact, there are an infinite number of such curves (whose corresponding surfaces all have the same flux). The possible $\unicode[STIX]{x1D6FE}$ depend on the gauge of $\hat{\boldsymbol{n}}\times \boldsymbol{A}_{p}$ , as sketched in figure 2(b) which shows a different curve $\unicode[STIX]{x1D6FE}^{\prime }$ arising in a different gauge  $\boldsymbol{A}^{\prime }$ . This loop encloses a different amount of magnetic flux, corresponding to the different value of ${{\mathcal{A}}_{R}}^{\prime }(\boldsymbol{x})$ in this new gauge. In summary, the gauge dependence of ${\mathcal{A}}_{R}$ may be viewed simply as a choice of how to close flux surfaces on the boundary. Whichever choice is made, ${\mathcal{A}}_{R}(\boldsymbol{x})$ is an ideal-MHD invariant.

Figure 2.

Interpretation of ${\mathcal{A}}_{R}(\boldsymbol{x})$ as a magnetic flux, for two different gauges in (a) and (b). The field lines $L(\boldsymbol{x})$ and $L_{p}(\boldsymbol{x})$ are the same in each case, but the possible curves $\unicode[STIX]{x1D6FE}$ and $\unicode[STIX]{x1D6FE}^{\prime }$ that complete the loop by linking $\boldsymbol{x}_{+}$ and $\boldsymbol{x}_{p+}$ on $\unicode[STIX]{x2202}V$ differ between the two gauges, so that ${\mathcal{A}}_{R}(\boldsymbol{x})\neq {{\mathcal{A}}_{R}}^{\prime }(\boldsymbol{x})$ .

At first sight, it may seem from (3.9) that we can choose any distribution of ${\mathcal{A}}_{R}$ that we like on field lines that intersect $\unicode[STIX]{x2202}V$ , simply by changing the gauge  $\unicode[STIX]{x1D719}$ . Indeed, this is true. However, choosing $\unicode[STIX]{x1D719}$ in this way requires you to first know ${\mathcal{A}}_{R}$ for $\boldsymbol{B}$ in some initial gauge. If we choose the gauge a priori, based only on  $\boldsymbol{B}_{p}$ , then the resulting ${\mathcal{A}}_{R}$ can give us meaningful information about  $\boldsymbol{B}$ .

In § 5, we will illustrate the effect on ${\mathcal{A}}_{R}$ of some different choices of $\boldsymbol{A}_{p}$ in particular examples. However, we will first propose a general candidate for the ‘best’ choice of gauge, that can be used to uniquely define  ${\mathcal{A}}_{R}$ .

3.4 Minimal gauge

Since ${\mathcal{A}}_{R}$ depends only on the distribution of $\hat{\boldsymbol{n}}\times \boldsymbol{A}_{p}$ on $\unicode[STIX]{x2202}V$ , our idea is to choose this gauge by minimising the integral $\oint _{\unicode[STIX]{x2202}V}|\hat{\boldsymbol{n}}\times \boldsymbol{A}_{p}|^{2}\,\text{d}^{2}x$ over the whole boundary. Intuitively, this will give us the smallest overall boundary contribution to  ${\mathcal{A}}_{R}$ , for a given reference field $\boldsymbol{B}_{p}$ (whether potential or not). In appendix A, we prove that this integral is minimised if

(3.10) $$\begin{eqnarray}\unicode[STIX]{x1D735}_{h}\boldsymbol{\cdot }\boldsymbol{A}_{p}=0\quad \text{on}~\unicode[STIX]{x2202}V,\end{eqnarray}$$

where $\unicode[STIX]{x1D735}_{h}\boldsymbol{\cdot }$ denotes the two-dimensional divergence in the plane of $\unicode[STIX]{x2202}V$ . This is the so-called ‘universal gauge’ condition of Hornig (Reference Hornig2006). Although this condition does not specify $\boldsymbol{A}_{p}$ uniquely within  $V$ , it does specify $\hat{\boldsymbol{n}}\times \boldsymbol{A}_{p}$ uniquely on $\unicode[STIX]{x2202}V$ , and hence uniquely specifies  ${\mathcal{A}}_{R}$ . To see this, suppose both $\boldsymbol{A}_{p}$ and $\boldsymbol{A}_{p}^{\prime }=\boldsymbol{A}_{p}+\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}$ satisfy (3.10). Then we must have $\unicode[STIX]{x1D6FB}_{h}^{2}\unicode[STIX]{x1D719}=0$ , but since $\unicode[STIX]{x2202}V$ is a closed surface, this implies that $\unicode[STIX]{x1D719}$ is constant on $\unicode[STIX]{x2202}V$ , so that $\hat{\boldsymbol{n}}\times \boldsymbol{A}_{p}^{\prime }=\hat{\boldsymbol{n}}\times \boldsymbol{A}_{p}$ . Note that ${\mathcal{A}}_{p}$ does not necessarily vanish in this gauge when $\boldsymbol{B}_{p}$ is a potential field – we will see this in § 5.

In the rest of this paper, we will denote vector potentials satisfying (3.10) by $\boldsymbol{A}^{\ast }$ or  $\boldsymbol{A}_{p}^{\ast }$ , the corresponding field-line helicities by ${\mathcal{A}}^{\ast }$ and ${\mathcal{A}}_{p}^{\ast }$ and the corresponding relative field-line helicity by  ${\mathcal{A}}_{R}^{\ast }$ . In § 5, we will compute ${\mathcal{A}}_{R}^{\ast }$ for some example magnetic fields, and compare it to ${\mathcal{A}}_{R}$ in some other gauges.

3.5 Summary of the proposed definition

To summarise, let $\boldsymbol{A}_{p}^{\ast }$ be a vector potential for $\boldsymbol{B}_{p}$ satisfying

(3.11) $$\begin{eqnarray}\unicode[STIX]{x1D735}_{h}\boldsymbol{\cdot }\boldsymbol{A}_{p}^{\ast }=0\quad \text{on}~\unicode[STIX]{x2202}V,\end{eqnarray}$$

and let $\boldsymbol{A}^{\ast }$ be a vector potential for $\boldsymbol{B}$ satisfying

(3.12) $$\begin{eqnarray}\hat{\boldsymbol{n}}\times \boldsymbol{A}^{\ast }=\hat{\boldsymbol{n}}\times \boldsymbol{A}_{p}^{\ast }\quad \text{on}~\unicode[STIX]{x2202}V.\end{eqnarray}$$

Then our proposed definition of relative field-line helicity at $\boldsymbol{x}\in \unicode[STIX]{x2202}V$ is the difference

(3.13) $$\begin{eqnarray}{\mathcal{A}}_{R}^{\ast }(\boldsymbol{x}):={\mathcal{A}}^{\ast }(\boldsymbol{x})-{\mathcal{A}}_{p}^{\ast }(\boldsymbol{x}),\end{eqnarray}$$

where ${\mathcal{A}}^{\ast }$ , ${\mathcal{A}}_{p}^{\ast }$ are the field-line helicities of $\boldsymbol{B}$ and $\boldsymbol{B}_{p}$ in these two gauges, defined along their respective magnetic field lines.

4.1 Reference potential field

In a Cartesian box, the potential reference field $\boldsymbol{B}_{p}$ that matches $B_{n}$ on all six boundary faces is straightforward to compute. We write $\boldsymbol{B}_{p}=\unicode[STIX]{x1D735}u$ leading to the three-dimensional Laplace equation $\unicode[STIX]{x1D6FB}^{2}u=0$ . On the staggered mesh, the potential $u$ is located at the centre of each three-dimensional cell, so that each component of $\boldsymbol{B}_{p}$ may be computed from $u$ by a central difference. We approximate the Laplacian operator $\unicode[STIX]{x1D6FB}^{2}$ with a second-order (7-point) finite-difference stencil, for which the solution is efficiently and simply obtained using a fast-Poisson technique with Neumann boundary conditions (Press et al. Reference Press, Teukolsky, Vetterling and Flannery1992). With this technique, the current density $\unicode[STIX]{x1D735}\times \boldsymbol{B}_{p}$ vanishes to machine precision when computed from $\boldsymbol{B}_{p}$ with central differences.

4.2 Minimal gauge vector potential

To find $\boldsymbol{A}^{\ast }$ (or $\boldsymbol{A}_{p}^{\ast }$ ), we have developed a numerical routine that starts from $\boldsymbol{A}$ (or  $\boldsymbol{A}_{p}$ ) in an arbitrary gauge, and performs the necessary gauge transformation. In particular, for the case of  $\boldsymbol{A}$ , we seek a scalar function $\unicode[STIX]{x1D712}(x,y,z)$ on $V$ such that $\boldsymbol{A}^{\ast }=\boldsymbol{A}+\unicode[STIX]{x1D735}\unicode[STIX]{x1D712}$ . The condition $\unicode[STIX]{x1D735}_{h}\boldsymbol{\cdot }\boldsymbol{A}^{\ast }=0$ on $\unicode[STIX]{x2202}V$ requires that $\overline{\unicode[STIX]{x1D712}}:=\unicode[STIX]{x1D712}|_{\unicode[STIX]{x2202}V}$ satisfy the two-dimensional Poisson equation

(4.1) $$\begin{eqnarray}\unicode[STIX]{x1D6FB}_{h}^{2}\overline{\unicode[STIX]{x1D712}}=-\unicode[STIX]{x1D735}_{h}\boldsymbol{\cdot }\boldsymbol{A}\end{eqnarray}$$

over this whole boundary $\unicode[STIX]{x2202}V$ . Since $\unicode[STIX]{x2202}V$ is a closed surface, the right-hand side necessarily satisfies the compatibility condition $\oint _{\unicode[STIX]{x2202}V}\unicode[STIX]{x1D735}_{h}\boldsymbol{\cdot }\boldsymbol{A}\,\text{d}^{2}x=0$ , so that the solution for $\overline{\unicode[STIX]{x1D712}}$ exists and is unique up to an additive constant  $\overline{\unicode[STIX]{x1D712}}_{0}$ .

Equation (4.1) is solved numerically by a finite-difference method, approximating the Laplacian operator with the standard 5-point stencil. For our Cartesian domain, care must be taken to couple values of $\overline{\unicode[STIX]{x1D712}}$ on neighbouring faces so that the normal component of $\unicode[STIX]{x1D735}_{h}\overline{\unicode[STIX]{x1D712}}$ is continuous. (Strictly speaking, we are approximating a weak solution to (4.1), since the differential operators in (4.1) are not defined on the edges of the cube.) This coupling between faces precludes the use of a fast-Poisson solver, but a direct solution of the resulting (sparse) linear system remains practical for this two-dimensional problem. Although $\boldsymbol{A}^{\ast }$ does not depend on the additive constant  $\overline{\unicode[STIX]{x1D712}}_{0}$ , we ensure that the linear system is invertible by fixing $\overline{\unicode[STIX]{x1D712}}=0$ at one vertex of $\unicode[STIX]{x2202}V$ .

Given the solution for $\overline{\unicode[STIX]{x1D712}}$ , we can then choose any arbitrary extension of $\unicode[STIX]{x1D712}$ into the interior of  $V$ , since this choice does not affect  ${\mathcal{A}}^{\ast }$ . It is simplest to take $\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D712}=0$ in the volume, and solve this three-dimensional Laplace equation with Dirichlet boundary conditions $\unicode[STIX]{x1D712}=\overline{\unicode[STIX]{x1D712}}$ on $\unicode[STIX]{x2202}V$ . This can be solved using a standard fast-Poisson method.

4.3 Matching the reference gauge

Although not required for computing ${\mathcal{A}}_{R}^{\ast }$ , we have also implemented a numerical routine for changing the gauge of $\boldsymbol{A}$ to $\boldsymbol{A}^{\prime }$ so that $\hat{\boldsymbol{n}}\times \boldsymbol{A}^{\prime }=\hat{\boldsymbol{n}}\times \boldsymbol{A}_{p}$ on $\unicode[STIX]{x2202}V$ , given some  $\boldsymbol{A}_{p}$ . This will be used in § 5 to illustrate the computation of ${\mathcal{A}}_{R}$ in gauges other than the minimal gauge.

Writing $\boldsymbol{A}^{\prime }=\boldsymbol{A}+\unicode[STIX]{x1D735}\unicode[STIX]{x1D712}$ , we need to compute $\unicode[STIX]{x1D712}(x,y,z)\in V$ such that $\hat{\boldsymbol{n}}\times \boldsymbol{A}^{\prime }=\hat{\boldsymbol{n}}\times \boldsymbol{A}_{p}$ on $\unicode[STIX]{x2202}V$ . This may be done by first solving for $\unicode[STIX]{x1D712}|_{\unicode[STIX]{x2202}V}=\overline{\unicode[STIX]{x1D712}}$ on $\unicode[STIX]{x2202}V$ , then extending to a solution on the interior of $V$ . We need $\overline{\unicode[STIX]{x1D712}}$ to satisfy $\unicode[STIX]{x1D735}_{h}\overline{\unicode[STIX]{x1D712}}=\boldsymbol{A}_{ph}-\boldsymbol{A}_{h}$ , so take the divergence and solve the resulting two-dimensional Poisson equation

(4.2) $$\begin{eqnarray}\unicode[STIX]{x1D6FB}_{h}^{2}\overline{\unicode[STIX]{x1D712}}=\unicode[STIX]{x1D735}_{h}\boldsymbol{\cdot }(\boldsymbol{A}_{p}-\boldsymbol{A})\end{eqnarray}$$

on $\unicode[STIX]{x2202}V$ . As in § 4.2, a solution exists because the compatibility condition is satisfied. Since the Neumann boundary condition $\hat{\boldsymbol{n}}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\overline{\unicode[STIX]{x1D712}}=\hat{\boldsymbol{n}}\boldsymbol{\cdot }(\boldsymbol{A}_{p}-\boldsymbol{A})$ is fixed on the edges of each face, it is possible to solve (4.2) separately on each face. This means that a fast-Poisson solver can readily be used, unlike in § 4.2. This gives $\overline{\unicode[STIX]{x1D712}}$ on each face $S_{i}$ up to a constant  $\overline{\unicode[STIX]{x1D712}}_{i}$ , for $i=1,\ldots ,6$ . One of these constants may be chosen arbitrarily, and the other five are then easily determined by imposing continuity of $\overline{\unicode[STIX]{x1D712}}$ at each edge. (Although there are 12 edges, the corresponding jumps in $\overline{\unicode[STIX]{x1D712}}$ are not independent, so the problem is well posed.)

Once the continuous solution for $\overline{\unicode[STIX]{x1D712}}$ is obtained, this may be extended to $\unicode[STIX]{x1D712}$ in the interior of $V$ in the same way as § 4.2, namely solving the three-dimensional Laplace equation $\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D712}=0$ in $V$ with Dirichlet boundary conditions $\unicode[STIX]{x1D712}=\overline{\unicode[STIX]{x1D712}}$ on $\unicode[STIX]{x2202}V$ . Again, this choice of interior distribution has no effect on the field-line helicity.

5.1 Low–Lou equilibrium

First, we consider one of the class of nonlinear force-free equilibria derived by Low & Lou (Reference Low and Lou1990), which is given analytically except for the solution of a single ordinary differential equation. Our calculations start from a three-dimensional datacube of  $\boldsymbol{B}$ , computed by Valori et al. (Reference Valori, Pariat, Anfinogentov, Chen, Georgoulis, Guo, Liu, Moraitis, Thalmann and Yang2016) and kindly shared by these authors. This corresponds to the specific solution with parameters $n=1$ , $l=0.3$ , $\unicode[STIX]{x1D719}=\unicode[STIX]{x03C0}/4$ , in a Cartesian box $V=[-1,1]\times [-1,1]\times [0,2]$ . We used four datacubes with uniform mesh spacing  $\unicode[STIX]{x0394}x=2/32$ , $2/64$ , $2/128$ and $2/256$ . As described in § 4, our numerical code uses a staggered mesh, so we first interpolated the original $\boldsymbol{B}$ components from cell vertices to the centres of the cell faces, using trilinear interpolation. We also tested tricubic interpolation here but found it unnecessary.

Figure 3 shows the magnetic field lines for this magnetic field $\boldsymbol{B}$ (a) and for the corresponding reference potential field $\boldsymbol{B}_{p}$ (b). The latter was computed numerically to match $B_{n}$ on all six boundaries, using the method described in § 4.1. In figure 3, each field line $L(\boldsymbol{x})$ is coloured by its integrated parallel current density, defined as

(5.1) $$\begin{eqnarray}J_{\Vert }(\boldsymbol{x})=\int _{L(\boldsymbol{x})}\unicode[STIX]{x1D735}\times \boldsymbol{B}\boldsymbol{\cdot }\,\text{d}\boldsymbol{l}\end{eqnarray}$$

and normalised by $B_{0}$ , the mean value of $|B_{z}(x,y,0)|$ . Maps of $J_{\Vert }/B_{0}$ for all field lines seeded from the $z=0$ boundary are shown in (c,d). In the Low–Lou equilibrium, the electric currents are smoothly distributed within the arcade of magnetic loops surrounding the central polarity-inversion line. For  $\boldsymbol{B}_{p}$ , we have $\unicode[STIX]{x1D735}\times \boldsymbol{B}_{p}\equiv \mathbf{0}$ and hence $J_{\Vert }$ vanishes (to machine precision).

Next, we have computed maps of the field-line helicities ${\mathcal{A}}$ , ${\mathcal{A}}_{p}$ and the relative field-line helicity ${\mathcal{A}}_{R}:={\mathcal{A}}-{\mathcal{A}}_{p}$ using four different gauges. These include the minimal gauge  $\boldsymbol{A}_{p}^{\ast }$ , but also three others, denoted $\boldsymbol{A}_{p}^{a}$ , $\boldsymbol{A}_{p}^{b}$ , $\boldsymbol{A}_{p}^{c}$ and defined as follows. All three satisfy the DeVore (Reference DeVore2000) condition that $A_{z}=0$ . The first is integrated upward in from the lower $z$ boundary,

(5.2) $$\begin{eqnarray}\boldsymbol{A}_{p}^{a}(x,y,z)=\boldsymbol{A}_{0}^{a}(x,y)-\hat{\boldsymbol{z}}\times \int _{0}^{z}\boldsymbol{B}(x,y,z^{\prime })\,\text{d}z^{\prime },\end{eqnarray}$$

and uses the Coulomb gauge form $\boldsymbol{A}_{0}^{a}=\unicode[STIX]{x1D735}\times (\unicode[STIX]{x1D6F9}\hat{\boldsymbol{z}})$ for the boundary term, with $\unicode[STIX]{x1D6F9}$ computed by solving the Poisson equation $\unicode[STIX]{x1D6FB}_{h}^{2}\unicode[STIX]{x1D6F9}=-B_{z}(x,y,0)$ (we separate the monopole component and solve for the remainder using periodic boundary conditions). The second gauge is similar but integrated down from the upper $z$ boundary, so that

(5.3) $$\begin{eqnarray}\boldsymbol{A}_{p}^{b}(x,y,z)=\boldsymbol{A}_{0}^{b}(x,y)+\hat{\boldsymbol{z}}\times \int _{z}^{2}\boldsymbol{B}(x,y,z^{\prime })\,\text{d}z^{\prime }\end{eqnarray}$$

and $\boldsymbol{A}_{0}^{b}$ is now computed on $z=2$ from $B_{z}(x,y,2)$ . The third is an alternative DeVore gauge that still has the form

(5.4) $$\begin{eqnarray}\boldsymbol{A}_{p}^{c}=\boldsymbol{A}_{0}^{c}(x,y)-\hat{\boldsymbol{z}}\times \int _{0}^{z}\boldsymbol{B}(x,y,z^{\prime })\,\text{d}z^{\prime },\end{eqnarray}$$

but uses the alternative boundary term

(5.5) $$\begin{eqnarray}\boldsymbol{A}_{0}^{c}=-\frac{1}{2}\left(\int _{-1}^{y}B_{z}(x,y^{\prime },0)\,\text{d}y^{\prime }\right)\hat{\boldsymbol{x}}+\frac{1}{2}\left(\int _{-1}^{x}B_{z}(x^{\prime },y,0)\,\text{d}x^{\prime }\right)\hat{\boldsymbol{y}}.\end{eqnarray}$$

There is nothing special about these chosen gauges, other than their simplicity of computation. The latter has led to the use of $\boldsymbol{A}_{p}^{a}$ and $\boldsymbol{A}_{p}^{b}$ for computation of $H_{R}$ (e.g. Valori et al. Reference Valori, Pariat, Anfinogentov, Chen, Georgoulis, Guo, Liu, Moraitis, Thalmann and Yang2016), where, of course, the choice of gauge does not affect the result, except for numerical error.

Figure 4 shows maps of ${\mathcal{A}}$ , ${\mathcal{A}}_{p}$ , and ${\mathcal{A}}_{R}$ for each of the four gauges, computed at resolution $256\times 256$ from the datacube with resolution $\unicode[STIX]{x0394}x=2/128$ . These maps show only the $z=0$ boundary. For the first three gauges, $\hat{\boldsymbol{n}}\times \boldsymbol{A}$ has been matched to $\hat{\boldsymbol{n}}\times \boldsymbol{A}_{p}$ using the gauge-matching procedure described in § 4. For  $\boldsymbol{A}_{p}^{\ast }$ , this is not necessary as both $\boldsymbol{A}$ and $\boldsymbol{A}_{p}$ are independently matched to the minimal gauge using the procedure in § 4.2.

The first important observation is that all three distributions depend on the gauge choice $\hat{\boldsymbol{n}}\times \boldsymbol{A}_{p}$ . In fact, the results from $\boldsymbol{A}_{p}^{a}$ and especially $\boldsymbol{A}_{p}^{b}$ are quite close to $\boldsymbol{A}_{p}^{\ast }$ in this case. This is because, for this specific  $\boldsymbol{B}$ , the majority of $|B_{n}|$ is located on the lower boundary, where $\boldsymbol{A}_{0}^{a}$ and $\boldsymbol{A}_{0}^{b}$ satisfy the solenoidal condition $\unicode[STIX]{x1D735}_{h}\boldsymbol{\cdot }\boldsymbol{A}=0$ . The gauge $\boldsymbol{A}_{p}^{c}$ (figure 4 g–i) lacks this property, and gives more radically different results from  $\boldsymbol{A}_{p}^{\ast }$ . Notice that $\boldsymbol{A}_{p}^{c}$ leads to larger maximum values of  $|{\mathcal{A}}_{R}|$ , and that these are more localised than either the regions of strongest $|{\mathcal{A}}_{R}^{\ast }|$ or of strongest $J_{\Vert }$ (figure 3 c).

Figure 4.

Field-line helicity on $z=0$ for the Low–Lou example in four different gauges (top to bottom). The left and middle panels show ${\mathcal{A}}$ for the original and potential fields, respectively, while the right panel shows their difference ${\mathcal{A}}_{R}$ in each case, which is the relative field-line helicity. The definitions of the vector potentials are given in the text, but note that (j–l) show the minimal gauge  $\boldsymbol{A}^{\ast }$ . Dashed lines show contours of $B_{z}=B_{pz}$ for context. All of the field-line helicities are normalised by the boundary flux $\unicode[STIX]{x1D6F7}_{0}=(1/2)\int _{z=0}|B_{n}|\,\text{d}^{2}x$ .

Notice that the field-line helicity of the reference field, ${\mathcal{A}}_{p}$ , is weaker than ${\mathcal{A}}$ in both $\boldsymbol{A}^{a}$ and  $\boldsymbol{A}^{b}$ , and even more so in  $\boldsymbol{A}^{\ast }$ . This is a good indication that these are ‘sensible’ gauge choices, since $\boldsymbol{B}_{p}$ is the minimum-energy field and cannot have significantly twisted field lines (it can, however, still exhibit topological structure due to the non-uniform boundary conditions $B_{pn}$ , cf. Prior & Yeates Reference Prior and Yeates2014). Similarly, although ${\mathcal{A}}^{\ast }$ is primarily located within the main current-carrying region (cf. figure 3 c), it does not completely vanish outside of this region. The boundary of this current-carrying region appears as a sharp separatrix between field lines with both ends on the lower boundary  $z=0$ and field lines with only one end on this boundary. Notice that the location of this separatrix differs between $\boldsymbol{B}$ and  $\boldsymbol{B}_{p}$ , so that the relative field-line helicity map ${\mathcal{A}}_{R}$ inherits discontinuous jumps from both ${\mathcal{A}}$ and ${\mathcal{A}}_{p}$ .

As a validation exercise, we have verified that the all of the ${\mathcal{A}}_{R}$ in figure 4 integrate to give the same  $H_{R}$ , as they should by definition. Figure 5 shows the resulting $H_{R}$ for each of the four gauges, as a function of the mesh spacing $\unicode[STIX]{x0394}x$ for the four datacubes. Solid curves show volume integrals of $(\boldsymbol{A}+\boldsymbol{A}_{p})\boldsymbol{\cdot }(\boldsymbol{B}-\boldsymbol{B}_{p})$ , using the composite trapezium rule, whereas dashed curves show estimates of $H_{R}$ from integrating ${\mathcal{A}}_{R}$ using (3.1). For the latter computations, field lines were traced from a mesh with spacing $\unicode[STIX]{x0394}x\times \unicode[STIX]{x0394}x$ on each of the six boundary faces, and the integral then estimated with the composite trapezium rule. It is clear from figure 5 that taking a finer mesh leads to convergence to a common value of $H_{R}$ between the different gauges, whether computed by volume integration or from  ${\mathcal{A}}_{R}$ . Moreover, this common value is consistent with the computations of Valori et al. (Reference Valori, Pariat, Anfinogentov, Chen, Georgoulis, Guo, Liu, Moraitis, Thalmann and Yang2016, see their figure 7a).

Figure 5.

Numerically estimated relative helicity $H_{R}$ for the Low–Lou example, as a function of the numerical mesh spacing, using four different gauges as shown in the legend. In each case, solid lines/circles show direct volume-integrated $H_{R}$ using the composite trapezium rule, while dashed lines/pluses are computed by integrating ${\mathcal{A}}_{R}|B_{n}|/2$ over all six boundary faces. The value of $H_{R}$ is normalised by the square of the boundary flux  $\unicode[STIX]{x1D6F7}_{0}$ .

5.2 Titov–Démoulin equilibria

More realistic models of non-potential solar active regions are given by the Titov–Démoulin family of equilibria (Titov & Démoulin Reference Titov and Démoulin1999), in which the free magnetic energy is localised to a toroidal current channel contained within a surrounding potential field. We utilise datacubes of four such equilibria computed and kindly shared by Valori et al. (Reference Valori, Pariat, Anfinogentov, Chen, Georgoulis, Guo, Liu, Moraitis, Thalmann and Yang2016). For details of the specific construction of these equilibria, see that paper. Here we just note that the domain is the Cartesian volume $V=[-3.18,3.18]\times [-5.10,5.10]\times [0,4.56]$ , and that the primary difference between the four cases is the twist of magnetic field lines in the current channel. The cases are denoted $N=0.1,0.5,1,3$ , where $N$ is the (approximate) number of twists made by a field line along the current channel. Corresponding magnetic field-line plots are shown in figure 6(a,d,g,j). Note that the boundary conditions $B_{n}$ and potential magnetic energy also differ between the cases. The relative free energies (ratio of non-potential to total magnetic energy) are ${\approx}0.002$ , 0.03, 0.20 and 0.15, respectively. In all computations here, we use datacubes with $\unicode[STIX]{x0394}x=0.06$ .

Figure 6.

The Titov–Démoulin solution for increasing (top to bottom) values of the twist parameter  $N$ , showing the magnetic field lines, field-line integrated parallel current $J_{\Vert }/B_{0}$ and relative field-line helicity ${\mathcal{A}}_{R}^{\ast }$ in each case. The field lines are coloured by  $J_{\Vert }$ , with the same colour scale as the middle column. The colour scales are fixed; actual maxima of $|J_{\Vert }/B_{0}|$ are 4.6, 21.2, 53.3, 199.4, and of $|{\mathcal{A}}_{R}^{\ast }/\unicode[STIX]{x1D6F7}_{0}|$ are 0.03, 0.13, 0.31, 0.32, respectively.

Figure 7.

Field-line helicity on $z=0$ for the $N=0.5$ Titov–Démoulin model in four different gauges (top to bottom). The format is the same as figure 4.

We begin our analysis with the same computation of ${\mathcal{A}}$ , ${\mathcal{A}}_{p}$ and ${\mathcal{A}}_{R}$ in four different gauges as for the Low–Lou equilibrium. These are shown in figure 7 for one of the datacubes: $N=0.5$ . Once again, the results are clearly dependent on the gauge. However, in this case, the downward DeVore–Coulomb gauge ( $\boldsymbol{A}_{p}^{b}$ ) gives very similar results to the minimal gauge ( $\boldsymbol{A}_{p}^{\ast }$ ), whereas the upward DeVore–Coulomb gauge ( $\boldsymbol{A}_{p}^{a}$ ) gives rather different results, along with  $\boldsymbol{A}_{p}^{c}$ . In all four cases, ${\mathcal{A}}$ agrees well with ${\mathcal{A}}_{p}$ outside of the toroidal current channel (which appears as two tear-drop shaped footpoints on this boundary). However, for $\boldsymbol{A}_{p}^{b}$ and $\boldsymbol{A}_{p}^{\ast }$ , the current channel has a much higher strength of $|{\mathcal{A}}|$ than outside, whereas for $\boldsymbol{A}_{p}^{a}$ and $\boldsymbol{A}_{p}^{c}$ , comparable strengths of $|{\mathcal{A}}|$ are found at locations both inside and outside the channel.

In this example, the relative field-line helicity ${\mathcal{A}}_{R}^{\ast }$ (and also ${\mathcal{A}}_{R}^{b}$ ) is much more concentrated than in the Low–Lou equilibrium. This holds for all  $N$ , as shown in figure 6. It is strongest in the tear-drop shaped footpoints of the current channel, although diffuse non-zero values are found throughout the ‘closed field’ region (i.e. the region between the diagonal separatrices within which field lines have both ends on $z=0$ and cross the central polarity-inversion line). There are also thin ridges in ${\mathcal{A}}_{R}^{\ast }$ along the separatrices; these arise because the separatrices are in slightly different locations in $\boldsymbol{B}$ and  $\boldsymbol{B}_{p}$ . Outside of these separatrices, ${\mathcal{A}}_{R}^{\ast }$ is close to zero, reflecting the fact that these field lines are close to potential. Figure 8 shows a comparison of field lines of $\boldsymbol{B}$ and $\boldsymbol{B}_{p}$ for two of the datacubes, showing this effect.

Figure 8.

Field lines of $\boldsymbol{B}$ (purple) and $\boldsymbol{B}_{p}$ (green), traced from the same starting points in $y>0$ (to left in this view), for the Titov–Démoulin examples with (a)  $N=0.1$ and (b)  $N=1$ .

If we compare the maximum magnitude of $|{\mathcal{A}}_{R}^{\ast }|$ as $N$ is increased (figure 6), we see that the maximum $|{\mathcal{A}}_{R}^{\ast }|$ (relative to  $\unicode[STIX]{x1D6F7}_{0}$ ) increases from $N=0.1$ to 0.5 and again from $N=0.5$ to 1, but is comparable between $N=1$ and 3, despite the fact that the maximum $J_{\Vert }/B_{0}$ increases by a factor 4. This reflects the fact that ${\mathcal{A}}_{R}^{\ast }$ is a global topological measure that does not solely measure the local twist (or ‘self-helicity’) of magnetic field lines. However, it does highlight that $|{\mathcal{A}}_{R}^{\ast }|$ will likely not relate directly to the stability of a given configuration, since the $N=3$ case is expected to be kink unstable while the $N=1$ (and $N=0.1$ , 0.5) is expected to be stable (Valori et al. Reference Valori, Pariat, Anfinogentov, Chen, Georgoulis, Guo, Liu, Moraitis, Thalmann and Yang2016). The similar peak values of $|{\mathcal{A}}_{R}^{\ast }|$ located in a smaller footprint for $N=3$ (as seen in figure 6 k) is consistent with this case having a lower $|H_{R}|$ than $N=1$ , as found by Valori et al. (Reference Valori, Pariat, Anfinogentov, Chen, Georgoulis, Guo, Liu, Moraitis, Thalmann and Yang2016).

Finally, we repeat the validation exercise of computing $H_{R}$ , as for the Low–Lou model. Figure 9(a) shows the results in the same format as figure 5, except that the horizontal axis now labels the cases with different $N$ (measured by the end-to-end twist, to facilitate comparison with figure 3a of Valori et al. Reference Valori, Pariat, Anfinogentov, Chen, Georgoulis, Guo, Liu, Moraitis, Thalmann and Yang2016). As for the Low–Lou example, we find consistent results when $H_{R}$ is computed by volume integration. When computed from ${\mathcal{A}}_{R}$ via (3.1), we obtain consistent values for $H_{R}$ for $N=0.1$ , 0.5 and 3, but $|H_{R}|$ is underestimated by up to 5 % for $N=1$ , depending on the gauge. Further investigation shows that this is not due to the mesh resolution (of either the datacube or the mesh for computing ${\mathcal{A}}_{R}$ on the boundary), but is caused by the fact that the initial data do not perfectly satisfy the solenoidal condition on our staggered mesh. To show this, figure 9(b) shows the same computations of $H_{R}$ when we include an additional correction to $\boldsymbol{B}$ after computing the (initial) $\boldsymbol{A}$ on the staggered mesh. This correction recomputes $\boldsymbol{B}$ from the numerical curl of  $\boldsymbol{A}$ , so that the solenoidal condition is satisfied to machine precision, but at the expense of changing  $\boldsymbol{B}$ . With this correction, figure 9(b) shows that the estimate of $H_{R}$ from ${\mathcal{A}}_{R}$ then agrees with the volume integral, in all four gauges. Notice however that $|H_{R}|$ is then overestimated in  $\boldsymbol{A}^{b}$ , compared with the other gauges. This arises because the solenoidal correction has a greater effect on the current-carrying region in the downward DeVore gauge, since this region is located at lower  $z$ . Thus this affects the estimated helicity more. This illustrates the importance of the solenoidal condition in estimating helicity (Valori et al. Reference Valori, Pariat, Anfinogentov, Chen, Georgoulis, Guo, Liu, Moraitis, Thalmann and Yang2016).

Figure 9.

Numerically estimated relative helicity $H_{R}$ for the Titov–Démoulin example, as a function of the end-to-end twist (in radians). As in figure 5, colours denote different gauges, and solid lines/circles show direct volume integrated $H_{R}$ using the composite trapezium rule, while dashed lines/pluses are computed by integrating ${\mathcal{A}}_{R}|B_{n}|/2$ ${\mathcal{A}}_{R}$ over all six boundary faces. In (b), an additional divergence-correction step is applied, as described in the text.