© ESO 2020

1. Introduction

The concept of magnetic helicity for a divergence free field B is most commonly introduced as the following scalar integral quantity

$$\begin{array}{c}\hfill H={\displaystyle {\int }_V}\mathit{A}·\mathit{B}\mathrm{d}V,\end{array}$$(1)

where A is the vector potential of B (B = ∇ × A). This measure was originally introduced by Woltjer (1958), and given a topological definition by Moffatt (1969) as the linking of magnetic field lines (see also Arnol’d & Khesin 1998 when field lines do not form closed curves). If we decompose a magnetic field into distinct magnetic regions (by distinct, we mean that fieldlines do not cross the boundaries of the regions within the volume V), then helicity can be decomposed into the sum of self helicities of each region, and mutual helicities between regions (Berger 1999). For example, if the regions are flux ropes, then the self helicity can be described as the twist and writhe of individual ropes, while the mutual helicity arises from the linking or braiding of distinct ropes. This decomposition has also been applied to the studies of coronal loops: see Aschwanden (2019), where the authors investigate how the stability of coronal loops is associated with the braided linking number. Other shapes are possible: for example an arcade in the solar corona can be sheared (self helicity), and it can also envelop a flux rope (mutual helicity).

Magnetic helicity is a well conserved quantity in low resistivity magnetohydrodynamics (Taylor 1974; Moffatt 2018). The conservation is maintained in less ideal conditions, albeit to a weaker degree (Berger 1984), making it an ideal approximate invariant for investigation into complex magnetic field systems (Ji et al. 1995; Brandenburg 2009; Contopoulos et al. 2009; Russell et al. 2015; Zuccarello et al. 2018).

Magnetic helicity plays an important role in studies of magnetohydrodynamic (MHD) turbulence in general, and dynamo theories of magnetic energy generation in particular (e.g. Vishniac & Cho 2001; Blackman & Brandenburg 2003; Sur et al. 2007; Brandenburg et al. 2016). In a two scale kinematic dynamo, the large scale energy can increase exponentially. This poses a problem for magnetic helicity conservation. If the large scale magnetic helicity increases exponentially, then the small scale field must have an equal and opposite helicity which also blows up. Dissipation of the small scale helicity may not be physically feasible.

A solution to this problem lies in making the dynamo inhomogeneous – the dynamo operates in one region of space (e.g. the base of the convection zone) and excess magnetic helicity is carried away (Brandenburg 2009; Vishniac & Shapovalov 2013). However, to model this process properly, we need to be able to specify how helicity is spatially distributed. In other words we need to be able to locate where helicity resides more precisely than simply using complete flux ropes.

Another area where helicity localisation could be useful is in the study of solar activity. Many studies show how magnetic helicity can flow from the interior into active regions (e.g. Berger & Ruzmaikin 2000; Kusano et al. 2002; Pevtsov 2003; Park et al. 2008; Dalmasse et al. 2014; Prior & MacTaggart 2019). A knowledge of how this helicity is distributed within the active region may aid the understanding and prediction of flares and coronal mass ejections. Scale dependence of magnetic helicity can also help in understanding the evolution of turbulence in the solar wind (Brandenburg et al. 2011).

Localising helicity is difficult for a number of reasons. First, while one can attempt to define helicity density as the quantity

$$\begin{array}{c}\hfill H_{\mathrm{den}}=\mathit{A}·\mathit{B},\end{array}$$(2)

this is in no sense gauge invariant as gradient fields can be added to A without changing the magnetic field. Second, it does not represent spatially localised information, in a particular gauge such as Coulomb gauge, A is an integral over the magnetic field, and is thus non-local. This has mathematical grounding – helicity is associated with the Gauss linking number, for which we must take a double integral across all space. If we only have information about a small patch of field, there is no way of knowing how a field line goes on to twist and writhe around the rest of the field.

What about integrals of helicity as in Eq. (1)? For the expression (1) to be physically meaningful, V must either be unbounded space, or if V is finite with boundary S then S must be a either a magnetic surface $(\mathit{B}·\widehat{\mathit{n}}{|}_S=0$), or the field must have periodic boundary conditions on S. If, for example, there is net flux perpendicular to two periodic directions this can lead to unphysical effects involving magnetic helicity (Berger 1997; Watson & Craig 2001; Brandenburg & Subramanian 2004).

In Sect. 2 we formally introduce the gauge problem associated with magnetic helicity, and briefly describe relative helicity, which gives a gauge-invariant measure when the volume is not bounded by a magnetic surface. In general this measure is not additive in the sense that the helicity of all space may not equal the sum of helicities of subvolumes. We then discuss Fourier decompositions of helicity, which help to provide information on how helicity behaves on various scales. Under certain circumstances (isotropic turbulence) this approach can overcome both the gauge and localisation problems; but not in the inhomogeneous aperiodic cases we would like to generalise to. This includes a discussion on the transformed two-point correlation tensor which decomposes into helicity and energy. Finally we discuss field line helicities which measure how one chosen field line interlinks with all other lines. This quantity can be used to accurately quantify reconnection activity in magnetic fields (Prior & Yeates 2018), however, the decomposition of helicity into contributions from individual fieldlines is still not localised.

In Sect. 3 we present helicity densities expressed as two-point correlation functions which are the building blocks of linking and winding. This final measure is gauge invariant, even for fields with non-trivial boundary distributions, as it does not depend on the vector potential for its definition. We select this measure as our base definition and in subsequent sections show it can be used to overcome the spatial localisation problem.

Section 4 provides a background to wavelet transforms and multi-resolution analysis as a solution to the localisation problem, and how these can be applied to helicity integrals. In Sect. 5 we give a formal introduction of the full 3D wavelet transform and its application to the helicity integral. Section 6 provides examples of how the wavelet multi-resolution helicity formulation can be applied in practice. This includes a pair of twisted flux ropes which present a trivial (null) spectral decomposition; the multi-resolution helicity decomposition is shown to resolve the spatial separation of the system’s entanglement. A second example of a pair of interlinked twisted flux ropes demonstrates how the decomposition can separate out the contributions from large scale linking and smaller scale twisting, as well as correctly assess the localisation of the helicity in this system. In Sect. 7 we consider the application of the multi-resolution wavelet decomposition to our geometric two-point correlation definition of helicity. This is used to derive linear helicity-energy decompositions for both the helicity and the field line helicity. Section 8 considers an example of a reconnecting magnetic braid, based on the numerical experiments in Wilmot-Smith et al. (2009, 2011) and Russell et al. (2015). The field line helicity multi-resolution analysis is utilised here. In particular we show that the field’s twisted structure and its field line entanglement balance their helicity fluxes at differing spatial scales. Further we show that the growth and then decay in magnetic energy of this system highly correlated with the field line helicity relaxation at that the dominant spatial scales. In Sect. 9 we apply the multi-resolution decomposition to a flux transport model and finally conclude in Sect. 10.

2. Existing helicity decompositions

2.1. The gauge problem

Suppose the volume of interest V is neither bounded by a magnetic surface or of infinite extent. This introduces a gauge dependence to the helicity integral: given some function Φ we can let A_G = A + ∇Φ, which induces a change in helicity corresponding to

$$\begin{array}{c}\hfill H_{\mathrm{G}}=H+{\displaystyle {\int }_S}\mathrm{\Phi }\mathit{B}·\mathrm{d}\mathit{S}.\end{array}$$(3)

To overcome this problem, Berger & Field (1984) introduced a gauge invariant measure of helicity known as “relative helicity”, which measures the magnetic helicity of the magnetic field B relative to some secondary field B₀ by taking the difference between their helicities across all space. However, this quantity does not have the property of additivity; if V is decomposed into subvolumes, the relative helicity from each subvolume may not equal the total in V. Further, the reference field may not be smooth across boundaries of sub-domains.

2.2. Fourier spectra

The splitting of magnetic fields into different scales is already core to the study of many magnetohydrodynamical systems: Verma (2004) provides an in-depth review of turbulent magnetohydrodynamic fields, which have energy interchanges occurring across a spectrum of spatial scales. Following Blackman (2004, 2015), and Subramanian & Brandenburg (2005), we can write the magnetic energy spectrum as

$$\begin{array}{c}\hfill E_k={\displaystyle \int {|\stackrel{\sim }{\mathit{B}}|}^2}{k}^2\mathrm{d}{\mathrm{\Omega }}_k,\end{array}$$(4)

where k = ||k||, and a tilde represents the Fourier transform. Ω_k then represents a spherical shell in wave space, given by all wavenumbers k₋ ≤ k ≤ k₊, for which k_± = ||k|| ± 0.5. In Fourier space, we have the direct relation $\stackrel{\sim }{\mathit{A}}=-i\mathit{k}\times \stackrel{\sim }{\mathit{B}}/k^2$. Thus we can write

$$\begin{array}{c}\hfill H_k={\displaystyle \int }i(\mathit{k}\times {\stackrel{\sim }{\mathit{B}}}^{\ast }(\mathit{k}))·\stackrel{\sim }{\mathit{B}}(\mathit{k})\mathrm{d}{\mathrm{\Omega }}_k,\end{array}$$(5)

and as such we have a gauge invariant measure of magnetic helicity at scale L = 2π/k which has the property of additivity (see for example Moffatt 1978; Blackman & Brandenburg 2003; Démoulin 2007; Brandenburg et al. 2016). The gauge invariance follows from the fact that the surface integral of (3) vanishes for the periodic boundary conditions required for a Fourier representation of B.

It is important to note that the Fourier decomposition can produce spurious results: if we imagine an infinite system of alternately twisted flux tubes, the Fourier transform of magnetic helicity would be zero at all scales (Asgari-Targhi & Berger 2009). To see this, suppose that B_z is constant so only has power at k = 0, but B_x and B_y vary in x and y to make the oppositely twisted tubes. Then at any k >  0, both k and $\stackrel{\sim }{\mathit{B}}(\mathit{k})$ will be in the x − y plane. Thus the triple product must involve B_z; but B_z will be zero for k >  0.

Also the Fourier spectrum does not give information on spatial locality. The windowed Fourier transform can help. An envelope function with compact support is convolved on top of the infinite sinusoidal functions. Taking the Fourier transform using such a reduced analytic form gives an idea of the variations corresponding to scales at a given locality, but has two downsides (aside from the requirement of periodicity): firstly, the transformation does not provide an orthogonal basis, which is required to maintain additivity. Secondly, the window size is fixed, meaning we cannot separate intense fluctuations which are on smaller scales than the window size from weak contributions on the same scale as the window size.

2.3. Two-point correlation functions

An important consequence of a Fourier decomposition of a magnetic field is that the helicity H_k can be related to the magnetic energy E_k via the transform of the two-point correlation tensor M_ij:

$$\begin{array}{c}\hfill M_{\mathit{ij}}(\mathit{X},\mathit{x})=B_i(\mathit{X}-\mathit{x})B_j(\mathit{X}+\mathit{x}).\end{array}$$(6)

In a periodic domain one can transform this function over the displacement x to obtain a skew-symmetric tensor function M̃_ij(X,k) of both position and scale, and further, for isotropic turbulence, this can be decomposed (in three dimensions) as

$$\begin{array}{c}\hfill {\stackrel{\sim }{M}}_{\mathit{ij}}=[({\delta }_{\mathit{ij}}-k_i^u{k}_j^u)2E_k-\mathrm{i}{k}_l^u{ϵ}_{\mathit{ijl}}k{H}_k]/8k^2\pi ,\end{array}$$(7)

where $k_i^u$ is the i^th component of the unit vector of k and ϵ_ijl the alternating tensor (Roberts & Soward 1975; Brandenburg et al. 2016). So the energy is the trace of the tensor M̃_ij and the helicity represented by the off-diagonal components. This is a potentially powerful relation between the magnetic helicity and energy on a given Fourier shell at each point of space.

In this article, we intend to provide a similar decomposition of magnetic helicity which preserves additivity and scale dependence, whilst also providing information about the spatial locality of terms contributing to the power at each scale. Key to our study is the lack of any assumptions about the boundary conditions of the magnetic fields or any isotropic assumptions. One result is a variant of (7) which can retain information on the spatial distribution of the magnetic energy and helicity, even in highly inhomogeneous systems.

2.4. Fieldline helicity

Field line helicity is another tool that has become more popular in recent years. For a given field line γ we have (Berger 1988; Yeates & Hornig 2013; Prior & Yeates 2014; Yeates & Page 2018; Moraitis et al. 2019)

$$\begin{array}{c}\hfill \mathcal{A}(\gamma )={\displaystyle {\int }_{\gamma }}\mathit{A}·\mathit{T}\mathrm{d}s,\end{array}$$(8)

where T = B/||B|| is the unit tangent vector along the fieldline, and s is the arclength parameter of its curve. The fieldline helicity measures the average winding of all field lines around the field line under analysis (Prior & Yeates 2014). Field line helicity can be seen as the limit of the methodology of Pevtsov (2003), where each magnetic surface encloses exactly one field line. If we imagine tracing the field lines between two planes, the field line helicity associated with a field line starting at each point (x_l, y_l) on some initial plane (typically taken as the lower boundary of a system) gives a two-dimensional density.

The fieldline helicity is linear in A so has the property of additivity. It is not, however, gauge invariant, nor does it have the property locality as the expression for A involves an integration over at least one spatial dimension (see Prior & Yeates 2014). There is a relative field line helicity version of this quantity, whose definition comes attached with various technical complexities (Yeates & Page 2018; Moraitis et al. 2019), but is an invariant for each individual field line. Further, there is some remaining gauge dependence.

2.5. The need for a different baseline definition of helicity

These decompositions of helicity so far all have drawbacks; none combine the properties of gauge independence, spatial locality, and additivity we seek. To overcome this we begin by adopting a purely geometrical definition of helicity. We first show this has no gauge dependence, then in the following sections we show how it can be used to overcome the additivity and locality problems.

3. Magnetic winding and gauge independence

3.1. Helicity is winding

Given any integral representation for A, the helicity integral becomes a double integral involving B evaluated at two different points. For example, in Coulomb gauge with A expressed via the Biot-Savart law,

$$\begin{array}{c}\hfill H=-\frac{1}{4\pi }{\displaystyle {\int }_V}{\displaystyle {\int }_{V^{\prime }}}\frac{\mathit{B}·{\mathit{B}}^{\prime }\times \mathit{r}}{r^3}{\mathrm{d}}^{3}x{\mathrm{d}}^3{x}^{\prime }.\end{array}$$(9)

The integrand can be regarded as a two-point correlation function for the magnetic field (Subramanian & Brandenburg 2005).

A particular topologically meaningful choice for A is the winding gauge. Prior & Yeates (2014) considered the winding gauge A^w,

$$\begin{array}{cc}& {\mathit{A}}^w(x,y,z)=\frac{1}{2\pi }{\displaystyle {\int }_{S_z}}\frac{\mathit{B}(x^{\prime },y^{\prime },z)\times \mathit{r}}{r^2}\mathrm{d}{x}^{\prime }\mathrm{d}{y}^{\prime }\hfill \\ \hfill & \mathit{r}=(x-x^{\prime },y-y^{\prime },0),\hfill \end{array}$$(10)

where S_z is a plane of constant z value. With this choice the helicity can be written as

$$\begin{array}{c}\hfill H=\frac{1}{2\pi }{\displaystyle {\int }_{z_0}^{z_1}}{\displaystyle {\int }_{S_z}}{\displaystyle {\int }_{S_z^{\prime }}}\frac{\mathit{B}·\mathit{B}({\mathit{x}}^{\mathbf{\prime }},{\mathit{y}}^{\mathbf{\prime }},\mathit{z})\times \mathit{r}}{r^3}\mathrm{d}{x}^{\prime }\mathrm{d}{y}^{\prime }\mathrm{d}x\mathrm{d}y\mathrm{d}z,\end{array}$$(11)

for a two-point correlation within the plane S_z of constant z. Here the set of planes S_z, z ∈ [z₀, z₁] cover the whole domain V. Prior & Yeates (2014) showed that (11) is just the flux-weighted average winding of all pairs of field line of B with each other, that is

$$\begin{array}{cc}\hfill H(\mathit{B})=\frac{1}{2\pi }{\displaystyle {\int }_{z_0}^{z_1}}{\displaystyle {\int }_{S_z}}& {\displaystyle {\int }_{S_z^{\prime }}}{B}_z(\mathit{x})B_z({\mathit{x}}^{\mathbf{\prime }})\frac{\mathrm{d}}{\mathrm{d}z}\mathrm{\Theta }(\mathit{x},{\mathit{x}}^{\mathbf{\prime }}){\mathrm{d}}^2\mathit{x}{\mathrm{d}}^2{\mathit{x}}^{\mathbf{\prime }}\mathrm{d}z\hfill \end{array}$$(12)

In any planar slices S_z, $S_z^{\prime }$, Θ(x, x′) defines the “angle” between the two fieldlines centred at x and x′,

$$\begin{array}{c}\hfill \mathrm{\Theta }(\mathit{x},{\mathit{x}}^{\mathbf{\prime }})=arctan\left(\frac{y-y^{\prime }}{x-x^{\prime }}\right),\end{array}$$(13)

We demonstrate this pictorially in Fig. 1 for two curves γ and γ′. Similarly a visualisation of the topological nature of the integrand of (11), a two-point correlation function which measures the net winding, is shown in Fig. 2.

Fig. 1. Winding number interpretation of helicity. The winding is defined by the mutual angle Θ between two curves γ and $\stackrel{\sim }{\gamma }$. The yellow arrows depict a fixed direction and the black arrows the joining vector of the two curves in a given plane Sz used to define Θ.

Fig. 2. Two point correlation interpretation of helicity. Given two field lines, the product of Bz for the first field line, and $B_{\theta }^{\prime }$ (in the plane) for the second line, as would be in r × B′, gives a measure of their mutual winding.

We remark that this requires that the field can be composed of a set of planar surfaces 𝒱 = {S_z|z ∈ [z₀,z₁]} and that, if the volume is finite in x or y, then the field B be tangent on the side surfaces. Berger & Hornig (2018) showed that this relation can be obtained from a poloidal-toroidal decomposition, and extended it to more general domains which can be constructed from sets of simply connected surfaces. Further it was shown in Prior & Yeates (2014) that this choice is preferable in that any other choice of gauge or reference field gives a helicity measure which is equivalent to choosing to measure the angle Θ with respect to a varying direction, whose rotation is non physical in that it does not relate to the entanglement of the field itself.

Prior & Yeates (2014) also showed that the field line helicity can be written as

$$\begin{array}{cc}& \mathcal{A}(\gamma )=\frac{1}{2\pi }{\displaystyle {\int }_{\gamma }}{\displaystyle {\int }_{S_z^{\prime }}}\frac{\mathit{T}·\mathit{B}(x^{\prime },y^{\prime },z(s))\times {\mathit{r}}_{\gamma }}{r_{\gamma }^2}\mathrm{d}{x}^{\prime }\mathrm{d}{y}^{\prime }\mathrm{d}s,\hfill \\ \hfill & {\mathit{r}}_{\gamma }=(x^{\prime }-{\gamma }_x,y^{\prime }-{\gamma }_y,0),\hfill \end{array}$$(14)

if the winding gauge is chosen. This represents the entanglement of the field line γ with the rest of the field. In Sect. 7.1 we show that, using a wavelet decomposition of B, it can be represented as a spatial sum of skew symmetric tensors whose trace give the magnetic energy and off-diagonal elements give the helicity, similar to the two-point correlation Fourier transform relationship (7).

3.2. A gauge independent measure of magnetic helicity

The crucial point about (11) is that it gives a definition of a quantity which measures the field-weighted entanglement of the magnetic field, and depends only on the field B (as does (14) for the fieldline helicity). To be sure, we can relate it to the classical magnetic helicity definition (H = ∫_VA ⋅ B) via a vector potential, but the following properties can be ascribed to the quantity purely on the basis of its topological definition in terms of winding rate dΘ/dz. Firstly, it is invariant under ideal evolutions which vanish on the domain boundaries (Prior & Yeates 2014). Secondly, it is approximately conserved for low plasma β relaxations (Russell et al. 2015). Thirdly, the field line helicity density can be used to directly quantify magnetic reconnection (Prior & Yeates 2018), even for fields with normal boundary components.

These are all the properties that mark the magnetic helicity as a fundamental quantity in solar physics applications. None of them rely on a vector potential definition to be applicable (as demonstrated in the indicated references). In what follows we assert the two-point correlation function as our fundamental definition of the magnetic helicity. For the sake of clarity we formally define the following two-point correlation integral

$$\begin{array}{c}\hfill \mathit{C}(x,y,z)={\displaystyle {\int }_{S_z^{\prime }}}\frac{\mathit{B}(x^{\prime },y^{\prime },z)\times \mathit{r}}{r^2}\mathrm{d}{x}^{\prime }\mathrm{d}{y}^{\prime }.\end{array}$$(15)

The product B ⋅ C represents the total winding and field weighted correlation of the field at a point (x, y) in the plane S_z with the whole field in that plane (via Eq. (11)). If the field is tangent on the boundaries of the plane S_z then C = A^w, but as we have discussed above C is a meaningful topological quantity even if this is not so.

We then have the following gauge free definitions of the helicity and field line helicity which place no requirements on the system’s boundary conditions

$$\begin{array}{cc}& H(\mathit{B})={\displaystyle {\int }_0^h}{\displaystyle {\int }_{S_z}}\mathit{B}·\mathit{C}\mathrm{d}x\mathrm{d}y\mathrm{d}z,\hfill \end{array}$$(16)

$$\begin{array}{cc}& \mathcal{A}(\mathit{B})={\displaystyle {\int }_{\gamma }}\mathit{T}·\mathit{C}(\gamma )\mathrm{d}s,\hfill \end{array}$$(17)

where z ∈ [0, h] is the z-coordinate which labels the set of planes S_z of the Cartesian domain on which the field is defined.

What remains is to decompose these quantities such that we have additivity and spatial locality.

4. Spatial localisation and additivity of the helicity

4.1. Localised field decompostion

Here we very briefly introduce the ideas underpinning our spatial localisation technique to give some geometrical intuition as to its interpretation. In order to spatially decompose the helicity we need a representation of the field B in terms of functions with compact support. For example the box function

$$\begin{array}{c}\hfill \varphi (x)=\{\begin{array}{cc}1& \mathrm{if}0\le x\le 1,\hfill \\ 0& \mathrm{if}x>1\text{or}x<0,\hfill \end{array}\end{array}$$(18)

whose 3D composite

$$\begin{array}{c}\hfill {\mathrm{\Phi }}_{x_0{y}_0{z}_0}(x,y,z)=\varphi (x-x_0)\varphi (y-y_0)\varphi (z-z_0),\end{array}$$(19)

gives a box of compact support which is translated in space. By discretizing the domain as x₀ = iΔx,  i ∈ 1, …N, and suitably scaling (18), a set of box functions Φ_lmn is created to cover the domain in a non overlapping fashion (e.g. Fig. 3a). We can then approximate the mean field B as

$$\begin{array}{cc}& \mathit{B}(x,y,z)\approx {\displaystyle \sum _{l=1}^N\sum _{m=1}^N\sum _{n=1}^N}{\mathit{B}}_{\mathit{lmn}}{\mathrm{\Phi }}_{\mathit{lmn}}(x,y,z),\hfill \end{array}$$(20)

$$\begin{array}{cc}& {\mathit{B}}_{\mathit{lmn}}={\displaystyle {\int }_V}\mathit{B}{\mathrm{\Phi }}_{\mathit{lmn}}\mathrm{d}V.\hfill \end{array}$$(21)

Fig. 3. Localised field decompositions. Panel a: a cuboid domain is split into smaller sub-boxes, which each contain a contribution to the total field, as shown in panel b. By making the discretisation even more dense, as shown in panel c, we make our approximation more accurate.

Each coefficient B_lmn is representative of the average behaviour of the field in the box (l, m, n) (Fig. 3b).

We can do something similar for the correlation integral C. Using the fact that the function Φ_lmn has compact support, an approximation to magnetic helicity is then given by the sum

$$\begin{array}{c}\hfill H(\mathit{B})\approx {\displaystyle \sum _{i=1}^N}{\displaystyle \sum _{j=1}^N}{\displaystyle \sum _{k=1}^N}{\mathit{B}}_{\mathit{lmn}}·{\mathit{C}}_{\mathit{lmn}},\end{array}$$(22)

(think of this as a spatial decomposition of the constant part of the Fourier series). Each triplet (lmn) gives the average of the density C ⋅ B in a particular cube of the domain. To capture the local variations, we can use a function such as

$$\begin{array}{c}\hfill \psi (x)=\{\begin{array}{cc}1\hfill & \mathrm{if}0\le x\le 1/2,\hfill \\ -1\hfill & \mathrm{if}1/2<x\le 1,\hfill \\ 0\hfill & \mathrm{if}x>1\mathrm{or}x<0,\hfill \end{array}\end{array}$$(23)

defined for all real x (such that it has compact support). The coefficients of B with respect to this function can then be added to (20) to give a more accurate approximation of the field (this is a little like breaking the sinusoid of the Fourier transform into sub components). The smaller the discretisation size (the spatial scale) N the more accurate the approximation (the discretisation in Fig. 3c would be more accurate than that in a).

Of course there are multiple issues with such a decomposition. For example how do we choose the scale of decomposition? In fact (with regards to the varying component) we might want to choose multiple scales for fields which have both large and small scale variation. How might we then add up these scales whilst avoiding redundancy? A branch of wavelet analysis called Multi-resolution Analysis tells us exactly how to perform such a decomposition orthogonally, and combine it across multiple scales. We shall introduce it formally in Sect. 5: the localised functions used above are the so-called Haar Scaling Function (18) and Haar Wavelet (23), and the sum (20) forms one part of the decomposition; the method for combining the varying field components is a little more complex. We will see in Sect. 6 that this spatial scale decomposition (multi-resolution analysis) of the helicity for the two example fields discussed above show a non trivial (absolute) variation across spatial scales which naturally identifies to “size” and position of the helicity producing components of the field.

4.2. Two-point winding correlation localisation

A decomposition such as (22) is still not fully localised, since the correlation integral C at a point (x, y, z) involves integration across planes of constant z of the domain. Thus the coefficient

$$\begin{array}{c}\hfill {\mathit{C}}_{\mathit{lmn}}={\displaystyle {\int }_V}\mathit{C}{\mathrm{\Phi }}_{\mathit{lmn}}\mathrm{d}V,\end{array}$$(24)

will include integration across all planes S_z containing the points (x, y, z) which have compact support from the function Φ_lmn. As such the quantity C_lmn ⋅ B_lmn represents the winding correlation of the field in the box of compact support of Φ_lmn with the rest of the field in the planes containing that box, as indicated in Fig. 4a.

Fig. 4. Geometric interpretations of spatially decomposed helicity contributions. Panel a: geometrical interpretation of the spatial contribution Clmn ⋅ Blmn of a spatial (wavelet) decomposition of the helicity. The red box represents the spatial sub-domain given by the triplet lmn. Each point in this red domain contributes a winding with the rest of the field in the plane in which it is contained. Because Clmn ⋅ Blmn is a sum over the whole red domain, the set of planes containing the red domain provide winding contributions to the sum, as indicated in the figure. Panel b: winding contribution obtained by spatially decomposing the field B inside the two-point correlation function. It represents the winding of the red curve γ (localised in the plane as represented by the vector B at that point) and the field lines in the sub-plane indicated in green.

A finer localisation can found be inserting the decomposition of the field B directly into the function C. For example, since the integral C is defined in the plane S_z we could use a two dimensional decomposition, that is use functions Φ_lm which approximate the field B as

$$\begin{array}{c}\hfill \mathit{B}\approx {\displaystyle \sum _{l=1}^N\sum _{m=1}^N}{\mathit{B}}_{\mathit{lm}}(z){\mathrm{\Phi }}_{\mathit{lm}}.\end{array}$$(25)

Inserting this into the field line winding integrand density (17) would yield terms in the form

$$\begin{array}{c}\hfill \mathit{T}·\frac{{\mathit{B}}_{\mathit{lm}}(z)\times \mathit{r}}{r^2}·\end{array}$$(26)

This represents the winding of the curve at the point (x, y) (represented by the vector T) with a localised sub domain of curves corresponding to the compact support of Φ_lm, as indicated Fig. 4b. We use this decomposition to develop a spatial helicity energy decomposition for both the helicity and field line helicity in Sect. 7. We stress that this finer decomposition requires an explicit representation of the helicity in terms of B, which our gauge free, physically, and topologically meaningful definition fulfils. This will be our ultimate solution to both the additivity and localisation problems.

5. Helicity, wavelets and multi-resolution analysis

We first consider dimensional scalar functions f(x) as 3D wavelet representations are composed of combinations of one dimensional wavelet decompositions. We will focus on the set of wavelets known as discrete wavelets which form the discrete wavelet transform, in particular the multi-resolution representation of this transform (for details on the continuous wavelet transform see e.g. Grossmann et al. 1990). What follows is far from a comprehensive overview of the mathematics the multi-resolution analysis which is beyond the scope of this study, for a more detailed introduction see, for example Farge (1992) for a practical introduction in a fluid dynamic context and Jawerth & Sweldens (1994) for more details on the underlying mathematics.

5.1. The basic idea for Haar wavelets

In the previous section we introduced the Haar Scaling function (18) and varying Haar Wavelet (23) used to characterise the mean and varying behaviour of a function f(x) respectively, over a given subset of the domain. As discussed a systematic means of choosing the and combining the domains (boxes) of compact support is required. The basic idea of a multi-resolution analysis for some function defined on a domain [0, 1] (one can always scale this to any finite domain) is to choose the domain spatial scales as factors of two 2^s, …s ∈ 1, 2, …, as indicated in Table 1. Then, for a given choice of scale s the functions (23) and (18) are mutually orthogonal and orthonormal with each other if suitably dilated and translated, that is,

$$\begin{array}{cc}& {\displaystyle {\int }_0^1}\sqrt{2^s}\varphi (2^{s}x-i)\sqrt{2^s}\psi (2^{s}x-j)\mathrm{d}x=0,\hfill \end{array}$$(27)

$$\begin{array}{cc}& {\displaystyle {\int }_0^1}\sqrt{2^s}\varphi (2^{s}x-i)\sqrt{2^s}\varphi (2^{s}x-j)\mathrm{d}x\hfill \\ \hfill & ={\displaystyle {\int }_0^1}\sqrt{2^s}\psi (2^{s}x-i)\sqrt{2^s}\psi (2^{s}x-j)\mathrm{d}x={\delta }_{\mathit{ij}}.\hfill \end{array}$$(28)

A common notation for these dilation/translation combinations is to write

$$\begin{array}{c}\hfill {\varphi }_{\mathit{si}}(x)=2^{s/2}\varphi (2^{s}x-i)\text{and}{\psi }_{\mathit{si}}(x)=2^{s/2}\psi (2^{s}x-i),\end{array}$$(29)

so that

$$\begin{array}{cc}& {\displaystyle {\int }_0^1}{\varphi }_{\mathit{si}}(x){\psi }_{\mathit{sj}}(x)\mathrm{d}x=0,\hfill \\ \hfill & {\displaystyle {\int }_0^1}{\varphi }_{\mathit{si}}(x){\varphi }_{\mathit{sj}}(x)\mathrm{d}x={\displaystyle {\int }_0^1}{\psi }_{\mathit{si}}(x){\psi }_{\mathit{sj}}\mathrm{d}x={\delta }_{\mathit{ij}}.\hfill \end{array}$$(30)

One can also see some of these conditions can be extended for comparisons between scales,

$$\begin{array}{cc}& {\displaystyle {\int }_0^1}{\varphi }_{\mathit{si}}(x){\psi }_{s^{\prime }j}(x)=0,\forall s^{\prime }\ge s\text{and}\hfill \\ \hfill & i\in 0,1,\cdots 2^s-1,j\in 0,1,\cdots 2^{s^{\prime }}-1,\hfill \end{array}$$(31)

as well as

$$\begin{array}{c}\hfill {\displaystyle {\int }_0^1}{\psi }_{\mathit{si}}(x){\psi }_{s^{\prime }{i}^{\prime }}(x)={\delta }_{s{s}^{\prime }i{i}^{\prime }}.\end{array}$$(32)

Thus if we pick some base scale s_b the orthogonality conditions (30), (31), and (32) ensure it is possible to write

$$\begin{array}{cc}\hfill f(x)& ={\displaystyle \sum _{i=0}^{2^{s_b}-1}}\langle {\varphi }_{s_{b}i}|f\rangle {\varphi }_{s_{b}i}(x)\hfill \\ \hfill & +{\displaystyle \sum _{s=s_b}^{\infty }}{\displaystyle \sum _{i=0}^{2^s-1}}\langle {\psi }_{\mathit{si}}|f\rangle {\psi }_{\mathit{si}}(x),\hfill \end{array}$$(33)

where

$$\begin{array}{c}\hfill \langle g,f\rangle ={\displaystyle {\int }_0^1}fg\mathrm{d}x,\end{array}$$(34)

for square integrable functions on [0, 1] (Jawerth & Sweldens 1994). It is a celebrated result of Ingrid Daubechies (Daubechies et al. 1993) to demonstrate that there are various classes of functions ϕ and ψ with compact support which satisfy the conditions (30), (31), (32), such that (33) can be used to represent square integrable functions. The specific choice of ϕ and ψ can often be quite important (for discussions on the matter see e.g. Farge et al. 1996; Zhang et al. 2004). The example calculations detailed in this study were performed with various wavelet choices, but no significant differences were observed, so these comparative calculations were omitted for brevity. In what follows all example calculations use the Haar basis.

In practice the series (33) will be truncated at some maximum scale s_m and the most common choice is to have s_b = 0, which prioritises the number of spatial scales used in the expansion, so that we have the following multi-resolution approximation:

$$\begin{array}{c}\hfill f(x)\approx \langle {\varphi }_0|f\rangle {\varphi }_0(x)+{\displaystyle \sum _{s=0}^{s_m}\sum _{i=0}^{2^s-1}}\langle {\psi }_{\mathit{si}},f\rangle {\psi }_{\mathit{si}}(x).\end{array}$$(35)

In what follows we use the equality sign for series such as (35) on the assumption it is understood this is actually an approximation.

The approximation (35) implies the varying behaviour of the function f on scales 2^t with t >  s_m cannot be properly represented. More formally it can be shown that each s_m yields a set of representable functions which are a subset of the square integrable functions, and further that as s_m is increased each previous set of representable functions is a subset of the current set of representable functions obtained by increasing s_m (Jawerth & Sweldens 1994). That is to say increasing s_m will always increase the set of functions representable.

There is an analogy with the truncation of Fourier series which means one cannot represent frequencies below a given scale. This is particularly clear for the Haar wavelet. As with Fourier series truncation, the quality of representation is dependent on the degree of small scale behaviour in the function f. In practice, the target applications in solar physics would involve simulations on a numerical grid. Thus the minimum scale 2^s_m would be determined by the size of the grid on which the magnetic field is resolved, so the approximation would be able to resolve the same scale of function variation as the numerical grid used to resolve the field.

5.2. Three dimensional multi-resolution analysis

In a three-dimensional Cartesian domain V, we expand the function’s behaviour along each direction via a one dimensional multi-resolution expansion. We assume a 3D function H(x) can be written as H_x(x)H_y(y)H_z(z) (Jawerth & Sweldens 1994). By writing each function as a multi-resolution expansion one obtains encounter eight types of combinations (four in 2D) for each scale s:

$$\begin{array}{c}\hfill {\psi }_{\mathit{slmn}}^{\mu }(\mathit{x})=\{\begin{array}{cc}{\varphi }_{\mathit{sl}}(x){\varphi }_{\mathit{sm}}(y){\varphi }_{\mathit{sn}}(z)\hfill & \text{if}\mu =1,\hfill \\ {\psi }_{\mathit{sl}}(x){\varphi }_{\mathit{sm}}(y){\varphi }_{\mathit{sn}}(z)\hfill & \text{if}\mu =2,\hfill \\ {\varphi }_{\mathit{sl}}(x){\psi }_{\mathit{sm}}(y){\varphi }_{\mathit{sn}}(z)\hfill & \text{if}\mu =3,\hfill \\ {\varphi }_{\mathit{sl}}(x){\varphi }_{\mathit{sm}}(y){\psi }_{\mathit{sn}}(z)\hfill & \text{if}\mu =4,\hfill \\ {\psi }_{\mathit{sl}}(x){\varphi }_{\mathit{sm}}(y){\psi }_{\mathit{sn}}(z)\hfill & \text{if}\mu =5,\hfill \\ {\psi }_{\mathit{sl}}(x){\psi }_{\mathit{sm}}(y){\varphi }_{\mathit{sn}}(z)\hfill & \text{if}\mu =6,\hfill \\ {\varphi }_{\mathit{sl}}(x){\psi }_{\mathit{sm}}(y){\psi }_{\mathit{sn}}(z)\hfill & \text{if}\mu =7,\hfill \\ {\psi }_{\mathit{sl}}(x){\psi }_{\mathit{sm}}(y){\psi }_{\mathit{sn}}(z)\hfill & \text{if}\mu =8.\hfill \end{array}\end{array}$$(36)

Writing the respective coefficients as

$$\begin{array}{c}\hfill H_{\mathit{slmn}}^{\mu }={\displaystyle {\int }_{\mathit{V}}}H(\mathit{x}){\psi }_{\mathit{slmn}}^{\mu }(\mathit{x})\mathrm{d}\mathit{x},\end{array}$$(37)

the ensuing multi-resolution decomposition will be

$$\begin{array}{cc}\hfill H(\mathit{x})& =H_0^1{\psi }_0(\mathit{x})\hfill \\ \hfill & +{\displaystyle \sum _{s=0}^{s_m}\sum _{l=0}^{2^s-1}\sum _{m=0}^{2^s-1}\sum _{n=0}^{2^s-1}\sum _{\mu =2}^8}{H}_{\mathit{slmn}}^{\mu }{\psi }_{\mathit{slmn}}^{\mu }(\mathit{x}),\hfill \end{array}$$(38)

see for instance Farge et al. (1996).

5.2.1. Compacted notation

Throughout this study no particular attention was paid to the contributions of individual μ terms, thus for each l, m, n we shall assume the μ summation has been performed. Further, for brevity we define an index k which, when summed over will be assumed to indicate a sum over l, m and n (or just l and n in 2D). Thus we write (38) as

$$\begin{array}{c}\hfill H(\mathit{x})={\displaystyle \sum _{s=0}^{s_m}\sum _{k=0}^{2^s-1}}{H}_{\mathit{sk}}{\psi }_{\mathit{sk}}(\mathit{x}).\end{array}$$(39)

5.2.2. Relative scale contributions

For a function H one can define the total (relative) contribution Q_s(H) to the multi-resolution decomposition at a scale s as

$$\begin{array}{c}\hfill Q_s(H)={\displaystyle \sum _{k=0}^{2^s-1}}{H}_{\mathit{sk}}.\end{array}$$(40)

Similarly we define the relative power P_s(H) at scale s as:

$$\begin{array}{c}\hfill P_s(H)=\frac{{\sum }_{k=0}^{2^s-1}|H_{\mathit{sk}}|}{{\sum }_{s=0}^{s_m}{\sum }_{k=0}^{2^s-1}|H_{\mathit{sk}}|}·\end{array}$$(41)

For comparison, if the function has the required periodicity we can calculate the power contained at each Fourier scale k through the quantity

$$\begin{array}{c}\hfill H_k={\displaystyle \sum _{|\mathit{k}|=k}}{H}_k(\mathit{k}).\end{array}$$(42)

where the H_k are the coefficients of the appropriate Fourier series of H.

5.3. Helicity formulae

We consider multi-resolution expansions (39) for B (one per component) and the multi-resolution expansion of C,

$$\begin{array}{c}\hfill \mathit{C}={\displaystyle \sum _{s=0}^{s_m}\sum _{k=0}^{2^s-1}}{\mathit{C}}_{\mathit{sk}}{\psi }_{\mathit{sk}}(\mathit{x}),\end{array}$$(43)

and substitute them into the helicity integral ∫_V C·B dV. Using the orthogonality relationships (30), (31), and (32) we obtain a summation over the coefficients of the two series

$$\begin{array}{c}\hfill H={\displaystyle \sum _{s=0}^{s_m}\sum _{k=0}^{2^s-1}}{\mathit{C}}_{\mathit{sk}}·{\mathit{B}}_{\mathit{sk}}={\displaystyle \sum _{s=0}^{s_m}\sum _{k=0}^{2^s-1}}{H}_{\mathit{sk}}.\end{array}$$(44)

So H_sk is the helicity contribution at scale s at position k = lmn (summing over all directions μ). As indicated in Fig. 4a the geometrical interpretation of this coefficient corresponds to the winding of the compact domain of scale s, centred at the coordinates indicated by the triplet k, with the field in the z-slice of the volume containing this domain.

5.4. Decompositions of classical magnetic helicity

We wish to stress that the above expression (44) is equally valid for the classical definition of magnetic helicity:

$$\begin{array}{c}\hfill H={\displaystyle \sum _{s=0}^{s_m}\sum _{k=0}^{2^s-1}}{\mathit{A}}_{\mathit{sk}}·{\mathit{B}}_{\mathit{sk}}={\displaystyle \sum _{s=0}^{s_m}\sum _{k=0}^{2^s-1}}{H}_{\mathit{sk}}.\end{array}$$(45)

for a well-defined magnetic vector potential A. Similarly, we can define (following Finn & Antonsen 1985) a spatial decomposition of gauge-invariant relative helicity:

$$\begin{array}{cc}\hfill H_{\mathrm{R}}& ={\displaystyle \sum _{s=0}^{s_m}\sum _{k=0}^{2^s-1}}{(\mathit{A}+{\mathit{A}}_0)}_{\mathit{sk}}·{(\mathit{B}-{\mathit{B}}_0)}_{\mathit{sk}}\hfill \\ \hfill & ={\displaystyle \sum _{s=0}^{s_m}\sum _{k=0}^{2^s-1}}{H}_{\mathrm{R},sk},\hfill \end{array}$$(46)

where A₀ and B₀ are the so-called references fields, whose normal boundary components are the same as those of B and A. The physical intuition of local winding is, however, only retained when using the winding gauge A^w.

6. Multi-resolution analysis of magnetic helicity: illustrative examples

In this section we present two examples which illustrate the benefits of the spatial decomposition offered by a multi-resolution analysis of magnetic helicity. Unless otherwise stated, all quantities in this section have arbitrary units. Both examples are equally valid in either the winding correlation regime introduced here, or by using the classical definition of helicity (with the winding gauge), as the magnetic field B is tangent at the side boundaries.

6.1. Oppositely twisted flux tubes

The first of our illustrative examples is that of a pair of oppositely twisted flux tubes whose vector field takes the form

$$\begin{array}{cc}\hfill \mathit{B}& =\frac{100}{a}exp[\frac{-({(x+0.55)}^2+y^2)}{a^2}](-y,x,0)\hfill \\ \hfill & -\frac{100}{a}exp[\frac{-({(x-0.55)}^2+y^2)}{a^2}](-y,x,0)+(0,0,10),\hfill \end{array}$$(47)

where we have taken a = 0.2. This field (independent of z) is visualised for the domain [−1, 1] × [−1, 1], in Fig. 5. Making an assumption of periodicity (which can be interpreted as an infinite repeating pattern of the form shown here), Fourier analysis indicates that that this magnetic field has zero overall helicity at every scale, even when the sum over k = k is taken with absolute values, as shown in Fig. 6.

Fig. 5. Magnetic field vector plot of Eq. (47) at z = 0, red indicates positive twist, and blue indicates negative twist.

Fig. 6. Fourier decomposition Hk of magnetic helicity for the field (47).

By contrast, in Fig. 7, we plot the H_sk for the wavelet multi-resolution analysis of the magnetic helicity at spatial scales r = 0 → 6 (along with the associated power P_s(H)). The plotting style is that of a “bubblegram”: each sub-domain of helicity given by the multi-resolution analysis is allocated a three-dimensional sphere at its centre. The radius of this sphere is dependent upon the absolute magnitude of the helicity of the sub-domain, and its colour, red or blue, indicating a positive/negative sign respectively. In this case the magnetic field, and thus the magnetic helicity, is independent of z, but we provide the full z range as a visual indicator of the three–dimensional localising nature of wavelets.

Fig. 7. Hsk for s = 0 → 6 applied to (47). The largest scale is a measure of the numerical round-off. At the two smallest scales 2−5, 6, the visual appearance of the bubblegram is distorted by the frequency of data points. The magnetic field and its associated helicity is independent of z.

The bubblegrams indicate that the helicity is well localised in space in accordance with Fig. 5, presenting with the correct sign of twist. It can be seen that the total helicity Q_s(H) at each scale is zero. The absolute magnetic helicity power P_s(H) is well localised in scale, as indicated in Fig. 8. Peak magnetic helicity occurs at half the spatial scale of the domain, which is in good agreement with the distribution of the twist in the magnetic field itself.

Fig. 8. Ps(H) for s = 0 → 6 applied to (47).

6.2. Linked rings

The magnetic helicity associated with two flux tubes, with linking number ℒ, identical individual internal twists 𝒯 and magnetic fluxes Φ is

$$\begin{array}{c}\hfill H_L=2(\mathcal{L}+\mathcal{T}){\mathrm{\Phi }}^2,\end{array}$$(48)

following Berger (1999). A simple example of such linked rings, R₁ and R₂ can be parametrised as

$$\begin{array}{cc}& R_1(r,\theta ,\varphi )=(Rcos(\theta )+rcos(\theta )cos(\varphi ),\hfill \\ \hfill & Rsin(\theta )+rsin(\theta )cos(\varphi ),rsin(\varphi )),\hfill \end{array}$$(49)

and

$$\begin{array}{cc}& R_2(r,\theta ,\varphi )=(C_x,C_y,C_z)+(Rcos(\theta )+rcos(\theta )cos(\varphi ),\hfill \\ \hfill & -rsin(\varphi ),Rsin(\theta )+rsin(\theta )cos(\varphi )),\hfill \end{array}$$(50)

for major radius R, minor radius r ∈ [0, r_m], toroidal angle θ and poloidal angle ϕ. The set C_x, C_y, C_z denote the centre of R₂. An example with r_m = 0.3 and R = 1 is shown in Fig. 9. We define the magnetic fields B_Ri of each ring as the sum of toroidal B_Rit and poloidal B_Rip components, with

$$\begin{array}{cc}& {\mathit{B}}_{{\mathrm{R}}_1,t}(x,y,z)=B_0(-\frac{y}{\sqrt{(x^2+y^2)}},-\frac{x}{\sqrt{q_1}},0),\hfill \\ \hfill & {\mathit{B}}_{{\mathrm{R}}_1,p}(x,y,z)=\mathcal{T}{B}_0(\frac{\mathit{xz}}{r_{\mathit{xy}}},-\frac{\mathit{yz}}{q_1},1-\frac{R}{\sqrt{q_1}}),\hfill \\ \hfill & {\mathit{B}}_{{\mathrm{R}}_2,t}(x,y,z)=B_0(-\frac{z}{\sqrt{q_2}},0,-\frac{x+1}{\sqrt{q_2}}),\hfill \\ \hfill & {\mathit{B}}_{{\mathrm{R}}_2,p}(x,y,z)=\mathcal{T}{B}_0(\frac{(x+1)y}{q_2},-1+\frac{R}{\sqrt{q_2}},\frac{\mathit{yz}}{q_2}),\hfill \end{array}$$(51)

Fig. 9. Linked flux tubes R1 (red) and R2 (green).

where q₁ = (x² + y²) and q₂ = ((x + 1)² + z²).

We choose R = 1 and C_x = 1, C_y = C_z = 0. Such an arrangement has an associated linking number of ℒ = 1, and we assign 𝒯 = −5, B₀ = 7 and r_m = 0.3, giving total magnetic helicity

$$\begin{array}{c}\hfill H_{\mathcal{L}=1,\mathcal{T}=-5}=(2-10){\mathrm{\Phi }}^2=-31.3,\end{array}$$(52)

where Φ = 1.98. In Fig. 10, we plot the magnetic helicity coefficients H_4k. The bubblegram indicates a distribution of magnetic helicity in correspondence to the distribution of the magnetic fields themselves, which we can attribute to the magnetic twist.

Fig. 10. H4k for linked tubes with 𝒯 = −5.

In Fig. 11, we calculate the ratio of the multi-resolution expansion of helicity with that of the analytical result (from Eq. (48)) as a function of scales included, which we define by the measure

$$\begin{array}{c}\hfill N_s(H)=\frac{{\sum }_{s^{\prime }=0}^s{Q}_{s^{\prime }}(H)}{H}\end{array}$$(53)

(with and without internal magnetic twist). Naturally N_s → 1 in the limit s → ∞ as the multi-resolution analysis would be exactly the magnetic helicity. Here we see the scales at which N_s gets significantly close to 1 differ in the two cases. This is as we would expect due to the differing spatial scales between twist (small scale) and linking (large scale). We see in Fig. 12 the regions of compact support for the Haar wavelet’s at scales s = 1, 2. The s = 1 and s = 2 scales tend to cover both tubes to some degree whilst the scales s = 3 and higher generally only cover one tube. This is reflected in Fig. 11 where we see the 𝒯 = 0 field is dominated by scales s = 1, 2, as scales s = 3 and higher will reflect that on the single tube interior scale there is no complex topology. By contrast the 𝒯 = −5 case has a more balanced distribution across the scales.

Fig. 11. Geometrical extent of Hsk at various scales (Haar wavelet). Panels a: s = 1 and b: s = 2 indicate the overlap of the two tubes is contained within the region of compact support. Panel c: (s = 3) the region of compact support will generally only cover one tube.

Fig. 12. Ns(H) calculated for two linked flux tubes. Calculations are shown with either 𝒯 = −5 (blue) or without (𝒯 = 0) internal twist (red).

7. Helicity, energy and topology

We insert the full three-dimensional multi-resolution decomposition of the field B into the correlation function C (15) to obtain,

$$\begin{array}{cc}& \mathit{C}(\mathit{x})=\frac{1}{2\pi }{\displaystyle \sum _{s=0}^{s_m}\sum _{k=0}^{2^s-1}}{\mathit{B}}_{\mathit{sk}}\times {\displaystyle {\int }_{S_z^{\prime }}}\frac{\mathit{r}}{r^2}{\psi }_{\mathit{sk}}(x^{\prime },y^{\prime },z)\mathrm{d}{x}^{\prime }\mathrm{d}{y}^{\prime }\hfill \\ \hfill & \mathit{r}=(x-x^{\prime },y-y^{\prime },0),\hfill \end{array}$$(54)

where the parameter dependence of the wavelet function ψ indicates the integration is over only the in-plane functions of the 3D wavelets. In order to directly compare the helicity to the energy (helicity has an extra dimension of length) we note that if the planes S_z have x and y widths L and aL, respectively then we can write x = uL and y = avL (0 ≤ u, v ≤ 1) so that r_L = L(u − u′,a(v − v′), 0), then

$$\begin{array}{cc}\hfill \mathit{C}(\mathit{x})& =\frac{L}{2\pi }{\displaystyle \sum _{s=0}^{s_m}\sum _{k=0}^{2^s-1}}{\mathit{B}}_{\mathit{sk}}\hfill \\ \hfill & \times {\displaystyle {\int }_{U_z^{\prime }}}\frac{a{\mathit{r}}_L}{r_L^2}{\psi }_{\mathit{sk}}(x(u^{\prime }),y(v^{\prime }),z)\mathrm{d}{u}^{\prime }\mathrm{d}{v}^{\prime },\hfill \end{array}$$(55)

where U_z is the unit square in the x–y plane. Inserting this into (16) we obtain the helicity in terms of the multi-resolution expansion of the field B alone. This calculation is most parsimoniously represented as a quadratic form, so we introduce some notation. We assume a Cartesian domain U_z × [0, h], with U_z a unit plane at height z, then the following quantities are dependent only upon the chosen wavelet, not the magnetic field itself.

$$\begin{array}{cc}\hfill W_i^{s^{\prime }s{k}^{\prime }k}=\frac{L}{2\pi }{\displaystyle {\int }_0^h}& {\displaystyle {\int }_{S_z\times {\displaystyle U_z^{\prime }}}}\frac{a{r}_L^i}{r_L^2}{\psi }_{s^{\prime }{k}^{\prime }}(x(u^{\prime }),y(v^{\prime }),z)\hfill \\ \hfill & {\psi }_{\mathit{sk}}(x,y,z)\mathrm{d}{u}^{\prime }\mathrm{d}{v}^{\prime }\mathrm{d}x\mathrm{d}y\mathrm{d}z,\hfill \end{array}$$(56)

The cross-product in (54) can be represented using a skew-symmetric matrix ${\mathbf{M}}_{\mathit{ij}}^{s^{\prime }s{k}^{\prime }k}$ which takes the form

$$\begin{array}{c}\hfill {\mathbf{M}}_{\mathit{ij}}^{s^{\prime }s{k}^{\prime }k}=\left(\begin{array}{ccc}0& 0& -W_2^{s^{\prime }s{k}^{\prime }k}\\ 0& 0& W_1^{s^{\prime }s{k}^{\prime }k}\\ W_2^{s^{\prime }s{k}^{\prime }k}& -W_1^{s^{\prime }s{k}^{\prime }k}& 0\end{array}\right).\end{array}$$(57)

Then, using the Einstein summation convention we have

$$\begin{array}{c}\hfill H={\displaystyle {\int }_V}\mathit{C}·\mathit{B}\mathrm{d}V={\mathbf{M}}_{\mathit{ij}}^{s{s}^{\prime }k{k}^{\prime }}{B}_{s^{\prime }{k}^{\prime }}^j{B}_{\mathit{sk}}^i.\end{array}$$(58)

We note that in general the wavelet orthogonality relationships cannot be applied to (56) as the in-plane integrals are over different copies of U_z. However, the z integration is over the same domain so $W_i^{s^{\prime }s{k}^{\prime }k}$ will vanish if n′≠n (from the vectors k = l′m′n′ and k = lmn).

7.1. Helicity as a skew symmetric operator

We note that the helicity is being represented as a product of the field at differing positions and scales through a skew-symmetric operator M. This is analogous to the result that the helicity in periodic domains can be represented as the skew symmetric part of the Fourier transform, as discussed in Sect. 2.3. In this case we use the decomposition ${\mathbf{M}}_{\mathit{ij}}^{s^{\prime }s{k}^{\prime }k}=L{\mathbf{I}}_{\mathit{ij}}^{s^{\prime }s{k}^{\prime }k}+{\mathbf{O}}_{\mathit{ij}}^{s^{\prime }s{k}^{\prime }k}$, where ${\mathbf{I}}_{\mathit{ij}}^{s^{\prime }s{k}^{\prime }k}$ (the superscript labelling is for notational convenience) is the identity matrix (one such matrix for each ss′kk′) and

$$\begin{array}{c}\hfill {\mathbf{O}}^{s^{\prime }s{k}^{\prime }k}=\left(\begin{array}{ccc}-L& 0& -W_2^{s^{\prime }s{k}^{\prime }k}\\ 0& -L& W_1^{s^{\prime }s{k}^{\prime }k}\\ W_2^{s^{\prime }s{k}^{\prime }k}& -W_1^{s^{\prime }s{k}^{\prime }k}& -L\end{array}\right),\end{array}$$(59)

so that

$$\begin{array}{c}\hfill H=L{\mathbf{I}}_{\mathit{ij}}^{s^{\prime }s{k}^{\prime }k}{B}_{s^{\prime }{k}^{\prime }}^j{B}_{\mathit{sk}}^i+{\mathbf{O}}_{\mathit{ij}}^{s^{\prime }s{k}^{\prime }k}{B}_{s^{\prime }{k}^{\prime }}^j{B}_{\mathit{sk}}^i.\end{array}$$(60)

The sum of contributions to the first term for which (s′,k′) = (s, k) give the energy of the field (the multi-resolution approximation of the energy to be precise), and we can thus decompose the sum as follows

$$\begin{array}{cc}& H=LE+N,\hfill \end{array}$$(61)

$$\begin{array}{cc}& N=L{\mathit{I}}_{\mathit{ij}}(1-{\delta }^{s^{\prime }s}{\delta }^{k^{\prime }k})B_{s^{\prime }{k}^{\prime }}^j{B}_{\mathit{sk}}^i+{\mathit{O}}_{\mathit{ij}}^{s^{\prime }s{k}^{\prime }k}{B}_{s^{\prime }{k}^{\prime }}^j{B}_{\mathit{sk}}^i,\hfill \end{array}$$(62)

where δ^s′s is the Kronecker delta function. The operator N contains additional topological information which constitutes the helicity. In the limit which the maximum scale parameter s_m (i.e. the smallest spatial scale) tends to ∞ this relationship is exact so there is a linear sum

$$\begin{array}{c}\hfill H(\mathit{B})=LE(\mathit{B})+N(\mathit{B}),\end{array}$$(63)

where N is the multi-resolution representation of a functional of the field which contains the topological information through the quantities $W_i^{s^{\prime }s{k}^{\prime }k}$. A field evolution for which H is conserved requires that

$$\begin{array}{c}\hfill \frac{\mathrm{d}E}{\mathrm{d}t}=-\frac{1}{L}\frac{\mathrm{d}N}{\mathrm{d}t}\end{array}$$(64)

which would (approximately) apply in significantly low plasma β resistive MHD simulations.

7.2. Field line helicity

Using (16) and (55) the fieldline helicity of a field line γ at scale s and position k = lm can be written as

$$\begin{array}{cc}\hfill \mathcal{A}(\gamma )& =\frac{L}{2\pi }{\displaystyle \sum _{s=0}^{s_m}\sum _{k=0}^{2^s-1}{\int }_{\gamma }}\frac{\mathit{B}}{|\mathit{B}|}·{\mathit{B}}_{\mathit{sk}}(z(s))\hfill \\ \hfill & \times {\displaystyle {\int }_{U_z^{\prime }}}\frac{a{\mathit{r}}_L}{r_L^2}{\psi }_{\mathit{sk}}(x(u^{\prime }),y(v^{\prime }),z)\mathrm{d}{u}^{\prime }\mathrm{d}{v}^{\prime }\mathrm{d}s,\hfill \end{array}$$(65)

where the summation over k implies a 2D multi-resolution decomposition, which is why the coefficient B_sk(z) of the multi-resolution expansion has z dependence. The contribution to 𝒜(γ) from one individual scale is then denoted 𝒜_s(γ), which is then further decomposed to individual localities by 𝒜_sk(γ).

Under ideal evolutions 𝒜(γ) is preserved so the sum of 𝒜_s(γ) must be preserved and changes in 𝒜_s(γ) must be balanced across the scales.

7.3. Helicity preserving field evolution

A particular class of fields of significant interest in the solar physics community are braided fields for which B_z >  0,  ∀x ∈ V, and hence all field lines pass through the domain from the bottom to top boundary. In such cases each field lines γ can represented by the points x₀ ∈ S₀ where they are rooted, such that 𝒜_sk(γ)≡𝒜_sk(x₀) and

$$\begin{array}{cc}\hfill H(\mathit{B})& ={\displaystyle {\int }_{S_0}}\mathcal{A}({\mathit{x}}_0)\mathrm{d}x\mathrm{d}y\hfill \\ \hfill & ={\displaystyle \sum _{s=0}^{s_m}\sum _{k=0}^{2^s-1}{\int }_{S_0}}{\mathcal{A}}_{\mathit{sk}}({\mathit{x}}_0)\mathrm{d}x\mathrm{d}y.\hfill \end{array}$$(66)

If the evolution is not ideal but such that the helicity is conserved (low plasma β resistive relaxations) the distribution 𝒜(γ) changes but the summation (66) must be preserved. In particular we have an alternative means of calculating the value of the operator N(B),

$$\begin{array}{cc}\hfill N(\mathit{B})& =H(\mathit{B})-LE(\mathit{B})\hfill \\ \hfill & ={\displaystyle \sum _{s=0}^{s_m}\sum _{k=0}^{2^s-1}{\int }_{S_0}}{\mathcal{A}}_{\mathit{sk}}({\mathit{x}}_0)\mathrm{d}x\mathrm{d}y-LE.\hfill \end{array}$$(67)

The advantage is that the field line helicity representation of N is linear in both s and k so, for example, we can decompose the contributions to N as the difference H − LE at each scale s, and this decomposition is orthogonal. It is this form which we choose to utilise in this study.

8. Fieldline helicity example: analytical magnetic reconnection (the Dundee braid)

Following the resistive MHD based braiding experiments in Wilmot-Smith et al. (2009, 2011) and Russell et al. (2015) we define a field composed of exponential twist units B_t(b₀, k, a, l, x_c, y_c, z_c) given by

$$\begin{array}{cc}\hfill {\mathit{B}}_t(b_0,k,a,l,x_c,y_c,z_c)& =\frac{2b_{0}k}{a}\exp (-\frac{{(x-x_c)}^2+{(y-y_c)}^2}{a^2}\hfill \\ \hfill & -\frac{{(z-z_c)}^2}{l^2})\mathit{R},\hfill \end{array}$$(68)

where

$$\begin{array}{c}\hfill \mathit{R}=(-(y-y_c),x-x_c,0).\end{array}$$(69)

The parameter b₀ 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_c, y_c, z_c). The braided field is then defined as a superposition of n pairs of positive and negative twists and a uniform vertical background field

$$\begin{array}{cc}\hfill {\mathit{B}}_b(1,a,l,d,z_0,s_d,n)& ={\displaystyle \sum _{i=1}^N}[{\mathit{B}}_t(1,1,a,l,-d,0,z_0+s_d(i-1))\hfill \\ \hfill & +{\mathit{B}}_t(1,-1,a,l,0,d,0,z_0+s_d(i-1))]+\widehat{z},\hfill \end{array}$$(70)

where d is the offset from the central axis, and s_d is the vertical spacing between consecutive twists (of the same sign) and z₀ the height of the first twist unit. In the literature, this is sometimes referred to as a “Dundee Braid”. By altering the extent of the twist units (the parameters a and l) one can control the overlap. The field lines in the region of overlap show significant entanglement (Fig. 13a) a property very well captured by the field line helicity distribution 𝒜(γ) (Fig. 13b). The helicity of this field is (with a suitable choice of parameters) essentially zero owing to the balance of positive and negative twisting. It was found that under a high magnetic Reynold’s number resistive MHD relaxation, under which the helicity is approximately conserved (Wilmot-Smith et al. 2011; Russell et al. 2015), that the field was able to simplify via localised reconnection into (roughly) a pair of oppositely twisted flux ropes.

Fig. 13. Topological measures of (70). Panel a: subset of the field lines in the region where the fields opposing twist units overlap. The field line helicity of the green field line indicated would have contributions due to its own complex geometry as well as its entanglement with the field. Panel b: field line helicity distribution (calculated using the code used in Prior & Yeates 2018) of (71) with t = 0, there is significant small scale structure indicating the field’s complex entanglement.

To keep matters simple in this first application of the multi-resolution decomposition 𝒜_sk, we define a rough analytic approximation of this relaxation process with the following parameterised magnetic field:

$$\begin{array}{c}\hfill \mathit{B}={\mathit{B}}_b(1,D_1(t),D_2(t),1,-20,8,3)\end{array}$$(71)

where

$$\begin{array}{cc}& D_1(t)=\sqrt{2(1-t)},\hfill \\ \hfill & D_2(t)=2(1+2t).\hfill \end{array}$$

This field is considered in a domain x, y ∈ [−4, 4], z ∈ [−24, 24], which are the dimensions (and parameters for t = 0) used in Wilmot-Smith et al. (2009, 2011) and Russell et al. (2015). As t increases the twisted units become more and more separated in the horizontal direction, as shown in Fig. 14. The twist units (with the same sign) also merge vertically to form two non overlapping twisted flux tubes at t = 1. The decrease in overlap between the oppositely twisted units also tends to reduce the complex field entanglement (as we shall shortly see this is not true for low t). It was checked numerically that the total helicity H(B, x, y, z, t) (essentially) remains zero for all t, a property designed to approximate the numerically observed conservation of helicity in the low plasma β MHD simulations. The Fourier expansion of the magnetic helicity of this field is zero throughout (even when an absolute magnitude sum is used). In Fig. 15, we present the field line distribution of the scale decomposed field line helicity decomposition 𝒜_s(x₀), which involves a summation over the spatial parameter k. We remind the reader that for the field line helicity there is one such summation for each point x₀ (i.e. each field line) – hence this is still a spatial distribution. The evolution of these distributions is shown at times t = 0, 0.2, 0.4, 0.6, 0.8, 0.95.

Fig. 14. Vector plots of (71). From top left: four time steps t = 0, 0.3, 0.6, 0.9 at z = 0. Red (blue) denotes the positively (negatively) twisted regions.

Fig. 15. Field line helicity distributions 𝒜s(x0) of (71). The columns from left to right represent time steps t = 0, 0.2, 0.4, 0.6, 0.8, 0.95 integrated on a domain [−4, 4]2 in x, y and [−24, 24] in z, with 400 × 400 field lines. The vertical direction indicates increasing scale s (except for the top line which is the total sum over s (it is the actual distribution 𝒜)).

A couple of observations are worth making. Firstly, at t = 0 all scales 𝒜_s(x₀) show (to varying degrees) the complex mixing pattern present in the full distribution. This is a result of the field line geometry (i.e. the geometry of the green curve in Fig. 13a). Eventually this pattern disappears as the field lines reconnect and disentangle, again this is true of all scales. Secondly, there is a surrounding distribution which is most clear at the scales s = 1, 2; this persists throughout the relaxation. This is the twisted field structure of the field itself; as indicated in the twisted tube example of Sect. 6.1 twisted tube structures (which always compose the field in some manner) are dominated by contributions at these scales. Over the whole sum (over s at each t) these contributions cancel.

To quantify the entanglement variation highlighted in the first point we define a mixing parameter M as

$$\begin{array}{c}\hfill M=(|\frac{\partial {\mathcal{A}}_s({\mathit{x}}_0)}{\partial x}{|}^2+|\frac{\partial {\mathcal{A}}_s({\mathit{x}}_0)}{\partial y}{|}^2{)}^{\frac{1}{2}},\end{array}$$(72)

which will highlight the regions in which we see a rapid change in sign between positive and negative field line helicity 𝒜_s(x₀). Admittedly this will also capture simpler radial decay, but such contributions should be sufficiently weaker. The mixing associated with each scale, in the style of Fig. 15, is shown in Fig. 16. There are two observations. First that the mixing actually increases at first up to t = 0.4 then it decays. Second that the decay is more pronounced at larger length scales (smaller s).

Fig. 16. Mixing distributions M corresponding to Fig. 15. The distributions are shown at time steps t = 0, 0.2, 0.4, 0.6, 0.8, 0.95 on a domain [−4, 4]2 in x, y and [−24, 24] in z, with 400 × 400 field lines. The vertical direction indicates increasing scale s (except for the top line which is the total sum over s (it is the actual distribution 𝒜)).

In Fig. 17 we plot the total signed contribution per scale

$$\begin{array}{c}\hfill Q_S(\mathcal{A})={\displaystyle {\int }_{S_0}}{\mathcal{A}}_s({\mathit{x}}_0)\mathrm{d}{x}_0\mathrm{d}{y}_0,\end{array}$$(73)

spatial integrals over the distributions shown in Fig. 14. Note that we use the Q_s notation used to indicate spatial summation earlier, here it includes the spatial integration over all fieldlines, that is all x₀. There is always (approximately) as much negative as positive contribution, reflecting the total helicity conservation of the field. These values are dominated by the lower scale. Their relative magnitudes increase up to about t = 0.4 then decrease over time. It is interesting that the balance of positive and negative values is always maintained by the same scales (albeit with decreasing magnitudes).

Fig. 17. Total fieldline helicity power QS(𝒜) attributed to each spatial scale. The distribution for s ∈ [0, 6] is shown for time periods t = 0 to t = 0.95 for the Dundee braid relaxation.

In Fig. 18 we plot the absolute power P_s(𝒜):

$$\begin{array}{c}\hfill P_S(\mathcal{A})=\frac{{\int }_{S_0}|{\mathcal{A}}_s({\mathit{x}}_0)|\mathrm{d}{x}_0\mathrm{d}{y}_0}{{\sum }_{s=0}^{s_m}{\int }_{S_0}|{\mathcal{A}}_s({\mathit{x}}_0)|\mathrm{d}{x}_0\mathrm{d}{y}_0},\end{array}$$(74)

associated with each spatial scale for time steps t = 0 to t = 0.95. For early times the values (mostly) decrease with s. However, as the twisting units separate and merge the scale s = 2 becomes more prevalent, reflecting the coherent development of the twisted flux ropes. In Figs. 19 and 20 we see the total power normalised power across all scales of both the field line helicity 𝒜 and the mixing M as a function of time, given by

$$\begin{array}{c}\hfill P_T(H)=\frac{P_s(H)}{{\max }_t{P}_s(H)}·\end{array}$$(75)

Qualitatively the plots are very similar, showing a peak around 0.35 and then a relatively large drop as the twist units properly separate.

We can also directly compare the evolution of the scaled-fieldline helicity with that of magnetic energy, as shown in Fig. 21, where we plot the absolute normalised magnetic energy against that of fieldline helicity, normalised within each scale (P_T(𝒜) versus P_T(E)).

Fig. 18. Total fieldline helicity absolute normalised power PS(𝒜) attributed to each spatial scale. The distribution for s ∈ [0, 6] is shown for time periods t = 0 to t = 0.95 for the Dundee braid relaxation.

Fig. 19. Temporal evolution of PT(𝒜). Calculations are performed from t = 0 to t = 0.95 for the Dundee braid relaxation.

Fig. 20. Temporal evolution of PT(M). Calculations are performed from t = 0 to t = 0.95 for the Dundee braid relaxation.

Fig. 21. Temporal evolution of PT for 𝒜 (field line helicity) and E (magnetic energy). Each panel represents a scale s (except the first with is the total over all scales). The normalisation is scale dependent.

The correlation between the two time series is seen to decrease as the spatial scale decreases in size. Their relative decay is most strongly aligned at scales 2⁰ − 2⁻². Whilst the decay associated with fieldline helicity power is fairly consistent at all scales, the decay of magnetic energy is opposite to that of field line helicity at scales 2⁻⁵ and 2⁻⁶. It is no surprise that the scales s = 0, 2 are the most aligned. As we see in Fig. 17 these are the dominant contributors to the field line helicity variations in the field. As the magnitude of these peaks rise (up to t = 0.3) and fall t >  0.3 (Fig. 17) so concurrently does the energy. This is a potentially important observation: that the variations in the multi-resolution decomposition of the field line helicity 𝒜_sk are intimately correlated with the variations in energy in the field. In future studies it would be interesting to see whether this correlation is maintained in resistive relaxation simulations.

9. Flux of magnetic helicity

The flux of magnetic helicity through a surface is typically defined by

$$\begin{array}{cc}\hfill \frac{\mathrm{d}H}{\mathrm{d}t}& =-2{\displaystyle {\int }_V}\mathit{E}·\mathit{B}{\mathrm{d}}^3\mathrm{x}\hfill \\ \hfill & +{\displaystyle {\int }_S}(({\mathit{A}}_0·\mathit{v})\mathit{B}+({\mathit{A}}_0·\mathit{B})\mathit{v})·\widehat{\mathit{n}}{\mathrm{d}}^2\mathrm{x},\hfill \end{array}$$(76)

for a reference field A₀ (as in relative helicity) uniquely defined by the appropriate boundary conditions of magnetic field B, and velocity field v. The first term refers to dissipation within the volume, which has been shown to have an effective time scale less than energy dissipation, and we thus disregard it. The second expression can be interpreted as the sum of two individual fluxes: the effect of twisting motions on the boundary, and secondly the movement of magnetic field through the boundary.

Wavelet analysis allows us to define a fourth measure of helicity flux, giving an indication of how helicity moves spatially within the volume. An intuitive example of this could be a study of a coronal loop expanding through a simulated region, for which the twist associated with the flux rope would be seen to move spatially. multi-resolution analysis measures helicity as a set of coefficients H_sk attributed to a given scale and spatial domain (with compact support). We can then simply define

$$\begin{array}{c}\hfill \frac{\mathrm{d}{H}_{t,sk}}{\mathrm{d}t}=\frac{H_{t,sk}-H_{t-\delta t,sk}}{\delta t},\end{array}$$(77)

in the form of a finite difference approximation, for the multi-resolution analyses of two adjacent time snapshots.

Further, we can perform a direct (2D) multi-resolution analysis on each term of the analytical measure of flux. For instance,

$$\begin{array}{c}\hfill \frac{\mathrm{d}H}{\mathrm{d}{t}_{\mathit{sk}}}={\displaystyle {\int }_S}({\mathit{A}}_0·\mathit{v}){\psi }_{\mathit{sk}}(\mathit{x}){\mathrm{d}}^{2}x·{\displaystyle {\int }_S}{B}_r(\mathit{x}){\psi }_{\mathit{sk}}{\mathrm{d}}^{2}x,\end{array}$$(78)