2.1 Parametric curve and coordinate systems

The flux tube is constructed from a given continuous and differentiable closed curve $C$ without vanishing curvature in three-dimensional Euclidean space $\unicode[STIX]{x211D}^{3}$. The parametric equation of $C$ is $\boldsymbol{x}=\boldsymbol{c}(s)$, where $\boldsymbol{x}=(x,y,z)$ denotes spatial Cartesian coordinates and $s\in [0,L_{C})$ is the arclength parameter with the length $L_{C}$ of $C$. The Frenet–Serret frame on $C$ is

(2.1)$$\begin{eqnarray}\left.\begin{array}{@{}c@{}}{\displaystyle \frac{\text{d}\boldsymbol{T}}{\text{d}s}}=\unicode[STIX]{x1D705}\boldsymbol{N},\\ {\displaystyle \frac{\text{d}\boldsymbol{N}}{\text{d}s}}=-\unicode[STIX]{x1D705}\boldsymbol{T}+\unicode[STIX]{x1D70F}\boldsymbol{B},\\ {\displaystyle \frac{\text{d}\boldsymbol{B}}{\text{d}s}}=-\unicode[STIX]{x1D70F}\boldsymbol{N},\end{array}\right\}\end{eqnarray}$$

where $\boldsymbol{T}\equiv \text{d}\boldsymbol{c}/\text{d}s$ denotes the unit tangent, $\boldsymbol{N}\equiv (\text{d}\boldsymbol{T}/\text{d}s)/|\text{d}\boldsymbol{T}/\text{d}s|$ the unit normal, $\boldsymbol{B}\equiv \boldsymbol{T}\times \boldsymbol{N}$ the unit binormal, $\unicode[STIX]{x1D705}$ the curvature and $\unicode[STIX]{x1D70F}$ the torsion of the curve.

To facilitate the construction of vector field $\boldsymbol{{\mathcal{B}}}$ for the flux tube, we introduce a curved cylindrical coordinate system $(s,\unicode[STIX]{x1D70C},\unicode[STIX]{x1D703})$ in a tubular region surrounding curve $C$ (see Chui & Moffatt Reference Chui and Moffatt1995). The system satisfies

(2.2)$$\begin{eqnarray}\boldsymbol{x}=\boldsymbol{c}(s)+\unicode[STIX]{x1D70C}\cos \unicode[STIX]{x1D703}\,\boldsymbol{N}(s)+\unicode[STIX]{x1D70C}\sin \unicode[STIX]{x1D703}\,\boldsymbol{B}(s).\end{eqnarray}$$

As sketched in figure 1, $\boldsymbol{x}-\boldsymbol{c}(s)$ is perpendicular to $\boldsymbol{T}(s)$, $\unicode[STIX]{x1D70C}$ denotes the radial distance from $\boldsymbol{c}(s)$ and $\unicode[STIX]{x1D703}$ denotes the azimuthal angle from $\boldsymbol{N}(s)$ in the plane $S_{C}$ spanned by $\boldsymbol{N}(s)$ and $\boldsymbol{B}(s)$. We specify $\boldsymbol{{\mathcal{B}}}$ in terms of $(s,\unicode[STIX]{x1D70C},\unicode[STIX]{x1D703})$ in the tubular region,

(2.3)$$\begin{eqnarray}\unicode[STIX]{x1D6FA}_{C}=\left\{\boldsymbol{x}\,\Big|\,\boldsymbol{x}\in \unicode[STIX]{x211D}^{3},~\min _{s}|\boldsymbol{x}-\boldsymbol{c}(s)|<R_{v}\right\},\end{eqnarray}$$

where curve $C$ is the central axis, and $R_{v}$ denotes the radius of $\unicode[STIX]{x1D6FA}_{C}$ and is set to be small enough so that $s$ is uniquely determined by (2.2) (see Xiong & Yang Reference Xiong and Yang2019a).

Figure 1.

Schematic diagram of the relation between coordinates $(s,\unicode[STIX]{x1D70C},\unicode[STIX]{x1D703})$ and $(x,y,z)$. The flux tube is sketched by dashed curves, and its central axis is marked by the blue solid curve.

Within $\unicode[STIX]{x1D6FA}_{C}$, the coordinate systems $\boldsymbol{x}$ and $(s,\unicode[STIX]{x1D70C},\unicode[STIX]{x1D703})$ are convertible except for the $\unicode[STIX]{x1D703}$-component on $C$. Xiong & Yang (Reference Xiong and Yang2019a) proved that $(s,\unicode[STIX]{x1D70C})$ determined by (2.2) in $\unicode[STIX]{x1D6FA}_{C}$ satisfies

(2.4)$$\begin{eqnarray}\left.\begin{array}{@{}c@{}}s=\displaystyle \text{arg min}_{\unicode[STIX]{x1D701}}|\boldsymbol{x}-\boldsymbol{c}(\unicode[STIX]{x1D701})|,\\ \unicode[STIX]{x1D70C}=\displaystyle \text{min}_{\unicode[STIX]{x1D701}}|\boldsymbol{x}-\boldsymbol{c}(\unicode[STIX]{x1D701})|,\end{array}\right\}\end{eqnarray}$$

and the azimuthal coordinate can be expressed as

(2.5)$$\begin{eqnarray}\unicode[STIX]{x1D703}=\arccos \frac{[\boldsymbol{x}-\boldsymbol{c}(s)]\boldsymbol{\cdot }\boldsymbol{N}(s)}{\unicode[STIX]{x1D70C}}.\end{eqnarray}$$

2.2 Knotted vector field with tunable twist

We develop a general method for constructing a solenoidal vector field

(2.6)$$\begin{eqnarray}\boldsymbol{{\mathcal{B}}}=b^{s}(s,\unicode[STIX]{x1D70C},\unicode[STIX]{x1D703})\boldsymbol{e}_{s}+b^{\unicode[STIX]{x1D70C}}(s,\unicode[STIX]{x1D70C},\unicode[STIX]{x1D703})\boldsymbol{e}_{\unicode[STIX]{x1D70C}}+b^{\unicode[STIX]{x1D703}}(s,\unicode[STIX]{x1D70C},\unicode[STIX]{x1D703})\boldsymbol{e}_{\unicode[STIX]{x1D703}}\end{eqnarray}$$

in $\unicode[STIX]{x1D6FA}_{C}$ on the local frame spanned by unit vectors

(2.7)$$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\boldsymbol{e}_{s}=\boldsymbol{T},\\ \boldsymbol{e}_{\unicode[STIX]{x1D70C}}=\cos \unicode[STIX]{x1D703}\boldsymbol{N}+\sin \unicode[STIX]{x1D703}\boldsymbol{B},\\ \boldsymbol{e}_{\unicode[STIX]{x1D703}}=-\sin \unicode[STIX]{x1D703}\boldsymbol{N}+\cos \unicode[STIX]{x1D703}\boldsymbol{B},\end{array}\right\}\end{eqnarray}$$

for the knotted flux tube with finite thickness and tunable twist, where $b^{s}$, $b^{\unicode[STIX]{x1D70C}}$ and $b^{\unicode[STIX]{x1D703}}$ are the components of $\boldsymbol{{\mathcal{B}}}$ in the axial, radial and azimuthal directions, respectively. In particular, the component $b^{\unicode[STIX]{x1D703}}$ characterizes the helix of field lines around the central axis. This method is an extension of the constructions for straight helical tubes in Gold & Hoyle (Reference Gold and Hoyle1960) and for knotted tubes without twist in Xiong & Yang (Reference Xiong and Yang2019a).

Considering the geometry of knotted flux tubes, we propose three constraints as follows:

1. (i) The axial component of $\boldsymbol{{\mathcal{B}}}$ depends only on $\unicode[STIX]{x1D70C}$.

2. (ii) All the points of a field line of $\boldsymbol{{\mathcal{B}}}$ have the same distance from the central axis $C$.

3. (iii) On a field line of $\boldsymbol{{\mathcal{B}}}$, $\unicode[STIX]{x1D703}$ varies with $s$, and its rate of change depends only on $s$.

From constraint (i), we express the axial component of $\boldsymbol{{\mathcal{B}}}$ in (2.6) as

(2.8)$$\begin{eqnarray}b^{s}=\unicode[STIX]{x1D6E4}f(\unicode[STIX]{x1D70C}),\end{eqnarray}$$

where the constant $\unicode[STIX]{x1D6E4}$ denotes the strength of the toroidal flux along axis $C$, and the compactly supported kernel function $f(x)$ with $x\in [0,\infty )$ satisfies

(2.9)$$\begin{eqnarray}\left.\begin{array}{@{}c@{}}f(\unicode[STIX]{x1D706})=0,\quad \unicode[STIX]{x1D706}\geqslant R_{v},\\ 2\unicode[STIX]{x03C0}\displaystyle \int _{0}^{R_{v}}f(\unicode[STIX]{x1D706})\unicode[STIX]{x1D706}\,\text{d}\unicode[STIX]{x1D706}=1.\end{array}\right\}\end{eqnarray}$$

From constraint (ii), an arbitrary field line must be on a surface

(2.10)$$\begin{eqnarray}\unicode[STIX]{x1D70C}=C_{\unicode[STIX]{x1D70C}},\end{eqnarray}$$

with constant distance $C_{\unicode[STIX]{x1D70C}}$ from the central axis $C$.

From constraint (iii), the azimuth on a field line changes with $s$ at rate $\unicode[STIX]{x1D709}(s)$. Thus the number of turns of the field line around $C$ in the interval $[0,s]$ is

(2.11)$$\begin{eqnarray}\unicode[STIX]{x1D6EF}(s)=\frac{1}{2\unicode[STIX]{x03C0}}\int _{0}^{s}\unicode[STIX]{x1D709}(\unicode[STIX]{x1D706})\,\text{d}\unicode[STIX]{x1D706}.\end{eqnarray}$$

Namely, an arbitrary field line is on the surface

(2.12)$$\begin{eqnarray}\unicode[STIX]{x1D703}=2\unicode[STIX]{x03C0}\unicode[STIX]{x1D6EF}(s)+\unicode[STIX]{x1D703}_{0},\end{eqnarray}$$

where the constant $\unicode[STIX]{x1D703}_{0}$ represents the azimuth of the field line at $s=0$. As illustrated in figure 2, the field line passing through $(s,\unicode[STIX]{x1D70C},\unicode[STIX]{x1D703})=(0,C_{\unicode[STIX]{x1D70C}},\unicode[STIX]{x1D703}_{0})$ is on the intersection between surfaces (2.10) and (2.12).

Figure 2.

The field line can be expressed as the intersection between isosurfaces of (2.10) and (2.12).

Furthermore, the topology of twisted field lines is related to $\unicode[STIX]{x1D6EF}(L_{C})$. As illustrated in figure 3, the field line forms a closed curve on the isosurface of $\unicode[STIX]{x1D70C}$ for rational $\unicode[STIX]{x1D6EF}(L_{C})$, and it is densely distributed for irrational $\unicode[STIX]{x1D6EF}(L_{C})$.

Figure 3.

Field lines on the isosurface of $\unicode[STIX]{x1D70C}$ with (a) $\unicode[STIX]{x1D6EF}(s)=5s/L_{C}$ and (b) $\unicode[STIX]{x1D6EF}(s)=(2+\unicode[STIX]{x03C0})s/L_{C}$.

Substituting (2.10) and (2.12) into (2.2) yields the parametric equation for field lines of $\boldsymbol{{\mathcal{B}}}$ as

(2.13)$$\begin{eqnarray}\boldsymbol{x}(s)=\boldsymbol{c}(s)+C_{\unicode[STIX]{x1D70C}}\cos [2\unicode[STIX]{x03C0}\unicode[STIX]{x1D6EF}(s)+\unicode[STIX]{x1D703}_{0}]\boldsymbol{N}(s)+C_{\unicode[STIX]{x1D70C}}\sin [2\unicode[STIX]{x03C0}\unicode[STIX]{x1D6EF}(s)+\unicode[STIX]{x1D703}_{0}]\boldsymbol{B}(s),\end{eqnarray}$$

and the direction of $\boldsymbol{{\mathcal{B}}}$ is

(2.14)$$\begin{eqnarray}\frac{\text{d}\boldsymbol{x}}{\text{d}s}=\frac{\text{d}\boldsymbol{c}}{\text{d}s}+C_{\unicode[STIX]{x1D70C}}\cos \unicode[STIX]{x1D703}\left(\frac{\text{d}\boldsymbol{N}}{\text{d}s}+2\unicode[STIX]{x03C0}\frac{\text{d}\unicode[STIX]{x1D6EF}}{\text{d}s}\boldsymbol{B}\right)+C_{\unicode[STIX]{x1D70C}}\sin \unicode[STIX]{x1D703}\left(\frac{\text{d}\boldsymbol{B}}{\text{d}s}-2\unicode[STIX]{x03C0}\frac{\text{d}\unicode[STIX]{x1D6EF}}{\text{d}s}\boldsymbol{N}\right).\end{eqnarray}$$

Substituting (2.1), (2.10) and (2.11) into (2.14) yields

(2.15)$$\begin{eqnarray}\frac{\text{d}\boldsymbol{x}}{\text{d}s}=(1-\unicode[STIX]{x1D705}\unicode[STIX]{x1D70C}\cos \unicode[STIX]{x1D703})\boldsymbol{e}_{s}+\unicode[STIX]{x1D70C}(\unicode[STIX]{x1D709}+\unicode[STIX]{x1D70F})\boldsymbol{e}_{\unicode[STIX]{x1D703}}.\end{eqnarray}$$

Then using (2.8) and (2.15), we obtain

(2.16)$$\begin{eqnarray}b^{\unicode[STIX]{x1D70C}}=0\end{eqnarray}$$

and

(2.17)$$\begin{eqnarray}b^{\unicode[STIX]{x1D703}}=\unicode[STIX]{x1D6E4}\frac{f\unicode[STIX]{x1D70C}}{1-\unicode[STIX]{x1D705}\unicode[STIX]{x1D70C}\cos \unicode[STIX]{x1D703}}(\unicode[STIX]{x1D709}+\unicode[STIX]{x1D70F}).\end{eqnarray}$$

Finally substituting (2.8), (2.16) and (2.17) into (2.6), we derive the vector field

(2.18)$$\begin{eqnarray}\boldsymbol{{\mathcal{B}}}=\unicode[STIX]{x1D6E4}f(\unicode[STIX]{x1D70C})\left[\boldsymbol{e}_{s}+\frac{\unicode[STIX]{x1D70C}}{1-\unicode[STIX]{x1D705}\unicode[STIX]{x1D70C}\cos \unicode[STIX]{x1D703}}(\unicode[STIX]{x1D709}+\unicode[STIX]{x1D70F})\boldsymbol{e}_{\unicode[STIX]{x1D703}}\right]\end{eqnarray}$$

for coiled and twisted flux tubes. If the tube axis is specified as a knot, $\boldsymbol{{\mathcal{B}}}$ is knotted and is only non-trivial within the multiply connected domain $\unicode[STIX]{x1D6FA}_{C}$ with $\unicode[STIX]{x1D70C}<R_{v}$.

The axial part $\unicode[STIX]{x1D6E4}f\boldsymbol{e}_{s}$ in equation (2.18) corresponds to the vector field with a writhe-dependent helicity (see Xiong & Yang Reference Xiong and Yang2019a), and the azimuthal part$\unicode[STIX]{x1D6E4}f\unicode[STIX]{x1D70C}(\unicode[STIX]{x1D709}+\unicode[STIX]{x1D70F})\boldsymbol{e}_{\unicode[STIX]{x1D703}}/(1-\unicode[STIX]{x1D705}\unicode[STIX]{x1D70C}\cos \unicode[STIX]{x1D703})$ normal to the central axis represents the twist of field lines. The twist can be tuned by varying $\unicode[STIX]{x1D709}$ in (2.18) and characterized by the twist helicity discussed in § 2.3.

We remark that, if the curvature of the central axis is vanishing and the azimuth of field lines in the flux tube changes at a constant rate $\unicode[STIX]{x1D709}(s)=\unicode[STIX]{x1D6FE}$, equation (2.18) degenerates into

(2.19)$$\begin{eqnarray}\boldsymbol{{\mathcal{B}}}=\unicode[STIX]{x1D6E4}f(\unicode[STIX]{x1D70C})(\boldsymbol{e}_{s}+\unicode[STIX]{x1D70C}\unicode[STIX]{x1D6FE}\boldsymbol{e}_{\unicode[STIX]{x1D703}}),\end{eqnarray}$$

which is consistent with the form of $\boldsymbol{{\mathcal{B}}}$ for a straight flux tube in Gold & Hoyle (Reference Gold and Hoyle1960).

As proved below, the vector field constructed by (2.18) is solenoidal, which satisfies the definition of flux tubes for magnetic or vorticity fields.

Proof. From (2.1) and (2.2), the Jacobian matrix

(2.20)$$\begin{eqnarray}\unicode[STIX]{x1D645}\equiv \left(\begin{array}{@{}ccc@{}}{\displaystyle \frac{\unicode[STIX]{x2202}x}{\unicode[STIX]{x2202}s}} & {\displaystyle \frac{\unicode[STIX]{x2202}x}{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}}} & {\displaystyle \frac{\unicode[STIX]{x2202}x}{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}}}\\ {\displaystyle \frac{\unicode[STIX]{x2202}y}{\unicode[STIX]{x2202}s}} & {\displaystyle \frac{\unicode[STIX]{x2202}y}{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}}} & {\displaystyle \frac{\unicode[STIX]{x2202}y}{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}}}\\ {\displaystyle \frac{\unicode[STIX]{x2202}z}{\unicode[STIX]{x2202}s}} & {\displaystyle \frac{\unicode[STIX]{x2202}z}{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}}} & {\displaystyle \frac{\unicode[STIX]{x2202}z}{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}}}\end{array}\right)\end{eqnarray}$$

between coordinate systems $(s,\unicode[STIX]{x1D70C},\unicode[STIX]{x1D703})$ and $(x,y,z)$ is

(2.21)$$\begin{eqnarray}\unicode[STIX]{x1D645}=[(1-\unicode[STIX]{x1D705}\unicode[STIX]{x1D70C}\cos \unicode[STIX]{x1D703})\boldsymbol{T}-\unicode[STIX]{x1D70F}\unicode[STIX]{x1D70C}\sin \unicode[STIX]{x1D703}\,\boldsymbol{N}+\unicode[STIX]{x1D70F}\unicode[STIX]{x1D70C}\cos \unicode[STIX]{x1D703}\,\boldsymbol{B};\cos \unicode[STIX]{x1D703}\,\boldsymbol{N}+\sin \unicode[STIX]{x1D703}\,\boldsymbol{B};-\unicode[STIX]{x1D70C}\sin \unicode[STIX]{x1D703}\,\boldsymbol{N}+\unicode[STIX]{x1D70C}\cos \unicode[STIX]{x1D703}\,\boldsymbol{B}].\end{eqnarray}$$

Applying the inverse function theorem to (2.21), we derive

(2.22)$$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\unicode[STIX]{x1D735}s={\displaystyle \frac{\boldsymbol{T}}{1-\unicode[STIX]{x1D705}\unicode[STIX]{x1D70C}\cos \unicode[STIX]{x1D703}}},\\ \unicode[STIX]{x1D735}\unicode[STIX]{x1D70C}=\cos \unicode[STIX]{x1D703}\,\boldsymbol{N}+\sin \unicode[STIX]{x1D703}\,\boldsymbol{B},\\ \unicode[STIX]{x1D735}\unicode[STIX]{x1D703}={\displaystyle \frac{-\unicode[STIX]{x1D70F}}{1-\unicode[STIX]{x1D705}\unicode[STIX]{x1D70C}\cos \unicode[STIX]{x1D703}}}\boldsymbol{T}+{\displaystyle \frac{1}{\unicode[STIX]{x1D70C}}}(-\sin \unicode[STIX]{x1D703}\,\boldsymbol{N}+\cos \unicode[STIX]{x1D703}\,\boldsymbol{B}).\end{array}\right\}\end{eqnarray}$$

Taking the divergence of (2.6) yields

(2.23)$$\begin{eqnarray}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{{\mathcal{B}}}=\mathop{\sum }_{i=s,\unicode[STIX]{x1D70C},\unicode[STIX]{x1D703}}\unicode[STIX]{x1D735}b^{i}\boldsymbol{\cdot }\boldsymbol{e}_{i}+b^{i}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{e}_{i}.\end{eqnarray}$$

Substituting (2.8), (2.16) and (2.17) into (2.23) yields

(2.24)$$\begin{eqnarray}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{{\mathcal{B}}}=\unicode[STIX]{x1D6E4}\left[\unicode[STIX]{x1D735}f\boldsymbol{\cdot }\boldsymbol{e}_{s}+f\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{e}_{s}+\unicode[STIX]{x1D735}\frac{f\unicode[STIX]{x1D70C}(\unicode[STIX]{x1D709}+\unicode[STIX]{x1D70F})}{1-\unicode[STIX]{x1D705}\unicode[STIX]{x1D70C}\cos \unicode[STIX]{x1D703}}\boldsymbol{\cdot }\boldsymbol{e}_{\unicode[STIX]{x1D703}}+\frac{f\unicode[STIX]{x1D70C}(\unicode[STIX]{x1D709}+\unicode[STIX]{x1D70F})}{1-\unicode[STIX]{x1D705}\unicode[STIX]{x1D70C}\cos \unicode[STIX]{x1D703}}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{e}_{\unicode[STIX]{x1D703}}\right],\end{eqnarray}$$

and using (2.7) and (2.22), we have

(2.25)$$\begin{eqnarray}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{{\mathcal{B}}}=\unicode[STIX]{x1D6E4}\left[\unicode[STIX]{x1D735}\frac{f\unicode[STIX]{x1D70C}(\unicode[STIX]{x1D709}+\unicode[STIX]{x1D70F})}{1-\unicode[STIX]{x1D705}\unicode[STIX]{x1D70C}\cos \unicode[STIX]{x1D703}}\boldsymbol{\cdot }(-\sin \unicode[STIX]{x1D703}\,\boldsymbol{N}+\cos \unicode[STIX]{x1D703}\,\boldsymbol{B})+\frac{\unicode[STIX]{x1D705}f\unicode[STIX]{x1D70C}(\unicode[STIX]{x1D709}+\unicode[STIX]{x1D70F})\sin \unicode[STIX]{x1D703}}{(1-\unicode[STIX]{x1D705}\unicode[STIX]{x1D70C}\cos \unicode[STIX]{x1D703})^{2}}\right].\end{eqnarray}$$

Applying $(-\sin \unicode[STIX]{x1D703}\,\boldsymbol{N}+\cos \unicode[STIX]{x1D703}\,\boldsymbol{B})\boldsymbol{\cdot }\unicode[STIX]{x1D735}s=0$ and $(-\sin \unicode[STIX]{x1D703}\,\boldsymbol{N}+\cos \unicode[STIX]{x1D703}\,\boldsymbol{B})\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D70C}=0$ to (2.25) yields

(2.26)$$\begin{eqnarray}\displaystyle & & \displaystyle \unicode[STIX]{x1D735}\frac{f\unicode[STIX]{x1D70C}(\unicode[STIX]{x1D709}+\unicode[STIX]{x1D70F})}{1-\unicode[STIX]{x1D705}\unicode[STIX]{x1D70C}\cos \unicode[STIX]{x1D703}}\boldsymbol{\cdot }(-\sin \unicode[STIX]{x1D703}\,\boldsymbol{N}+\cos \unicode[STIX]{x1D703}\,\boldsymbol{B})\nonumber\\ \displaystyle & & \displaystyle \quad =\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}}\left(\frac{f\unicode[STIX]{x1D70C}(\unicode[STIX]{x1D709}+\unicode[STIX]{x1D70F})}{1-\unicode[STIX]{x1D705}\unicode[STIX]{x1D70C}\cos \unicode[STIX]{x1D703}}\right)\unicode[STIX]{x1D735}\unicode[STIX]{x1D703}\boldsymbol{\cdot }(-\sin \unicode[STIX]{x1D703}\,\boldsymbol{N}+\cos \unicode[STIX]{x1D703}\,\boldsymbol{B})=-\frac{\unicode[STIX]{x1D705}f\unicode[STIX]{x1D70C}(\unicode[STIX]{x1D709}+\unicode[STIX]{x1D70F})\sin \unicode[STIX]{x1D703}}{(1-\unicode[STIX]{x1D705}\unicode[STIX]{x1D70C}\cos \unicode[STIX]{x1D703})^{2}}.\qquad\end{eqnarray}$$

Finally, substituting (2.26) into (2.25) yields $\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{{\mathcal{B}}}=0$.☐

In the numerical simulation of flux tubes, we usually specify (2.18) in terms of Cartesian coordinates, so we need to transform $\boldsymbol{x}$ into $(s,\unicode[STIX]{x1D70C},\unicode[STIX]{x1D703})$ by (2.4) and (2.5). The numerical implementation of this transformation and the calculation of (2.18) is described in detail in appendix A. We remark that this numerical algorithm is critical for the smoothness of the constructed field and the accuracy of the helicity calculation.

2.3 Helicity

We demonstrate that the helicity of the constructed vector field $\boldsymbol{{\mathcal{B}}}$ in (2.18) has an analytic expression distinguishing contributions from the writhe, total torsion and intrinsic twist (see Moffatt & Ricca Reference Moffatt and Ricca1992). The gauge-invariant helicity is defined as

(2.27)$$\begin{eqnarray}H\equiv \int _{\unicode[STIX]{x1D6FA}}h\,\text{d}\unicode[STIX]{x1D6FA},\end{eqnarray}$$

where $h\equiv \boldsymbol{{\mathcal{A}}}\boldsymbol{\cdot }\boldsymbol{{\mathcal{B}}}$ is the helicity density, with the vector potential $\boldsymbol{{\mathcal{A}}}$ satisfying $\boldsymbol{{\mathcal{B}}}=\unicode[STIX]{x1D735}\times \boldsymbol{{\mathcal{A}}}$, and $\unicode[STIX]{x1D6FA}$ denotes the domain containing the knot. This domain is bounded by a surface with $\boldsymbol{{\mathcal{B}}}\boldsymbol{\cdot }\boldsymbol{n}=0$ (see Berger & Field Reference Berger and Field1984), where $\boldsymbol{n}$ denotes the surface normal.

For the knotted flux tubes consisting of nested toroidal magnetic surfaces, Chui & Moffatt (Reference Chui and Moffatt1995) proposed a feasible approach to derive $H$ based on a ‘zero-framed’ coordinate system

(2.28)$$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\unicode[STIX]{x1D70D}=s,\\ \unicode[STIX]{x1D712}=\unicode[STIX]{x1D719}_{v}(s,\unicode[STIX]{x1D70C},\unicode[STIX]{x1D703}),\\ \unicode[STIX]{x1D719}=\unicode[STIX]{x1D719}_{a}(s,\unicode[STIX]{x1D703}).\end{array}\right\}\end{eqnarray}$$

Here,

(2.29)$$\begin{eqnarray}\unicode[STIX]{x1D719}_{v}(s,\unicode[STIX]{x1D70C},\unicode[STIX]{x1D703})=\unicode[STIX]{x1D70C}/R_{v}\end{eqnarray}$$

is a normalized vector-surface field, whose isosurface is an integral surface of ${\mathcal{B}}$ consisting of field lines. We remark that one of coordinates in the system of Hamada (Reference Hamada1962) or Boozer (Reference Boozer1981) also has this property of $\unicode[STIX]{x1D719}_{v}$. The angular variable $\unicode[STIX]{x1D719}_{a}(s,\unicode[STIX]{x1D703})$ in (2.28) makes the ribbon-like isosurface of $\unicode[STIX]{x1D719}_{a}$ ‘untwisted’. For a smooth curve $C$, we take

(2.30)$$\begin{eqnarray}\unicode[STIX]{x1D719}_{a}(s,\unicode[STIX]{x1D703})=\unicode[STIX]{x1D703}+\frac{2\unicode[STIX]{x03C0}{\mathcal{L}}_{S}}{L_{C}}s,\end{eqnarray}$$

where the self-linking number (see Pohl Reference Pohl1968; Moffatt & Ricca Reference Moffatt and Ricca1992)

(2.31)$$\begin{eqnarray}{\mathcal{L}}_{S}=W_{r}+T_{t}\end{eqnarray}$$

is decomposed into the writhe

(2.32)$$\begin{eqnarray}W_{r}=\frac{1}{4\unicode[STIX]{x03C0}}\oint _{C}\oint _{C}\frac{[\boldsymbol{N}(s)\times \boldsymbol{N}(s^{\prime })]\boldsymbol{\cdot }[\boldsymbol{c}(s)-\boldsymbol{c}(s^{\prime })]}{|\boldsymbol{c}(s)-\boldsymbol{c}(s^{\prime })|^{3}}\,\text{d}s\,\text{d}s^{\prime }\end{eqnarray}$$

and the normalized total torsion

(2.33)$$\begin{eqnarray}T_{t}=\frac{1}{2\unicode[STIX]{x03C0}}\oint _{C}\unicode[STIX]{x1D70F}\,\text{d}s.\end{eqnarray}$$

Substituting (2.18) into the explicit expressions of toroidal and poloidal flux functions $T(x)$ and $P(x)$ in the zero-framed coordinate system (see Kruskal & Kulsrud Reference Kruskal and Kulsrud1958; Berger & Field Reference Berger and Field1984; Chui & Moffatt Reference Chui and Moffatt1995), after some algebra (see Xiong & Yang Reference Xiong and Yang2019a) we have

(2.34)$$\begin{eqnarray}\left.\begin{array}{@{}c@{}}T(x)=F(R_{v}x),\\ P(x)=[{\mathcal{L}}_{S}+\unicode[STIX]{x1D6EF}(L_{C})]F(R_{v}x),\end{array}\right\}\end{eqnarray}$$

with

(2.35)$$\begin{eqnarray}F(\unicode[STIX]{x1D70C})=2\unicode[STIX]{x03C0}\unicode[STIX]{x1D6E4}\int _{0}^{\unicode[STIX]{x1D70C}}f(\unicode[STIX]{x1D706})\unicode[STIX]{x1D706}\,\text{d}\unicode[STIX]{x1D706}.\end{eqnarray}$$

Additionally, the ‘zero-twist’ coordinate system (see Maggioni & Ricca Reference Maggioni and Ricca2009) also can be used to calculate $T(x)$ and $P(x)$.

Finally, substituting (2.34) into the general expression of helicity (see Chui & Moffatt Reference Chui and Moffatt1995)

(2.36)$$\begin{eqnarray}H=2\int _{0}^{1}T(x)\frac{\text{d}P(x)}{\text{d}x}\,\text{d}x\end{eqnarray}$$

in the zero-framed coordinate system, we derive an explicit expression for the helicity decomposition

(2.37)$$\begin{eqnarray}H=\unicode[STIX]{x1D6E4}^{2}(W_{r}+T_{w}).\end{eqnarray}$$

Here, the twist

(2.38)$$\begin{eqnarray}T_{w}=T_{t}+\unicode[STIX]{x1D6EF}(L_{C})\end{eqnarray}$$

is further decomposed into the total torsion $T_{t}$ and the intrinsic twist $\unicode[STIX]{x1D6EF}(L_{C})$. The geometric interpretation of the components $W_{r}$, $T_{t}$ and $\unicode[STIX]{x1D6EF}(L_{C})$ for a closed and self-linked ribbon is illustrated in figure 12 in Moffatt & Ricca (Reference Moffatt and Ricca1992), and these components can be computed from (2.32), (2.33) and (2.11) with given $\boldsymbol{c}(s)$ and $\unicode[STIX]{x1D709}(s)$.

Sometimes terms $\unicode[STIX]{x1D6E4}^{2}W_{r}$ and $\unicode[STIX]{x1D6E4}^{2}T_{w}$ in (2.37) are referred to as the ‘centreline helicity’ and ‘twist helicity’, respectively (e.g. Scheeler et al. Reference Scheeler, Kleckner, Proment, Kindlmann and Irvine2014). In the following numerical simulation, we study the effects of the sum of $T_{t}$ and $\unicode[STIX]{x1D6EF}(L_{C})$, namely the twist helicity, on the dynamical evolution of knotted structures. In addition, $H$ in (2.37) for a closed tube with constant $P(x)/T(x)$ is consistent with an experimentally accessible form $H=\unicode[STIX]{x1D6E4}\oint _{C}\boldsymbol{{\mathcal{A}}}\boldsymbol{\cdot }\text{d}\boldsymbol{l}$ in Bretherton (Reference Bretherton1970) and Scheeler et al. (Reference Scheeler, van Rees, Kedia, Kleckner and Irvine2017), where $\text{d}\boldsymbol{l}$ denotes a line element.

As another theoretical approach, the derivations of (2.18) and (2.37) are re-expressed using Clebsch potentials in appendix B, and then the present method for constructing $\boldsymbol{{\mathcal{B}}}$ is compared with other methods via Clebsch-like potentials in appendix C.

3.1 Numerical construction of trefoil flux tubes

Using the method proposed in § 2, we construct trefoil flux tubes with various $T_{w}$ from the parametric equation for $\boldsymbol{c}(\unicode[STIX]{x1D701})=(c_{x},c_{y},c_{z})$ as

(3.1)$$\begin{eqnarray}\left.\begin{array}{@{}c@{}}c_{x}(\unicode[STIX]{x1D701})=\unicode[STIX]{x03C0}+[1+0.5\cos (3\unicode[STIX]{x1D701})]\cos (2\unicode[STIX]{x1D701}),\\ c_{y}(\unicode[STIX]{x1D701})=\unicode[STIX]{x03C0}+[1+0.5\cos (3\unicode[STIX]{x1D701})]\sin (2\unicode[STIX]{x1D701}),\\ c_{z}(\unicode[STIX]{x1D701})=\unicode[STIX]{x03C0}-0.5\sin (2\unicode[STIX]{x1D701}).\end{array}\right\}\end{eqnarray}$$

Here, $\unicode[STIX]{x1D701}\in [0,2\unicode[STIX]{x03C0})$ is not an arclength parameter, but it does not affect the numerical construction in appendix A.

The writhe $W_{r}=3.52$ and the normalized total torsion $T_{t}=2.48$ of the curve (3.1) are calculated by (2.32) and (2.33) (see Xiong & Yang Reference Xiong and Yang2019a), respectively, and then the self-linking number ${\mathcal{L}}_{S}=6$ is obtained by (2.31).

In the corresponding magnetic field $\boldsymbol{{\mathcal{B}}}$ in (2.18), we set $\unicode[STIX]{x1D6E4}=1$ for the strength of magnetic flux and

(3.2)$$\begin{eqnarray}f(\unicode[STIX]{x1D706})=\frac{1}{2\unicode[STIX]{x03C0}\unicode[STIX]{x1D70E}^{2}}\exp \left(-\frac{\unicode[STIX]{x1D706}^{2}}{2\unicode[STIX]{x1D70E}^{2}}\right)\end{eqnarray}$$

for the flux distribution, with $\unicode[STIX]{x1D70E}=1/(16\sqrt{2\unicode[STIX]{x03C0}})\approx 0.025$. The radius of the flux tube is estimated as $R_{v}\approx 5\unicode[STIX]{x1D70E}$, and the flux tube with this small $R_{v}$ contains over 99.999 % of the toroidal flux in (2.18). We clip the value of $\boldsymbol{{\mathcal{B}}}$ outside the flux tube with $R_{v}$, so that $\boldsymbol{{\mathcal{B}}}$ can be considered as compactly supported.

We construct the trefoil flux tubes in a periodic box of sides $L=2\unicode[STIX]{x03C0}$. This computational domain $\unicode[STIX]{x1D6FA}$ encloses $\unicode[STIX]{x1D6FA}_{C}$, and is discretized on uniform grid points $N^{3}$ with $N=512$. The solenoidal vector potential with the Coulomb gauge is calculated by the Biot–Savart law with the spectral method

(3.3)$$\begin{eqnarray}\boldsymbol{{\mathcal{A}}}={\mathcal{F}}^{-1}\left(\frac{\text{i}\boldsymbol{k}\times \widehat{\boldsymbol{{\mathcal{B}}}}}{|\boldsymbol{k}|^{2}}\right),\end{eqnarray}$$

where ${\mathcal{F}}^{-1}$ denotes the operator of inverse Fourier transform, $\boldsymbol{k}$ the wavenumber in Fourier space, and $\widehat{\boldsymbol{{\mathcal{B}}}}={\mathcal{F}}(\boldsymbol{{\mathcal{B}}})$ the Fourier coefficient of $\boldsymbol{{\mathcal{B}}}$ with the Fourier transform operator ${\mathcal{F}}$.

We set up six cases to study the effect of twist on the evolution of magnetic knots. Since the total torsion $T_{t}$ for the same trefoil knot (3.1) is fixed, $T_{w}$ in (2.38) is tuned by varying the intrinsic twist as $\unicode[STIX]{x1D709}=-\unicode[STIX]{x1D70F}$, 0, $5\times (2\unicode[STIX]{x03C0}/L_{C})$, $10\times (2\unicode[STIX]{x03C0}/L_{C})$, $20\times (2\unicode[STIX]{x03C0}/L_{C})$ and $30\times (2\unicode[STIX]{x03C0}/L_{C})$, with the total length $L_{C}=15.95$ for (3.1). We list $H$, $T_{w}$, $\unicode[STIX]{x1D709}$ and $\unicode[STIX]{x1D6EF}(L_{C})$ of the six trefoil flux tubes in table 1, and we can see that $H$ grows with $T_{w}$ in the present study, as implied from (2.37). The helicity $H$ is calculated by (2.27) with the numerical integration of the product of $\boldsymbol{{\mathcal{A}}}$ computed by (3.3) and $\boldsymbol{{\mathcal{B}}}$ over the periodic box. This helicity is gauge-invariant owing to the vanishing flux of $\boldsymbol{{\mathcal{B}}}$ across the boundary of the integral domain (see Berger & Field Reference Berger and Field1984).

Table 1.

Parameters and relative errors in the numerical construction of six trefoil flux tubes.

It is noted that the effects of the geometry and topology of the central axis of knotted tubes on the helicity have been investigated in Kerr (Reference Kerr2018a) and Xiong & Yang (Reference Xiong and Yang2019a). In particular, the case with $T_{w}=0$ in table 1 corresponds to the knot construction in Xiong & Yang (Reference Xiong and Yang2019a) with $\boldsymbol{{\mathcal{B}}}=\unicode[STIX]{x1D6E4}f(\unicode[STIX]{x1D70C})\boldsymbol{e}_{s}$ and $H=\unicode[STIX]{x1D6E4}^{2}W_{r}$, where the intrinsic twist cancels out the total torsion in (2.38).

We assess the numerical construction and helicity calculation of trefoil flux tubes using (2.37). The relative numerical error in (2.37) is defined as

(3.4)$$\begin{eqnarray}\unicode[STIX]{x1D716}_{H}\equiv \frac{|H-\unicode[STIX]{x1D6E4}^{2}(W_{r}+T_{w})|}{|H|},\end{eqnarray}$$

and the negligible errors in table 1 indicate that the numerical construction is very accurate. In the numerical implementation, $\langle |\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{{\mathcal{B}}}|\rangle \sim O(10^{-7})$ computed from the six knotted fields are negligible, consistent with Theorem 1.

Figure 4 shows isosurfaces of $F(\unicode[STIX]{x1D70C})=0.9$ for three typical trefoil flux tubes in table 1. We observe that the magnetic field is smooth and compact, so the constructed trefoil flux tubes are suitable for initial conditions of DNS. The isosurfaces are colour coded by the helicity density, and all the flux tubes have $h>0$, indicating the same chirality of $\boldsymbol{{\mathcal{B}}}$ and $\boldsymbol{{\mathcal{A}}}$. Some field lines are integrated from points on the isosurfaces, and they are on the tubes implied by (B 5).

Figure 4.

Isosurfaces of $F(\unicode[STIX]{x1D70C})=0.9$ of initial knotted flux tubes with (a) $T_{w}=0$, (b) $T_{w}=12.48$ and (c) $T_{w}=32.48$. The isosurfaces are colour coded by the helicity density. Some field lines are integrated on the surfaces.

Figure 5 plots the spectrum $E_{b}(k)$ of $\boldsymbol{{\mathcal{B}}}$ for the initial trefoil flux tubes. The spectrum with $T_{w}=0$ decays smoothly at high wavenumbers owing to the Gaussian profile $f(\unicode[STIX]{x1D70C})$ in (2.18) (see Melander & Hussain Reference Melander and Hussain1989; Hussain & Duraisamy Reference Hussain and Duraisamy2011; Yao & Hussain Reference Yao and Hussain2020). The spectrum with $T_{w}=32.48$ has slight oscillations at high wavenumbers owing to the small-scale structures generated by helical field lines surrounding the tube axis in figure 4.

Figure 5.

Spectra of $\boldsymbol{{\mathcal{B}}}$ for trefoil flux tubes with different $T_{w}$.

3.2 DNS of MHD flows

The MHD equations for an incompressible flow of constant unit density in terms of the dimensionless velocity $\boldsymbol{u}$ and magnetic induction $\boldsymbol{b}$ read (e.g. Priest & Forbes Reference Priest and Forbes2000)

(3.5)$$\begin{eqnarray}\left.\begin{array}{@{}c@{}}{\displaystyle \frac{\unicode[STIX]{x2202}\boldsymbol{u}}{\unicode[STIX]{x2202}t}}+\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{u}=-\unicode[STIX]{x1D735}p+\boldsymbol{j}\times \boldsymbol{b}+\unicode[STIX]{x1D708}\unicode[STIX]{x1D6FB}^{2}\boldsymbol{u},\\ {\displaystyle \frac{\unicode[STIX]{x2202}\boldsymbol{b}}{\unicode[STIX]{x2202}t}}=\unicode[STIX]{x1D735}\times (\boldsymbol{u}\times \boldsymbol{b})+\unicode[STIX]{x1D708}_{m}\unicode[STIX]{x1D6FB}^{2}\boldsymbol{b},\end{array}\right\}\end{eqnarray}$$

together with $\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}=0$ and $\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{b}=0$, where $p$ denotes the pressure, $t$ the time, $\unicode[STIX]{x1D708}$ the kinematic viscosity, $\unicode[STIX]{x1D708}_{m}$ the magnetic diffusivity and $\boldsymbol{j}=\unicode[STIX]{x1D735}\times \boldsymbol{b}$ the current density. We set $\unicode[STIX]{x1D708}=\unicode[STIX]{x1D708}_{m}=0.01$. The magnetic Reynolds number is $Re_{m}\equiv \unicode[STIX]{x1D6E4}/\unicode[STIX]{x1D708}_{m}=100$, and the magnetic Prandtl number is $Pr_{m}\equiv \unicode[STIX]{x1D708}/\unicode[STIX]{x1D708}_{m}=1$.

The initial fluid is quiescent with

(3.6)$$\begin{eqnarray}\boldsymbol{u}_{0}(\boldsymbol{x})\equiv \boldsymbol{u}(\boldsymbol{x},t=0)=\mathbf{0},\end{eqnarray}$$

and the fluid motion is induced by the Lorentz force in (3.5). The initial magnetic induction is

(3.7)$$\begin{eqnarray}\boldsymbol{b}_{0}\equiv \boldsymbol{b}(\boldsymbol{x},t=0)=\boldsymbol{{\mathcal{B}}}\end{eqnarray}$$

for the trefoil flux tubes constructed in § 3.1.

The DNS of MHD flows is carried out in the periodic box $\unicode[STIX]{x1D6FA}$ discretized on uniform grid points $N^{3}=512^{3}$. A symmetric form of (3.5) with Elsässer variables $\boldsymbol{z}^{\pm }=\boldsymbol{u}\pm \boldsymbol{b}$ (see Elsässer Reference Elsässer1950) is solved using the pseudo-spectral method (see Aluie Reference Aluie2009; Hao, Xiong & Yang Reference Hao, Xiong and Yang2019). Aliasing errors are removed using the two-thirds truncation method with the maximum wavenumber $k_{max}\approx N/3$. The Fourier coefficient of the velocity is advanced in time using a second-order Runge–Kutta scheme, and the time step is chosen to ensure that the Courant–Friedrichs–Lewy number is less than 0.5 for numerical stability and accuracy.

In order to ensure that the present grid can fully resolve the flow evolution, we did a convergence test by conducting another DNS with the same flow parameters on grid points $N^{3}=1024^{3}$, and the statistics presented in § 4.1 are almost identical with those on $N^{3}=512^{3}$ (not shown).

4.1 Flow statistics

We present major flow statistics for the three typical cases $T_{w}=0$, 12.48 and 32.48 listed in table 1. Figure 6 plots the temporal evolution of the magnetic energy $E_{b}=\sum _{\boldsymbol{k}}|\widehat{\boldsymbol{b}}(\boldsymbol{k})|^{2}/2$ and kinetic energy $E_{u}=\sum _{\boldsymbol{k}}|\widehat{\boldsymbol{u}}(\boldsymbol{k})|^{2}/2$, with Fourier coefficients $\widehat{\boldsymbol{b}}={\mathcal{F}}(\boldsymbol{b})$ and $\widehat{\boldsymbol{u}}={\mathcal{F}}(\boldsymbol{u})$. We observe that $E_{b}$ monotonically decays in all the cases, because a part of $E_{b}$ is converted into $E_{u}$ through the magnetic relaxation. At the mean time, $E_{u}$ first increases owing to the energy conversion and peaks around $t=0.03$. Subsequently, the system gradually becomes stable and $E_{u}$ decays at late times owing to the viscous dissipation, indicating the end of the relaxation process.

Figure 6.

The temporal evolution of the (a) magnetic energy and (b) kinetic energy in the DNS of three trefoil flux tubes.

Figure 7 shows the evolution of the magnetic dissipation rate $\unicode[STIX]{x1D700}_{b}=\unicode[STIX]{x1D708}_{m}\sum _{\boldsymbol{k}}(|\boldsymbol{k}||\widehat{\boldsymbol{b}}(\boldsymbol{k})|)^{2}$ and magnetic helicity $H_{b}=\int _{\unicode[STIX]{x1D6FA}}\boldsymbol{a}\boldsymbol{\cdot }\boldsymbol{b}\,\text{d}\unicode[STIX]{x1D6FA}$ with $\boldsymbol{a}={\mathcal{F}}^{-1}(\text{i}\boldsymbol{k}\times \widehat{\boldsymbol{b}}/|\boldsymbol{k}|^{2})$. The maximum magnetic dissipation decreases with $T_{w}$, and the peaks of $\unicode[STIX]{x1D700}_{b}$ generally coincide with the important events of topological changes of magnetic fields, which will be explained in §§ 4.2 and 4.3.

Figure 7.

The temporal evolution of the (a) magnetic dissipation rate and (b) magnetic helicity in the DNS of three trefoil flux tubes.

Figure 8.

Isosurfaces of $|\boldsymbol{b}|$ in the temporal evolution of trefoil flux tubes with $T_{w}=0$, 12.48 and 32.48 (columns from left to right) at $t=0.006$, 0.02, 0.04 and 0.1 (rows from top to bottom). All the isosurfaces are colour coded by the helicity density. The isocontour values of $|\boldsymbol{b}|$ are (a) 65, (b) 65, (c) 60, (d) 35, (e) 26, (f) 65, (g) 30, (h) 16, (i) 50, (j) 13, (k) 10 and (l) 33.

The evolution of $H_{b}$ depends on the conversion between centreline and twist helicities in (2.37), and it can be non-monotonic owing to the cancellation of positive and negative $h$ in the integral of (2.27) (see Laing, Ricca & Sumners Reference Laing, Ricca and Sumners2015; Kerr Reference Kerr2018b). Despite the presence of dissipation effects, $H_{b}$ for the magnetic trefoil knot with $T_{w}=0$ is remarkably conserved, which is similar to the evolution of the trefoil vortex knot in viscous HD flows (see Scheeler et al. Reference Scheeler, Kleckner, Proment, Kindlmann and Irvine2014; Xiong & Yang Reference Xiong and Yang2019a). On the other hand, $H_{b}$ with imposed large $T_{w}=12.48$ and $32.48$ decays rapidly from the beginning.

4.2 Morphology of flux tubes

Next we discuss the relation between the major flow statistics and the morphology of the flux tubes. Under the induction of the Lorentz force, the initially quiescent flux tubes start to evolve, and their motion and morphology are significantly affected by the twisting degree $T_{w}$ of field lines. Figure 8 depicts the evolution of isosurfaces of $|\boldsymbol{b}|$ colour coded by $h$ at $t=0.006$, 0.02, 0.04 and 0.1. In the early stage at $t\leqslant 0.02$, the magnetic knot with $T_{w}=0$ splits into a pair of connected trefoil flux tubes. Similarly, some parts of the knot with $T_{w}=12.48$ split into three offspring tubes, but some parts with high $h$ are still connected. By contrast, the knot with $T_{w}=32.48$ does not split but only shrinks, similar to the relaxation of knotted magnetic filaments in ideal MHD flows (see Moffatt Reference Moffatt1990). Therefore, the twist of field lines with large $T_{w}$ appears to inhibit the splitting of flux tubes.

In addition, although the cases for highly twisted flux tubes (e.g. $T_{w}=32.48$) can exceed the threshold of the kink instability, such as the limit of Hood & Priest (Reference Hood and Priest1979), the evolution of magnetic knots is still stable owing to the numerical accuracy and symmetric configuration in our simulations without initial perturbation modes.

In figures 7(a) and 8(a,b), the occurrence of tube splitting with $T_{w}=0$ and $T_{w}=12.48$ coincides with the peaks of $\unicode[STIX]{x1D700}_{b}$ at $t\approx 0.003$, consistent with the coincidence between the vortex reconnection and the peak of the dissipation rate in the evolution of vortex knots in HD flows (see Xiong & Yang Reference Xiong and Yang2019a). For large $T_{w}=32.48$, $\unicode[STIX]{x1D700}_{b}$ of the knot without tube splitting in figure 8(c) has no obvious peak.

Figure 9 plots plane cuts of the contour of $|\boldsymbol{b}|$ at $x=\unicode[STIX]{x03C0}$. At $t=0$, the contour for $T_{w}=0$ in figure 9(a) is almost identical to those for other $T_{w}$. At $t=0.02$, two cross-sections of each flux tube at $t=0$ tend to approach to each other for all the cases, implying shrinking of the knots. In cases $T_{w}=0$ and 12.48, each cross-section splits into two and three parts with concentrated $|\boldsymbol{b}|$, respectively. In case $T_{w}=32.48$, the cross-section of $|\boldsymbol{b}|$ does not split, consistent with the observation in figure 8.

Figure 9.

Plane cuts of the contour of $|\boldsymbol{b}|$ at $x=\unicode[STIX]{x03C0}$: (a) $T_{w}=0$ at $t=0$; (b) $T_{w}=0$ at $t=0.02$; (c) $T_{w}=12.48$ at $t=0.02$; and (d) $T_{w}=32.48$ at $t=0.02$. The range of the contours is marked by the dashed line in figure 8(d).

At a later time $t=0.04$, separations of the split flux tubes with $T_{w}=0$ and $T_{w}=12.48$ are more notable, whereas the tube with $T_{w}=32.48$ remains the initial topology. Finally, at $t=0.1$, both $E_{b}$ and $E_{u}$ are low in figure 6, so dissipation processes with finite $\unicode[STIX]{x1D708}_{m}$ and $\unicode[STIX]{x1D708}$ dominate. The isosurfaces of $|\boldsymbol{b}|$ for all the cases deviate from flux tubes and their geometry becomes more complex.

We remark that the morphology of the isosurface of $|\boldsymbol{b}|$ can be different from that of field lines (see Xiong & Yang Reference Xiong and Yang2019b), and next we will discuss the detailed dynamics of field lines inside flux tubes.

4.3 Effects of twist on the splitting of flux tubes

We elucidate how the flux tubes split under the Lorentz force $\boldsymbol{f}_{L}\equiv \boldsymbol{j}\times \boldsymbol{b}$. Since $\boldsymbol{u}$ is very small at early times, we rewrite the momentum equation in (3.5) with only the leading-order terms of $\boldsymbol{u}$ as

(4.1)$$\begin{eqnarray}\frac{\unicode[STIX]{x2202}\boldsymbol{u}}{\unicode[STIX]{x2202}t}-\unicode[STIX]{x1D708}\unicode[STIX]{x1D6FB}^{2}\boldsymbol{u}=\boldsymbol{f}_{L}-\unicode[STIX]{x1D735}p.\end{eqnarray}$$

Taking the divergence of (4.1) for the incompressible flow, we have

(4.2)$$\begin{eqnarray}\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\,\boldsymbol{f}_{L}-\unicode[STIX]{x1D735}p)=0.\end{eqnarray}$$

Thus the force exerted on flux tubes can be interpreted as a projection of $\boldsymbol{f}_{L}$ onto divergence-free space.

The Lorentz force after the divergence-free projection can be expressed as

(4.3)$$\begin{eqnarray}\boldsymbol{f}_{p}\equiv \boldsymbol{f}_{L}-\unicode[STIX]{x1D735}\unicode[STIX]{x1D6F7},\end{eqnarray}$$

with $\unicode[STIX]{x1D6F7}=\unicode[STIX]{x0394}^{-1}\boldsymbol{f}_{L}$. In the numerical implementation, $\boldsymbol{f}_{p}$ is calculated by

(4.4)$$\begin{eqnarray}\boldsymbol{f}_{p}=\boldsymbol{f}_{L}-{\mathcal{F}}^{-1}\left(\frac{\widehat{\boldsymbol{f}}_{L}\boldsymbol{\cdot }\boldsymbol{k}}{|\boldsymbol{k}|^{2}}\boldsymbol{k}\right),\end{eqnarray}$$

with $\widehat{\boldsymbol{f}}_{L}={\mathcal{F}}(\,\boldsymbol{f}_{L})$.

Figure 10 plots the projection of integral lines of $\boldsymbol{f}_{p}$ on the plane at $x=\unicode[STIX]{x03C0}$ for the three magnetic knots at $t=0$. These force lines for $T_{w}=0$ form a counter-rotating, dipole-like pattern and drive a bundle of field lines inside the flux tube to split into two, also shown as the splitting tubes in figures 8(d) and 9(b).

Figure 10.

Projection of the integral lines of $\boldsymbol{f}_{p}$ onto the plane at $x=\unicode[STIX]{x03C0}$ for three magnetic knots with (a) $T_{w}=0$, (b) $T_{w}=12.48$ and (c) $T_{w}=32.48$ at $t=0$. The range of each panel is marked by the dotted box in figure 9(a). The grey dotted circle denotes the isosurface of $|\boldsymbol{b}|=5\,\%\,|\boldsymbol{b}|_{max}$, as the boundary of flux tubes. The yellow circle and the dashed curve, respectively, denote the spiral point and the limit cycle in the phase portrait for $\boldsymbol{f}_{p}$.

We find that $T_{w}$ has an impact on the phase portrait for the nonlinear characteristic equation of $\boldsymbol{f}_{L}$. Figure 10(a) for $T_{w}=0$ shows stable and unstable spiral points at the dipole centres on the phase plane at $x=\unicode[STIX]{x03C0}$, and the two critical points have nearly mirror symmetry. For $T_{w}>0$, the stable spiral point becomes a limit cycle in figure 10(b), and this Hopf bifurcation breaks the mirror symmetry of $\boldsymbol{f}_{p}$ in the phase portrait. Furthermore, the area of the limit cycle grows with $T_{w}$ in figure 10(c).

Implied from figures 9 and 10, the limit cycle appears to move faster than the stable spiral point and carry the major potion of high-intensity field lines in the evolution of flux tubes. The structural transition with $T_{w}$ in the phase portrait corresponds to the morphology of splitting/unsplitting flux tubes in figure 8(d–f).

Inside the flux tube, the spiral $\boldsymbol{f}_{L}$ can drive the initially parallel field lines into chaotic ones. For example, figure 11 depicts a highly chaotic and braided field line for the knot with $T_{w}=12.48$ at $t=0.04$, and the chaos of field lines can be further quantified (see Smiet et al. Reference Smiet, Candelaresi, Thompson, Swearngin, Dalhuisen and Bouwmeester2015). On the other hand, the field lines inside the limit cycle are very stable without chaos. Figure 12 depicts the coexistence of an inner invariant torus and an outer chaotic field line. Considering that field lines at high $Re_{m}$ move with the fluid from the Alfvén theorem, the mechanism for the transition to chaos is similar to the chaotic mixing for passive scalars (see Aref Reference Aref1984).

Figure 11.

A chaotic field line colour coded by $|\boldsymbol{b}|$ for the magnetic knot with $T_{w}=12.48$ at $t=0.04$.

Figure 12.

An inner field line (blue) on an invariant torus and an outer field line (red) in the chaotic region for the magnetic knot with $T_{w}=32.48$ at $t=0.04$.

The splitting of flux tubes and the chaotic rearrangement of field lines for $T_{w}=0$ and $T_{w}=12.48$ contribute to the rise of $\unicode[STIX]{x1D700}_{b}$ in figure 7(a) and the relatively low peaks of $E_{u}$ in figure 6(b). Therefore, increasing $T_{w}$ can achieve more stable field lines and less energy loss in the relaxation of flux tubes, which could be useful for the design of tokamaks (see Strauss Reference Strauss1976; Wesson Reference Wesson1976; Wingen et al. Reference Wingen, Wilcox, Seal, Unterberg, Cianciosa, Delgado-Aparicio, Hirshman and Lao1976).

Moreover, in contrast to the significant influence of $T_{w}$ on the evolution of magnetic knots, the twist with varying $\unicode[STIX]{x1D709}$ in table 1 has negligible effects on the evolution of vortex knots in HD flows, which is discussed in appendix D.