2.1 Three-dimensional hierarchical model

In this section, we introduce a 3-D hierarchical model to describe the merger of magnetic flux tubes. Our goal is to obtain explicit expressions for the evolution of quantities such as magnetic energy and the coherence length of magnetic field in both the merging plane and the direction perpendicular to it, as functions of time. Similar to its 2-D counterpart (Z19), the 3-D model is constructed in the framework of resistive incompressible RMHD, although the key ideas should remain valid in some collisionless plasma descriptions.

In our RMHD system, a constant strong background magnetic guide field $\boldsymbol{B}_{g}=B_{g}\hat{\boldsymbol{z}}$ is assumed, and the $xy$ plane is referred to as the perpendicular plane. We consider an ensemble of magnetic flux tubes of the same size and with the same magnetic field strength as the initial condition. The magnetic field in each flux tube can be decomposed into the constant guide field $\boldsymbol{B}_{g}$ and the perpendicular magnetic field $\boldsymbol{B}_{\bot }$, corresponding to an electric current parallel or anti-parallel to the guide field. For simplicity, each flux tube is assumed to be of cylindrical shape, with radius $R$ (in the $xy$ plane) and length $l$ in the $z$ direction. The lengths $R$ and $l$ represent the coherence lengths of the perpendicular magnetic field, $\boldsymbol{B}_{\bot }$, in the directions perpendicular and parallel, respectively, to the guide field. In the RMHD approximation, $B_{\bot }$ and $R$ are asymptotically smaller than $B_{g}$ and $l$ respectively. The magnetic field lines in each flux tube thus have either right- or left-handed helical shapes depending on the direction of the current, with asymptotically small pitch angles. The flux tubes are assumed to be randomly distributed in the system and volume filling, with half of them carrying currents parallel to the guide field, and the other half anti-parallel. Each flux tube is a system of concentric cylindrical magnetic flux surfaces nested around the tube’s single magnetic axis aligned with the guide field. It contains axial magnetic flux $\unicode[STIX]{x1D6F9}_{a}=B_{g}\unicode[STIX]{x03C0}R^{2}$ and poloidal (i.e. azimuthal around the magnetic axis) magnetic flux $\unicode[STIX]{x1D6F9}_{p}=B_{\bot }^{\text{tube}}Rl$, where $B_{\bot }^{\text{tube}}$ is the characteristic perpendicular magnetic field of each flux tube. In the perpendicular $xy$ plane, the magnetic potential of a flux tube is defined as $\unicode[STIX]{x1D713}^{\text{tube}}=B_{\bot }^{\text{tube}}R$, similar to its 2-D counterpart. Thus, $\unicode[STIX]{x1D6F9}_{p}=\unicode[STIX]{x1D713}^{\text{tube}}l$. figure 1 schematically represents the simplified geometry of a single tube as we have just described.

The magnetic energy associated with the perpendicular magnetic field contained in a flux tube is $\unicode[STIX]{x1D716}_{M}\simeq \unicode[STIX]{x03C0}R^{2}l(B_{\bot }^{\text{tube}})^{2}/(8\unicode[STIX]{x03C0})=(\unicode[STIX]{x1D713}^{\text{tube}})^{2}l/8$. Denoting by $N$ the number density of flux tubes per unit volume [$N\sim 1/(\unicode[STIX]{x03C0}R^{2}l)$], the total magnetic energy density (of the perpendicular magnetic field) of the system is, therefore, $E_{M}=\unicode[STIX]{x1D716}_{M}N=(B_{\bot }^{\text{tube}})^{2}/(8\unicode[STIX]{x03C0})$. In addition, we denote the number density of flux tubes per unit area in the $xy$ plane as $N_{\bot }$, where $N_{\bot }\sim 1/(\unicode[STIX]{x03C0}R^{2})$ and $N_{\bot }=Nl$.

A non-zero local magnetic helicity exists associated with each flux tube, defined as $h=\int \text{d}V_{\text{tube}}\,\boldsymbol{A}\boldsymbol{\cdot }\boldsymbol{B}$, where $\boldsymbol{B}=\unicode[STIX]{x1D735}\times \boldsymbol{A}$. In the RMHD limit, however, this expression simplifies as follows:

(2.1)$$\begin{eqnarray}\displaystyle h & = & \displaystyle \int \text{d}V_{\text{tube}}\,\boldsymbol{A}\boldsymbol{\cdot }\boldsymbol{B}=\int \text{d}V_{\text{tube}}(A_{z}B_{g}+A_{\bot }B_{\bot }^{\text{tube}})\nonumber\\ \displaystyle & {\approx} & \displaystyle -B_{g}\int \text{d}V_{\text{tube}}\unicode[STIX]{x1D713}^{\text{tube}}\approx -(B_{g}\unicode[STIX]{x03C0}R^{2})(B_{\bot }^{\text{tube}}Rl)=-\unicode[STIX]{x1D6F9}_{a}\unicode[STIX]{x1D6F9}_{p},\end{eqnarray}$$

represented by the product of axial flux and poloidal flux. In the derivation, the approximation is made based on the RMHD ordering $A_{\bot }\ll A_{z}$ and $B_{\bot }^{\text{tube}}\ll B_{g}$.

When two flux tubes with parallel axial currents are close to each other, they are unstable to the coalescence instability and approach each other (Finn & Kaw Reference Finn and Kaw1977), merging into a wider flux tube via reconnection of their poloidal fluxes. For simplicity, we assume the merging process to happen in discrete stages, denoted by index $n$. We assume the mergers to proceed in hierarchical fashion and, therefore, take all flux tubes to be identical at each stage. The index $n$ is thus added as a subscript to the above quantities to indicate the generation of mergers to which they belong. As the flux tubes merge, the above quantities evolve; i.e. they are functions of $n$ and so, implicitly, of time.

Figure 1.

Schematic of the flux-tube geometry considered in our analytical model. A flux tube is characterized by its radius $R$, length $l$, perpendicular magnetic field $B_{\bot }^{\text{tube}}$, axial flux $\unicode[STIX]{x1D6F9}_{a}$ and poloidal flux $\unicode[STIX]{x1D6F9}_{p}$.

We also assume mergers to be binary ($N_{\bot ,n+1}=N_{\bot ,n}/2$) and governed by basic conservation laws (similar arguments can be found in Medvedev et al. (Reference Medvedev, Fiore, Fonseca, Silva and Mori2005), Polomarov, Kaganovich & Shvets (Reference Polomarov, Kaganovich and Shvets2008), Lyutikov et al. (Reference Lyutikov, Sironi, Komissarov and Porth2017a) and Zrake & Arons (Reference Zrake and Arons2017)).

First, the merger of two flux tubes should conserve the axial magnetic flux: $\unicode[STIX]{x1D6F9}_{a,n+1}=2\unicode[STIX]{x1D6F9}_{a,n}$. Since the axial magnetic field $\boldsymbol{B}_{g}$ is a constant, it follows that the area of the perpendicular cross-section should be conserved and, therefore,

(2.2)$$\begin{eqnarray}R_{n+1}=\sqrt{2}R_{n}.\end{eqnarray}$$

Second, the perpendicular magnetic potential associated with each flux tube is not affected by reconnection and, as in the two-dimensional case (Z19), should remain the same: $\unicode[STIX]{x1D713}_{n+1}^{\text{tube}}=\unicode[STIX]{x1D713}_{n}^{\text{tube}}$. Combined with (2.2), we obtain the evolution of perpendicular magnetic field and its energy density

(2.3)$$\begin{eqnarray}B_{\bot ,n+1}^{\text{tube}}=B_{\bot ,n}^{\text{tube}}/\sqrt{2},\quad E_{M,n+1}=E_{M,n}/2.\end{eqnarray}$$

We assume that the parallel coherence length of each flux tube is determined by the critical balance condition, i.e. the statement that the pertinent parallel and perpendicular time scales should roughly be the same (Goldreich & Sridhar Reference Goldreich and Sridhar1995). Parallel dynamics is dictated by the propagation of shear Alfvén waves along the magnetic field. Therefore, the parallel time scale is $\unicode[STIX]{x1D70F}_{\Vert }\sim l/V_{A}$, with $V_{A}$ the Alfvén speed computed with the guide field. In the perpendicular plane, the relevant time scale is $\unicode[STIX]{x1D70F}_{\bot }\sim R/v_{A}$, with $v_{A}$ the Alfvén speed computed with $B_{\bot }^{\text{tube}}$. That this is the relevant perpendicular time scale to consider remains true even during mergers, since merging tubes are not static, but rather move in the perpendicular plane at roughly the ambient $v_{A}$. Note also that another relevant perpendicular time scale is the reconnection time – that is, the time it takes to merge two tubes through reconnection. This is always slower than Alfvénic, implying that critical balance can always establish itself dynamically.

We thus find

(2.4)$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D70F}_{\Vert }\sim \unicode[STIX]{x1D70F}_{\bot }\Rightarrow l/B_{g}\sim R/B_{\bot }^{\text{tube}}. & & \displaystyle\end{eqnarray}$$

Plugging equations (2.2)–(2.3) into this relation, we obtain the evolution of coherence length of magnetic field in the parallel direction

(2.5)$$\begin{eqnarray}l_{n+1}=2l_{n}.\end{eqnarray}$$

That is, $l\sim R^{2}$, i.e. the magnetic flux tubes grow faster in the parallel direction than on perpendicular planes. Knowing the magnetic potential remains the same during the merging process, the poloidal flux ($\unicode[STIX]{x1D6F9}_{p}=\unicode[STIX]{x1D713}^{\text{tube}}l$) increases with the length of the flux tubes: $\unicode[STIX]{x1D6F9}_{p,n+1}=2\unicode[STIX]{x1D6F9}_{p,n}$. As a result of the increasing flux contained in each generation of flux tubes as they merge and grow, the magnetic helicity associated with each flux tube also increases: $h_{n+1}=4h_{n}$.

Our next step in order to obtain a continuous-time description of the evolution is to compute the lifetime of each generation; this entails a discussion of the reconnection process mediating the mergers.

The Lundquist number corresponding to the merging tubes at any given $n$th stage is defined as $S_{n}\equiv R_{n}v_{A,n}/\unicode[STIX]{x1D702}$, where $\unicode[STIX]{x1D702}$ is the (constant) magnetic diffusivity. Importantly, $S\propto RB_{\bot }^{\text{tube}}/\unicode[STIX]{x1D702}\propto \unicode[STIX]{x1D713}^{\text{tube}}/\unicode[STIX]{x1D702}$. Therefore, and as we pointed out in our earlier work (Z19), the constancy of the magnetic potential $\unicode[STIX]{x1D713}^{\text{tube}}$ during the mergers implies the same for the Lundquist number: $S_{n+1}=S_{n}$. In particular, in MHD, if $S\lesssim 10^{4}$, reconnection proceeds in the Sweet–Parker (SP) regime (Parker Reference Parker1957; Sweet Reference Sweet and Lehnert1958). In this regime, the dimensionless reconnection rate, defined as $\unicode[STIX]{x1D6FD}_{\text{rec}}\equiv v_{\text{rec}}/v_{A}$, scales with the Lundquist number as $\unicode[STIX]{x1D6FD}_{\text{rec}}\simeq S^{-1/2}$. In contrast, if $S\gtrsim 10^{4}$, then magnetic reconnection proceeds in the plasmoid-dominated regime (Biskamp Reference Biskamp1986; Loureiro, Schekochihin & Cowley Reference Loureiro, Schekochihin and Cowley2007; Lapenta Reference Lapenta2008; Bhattacharjee et al. Reference Bhattacharjee, Huang, Yang and Rogers2009; Samtaney et al. Reference Samtaney, Loureiro, Uzdensky, Schekochihin and Cowley2009; Huang & Bhattacharjee Reference Huang and Bhattacharjee2010; Uzdensky, Loureiro & Schekochihin Reference Uzdensky, Loureiro and Schekochihin2010; Loureiro et al. Reference Loureiro, Samtaney, Schekochihin and Uzdensky2012; Loureiro & Uzdensky Reference Loureiro and Uzdensky2016) with $\unicode[STIX]{x1D6FD}_{\text{rec}}\simeq 0.01$. Therefore, we reach the non-trivial conclusion that the normalized reconnection rate, $\unicode[STIX]{x1D6FD}_{\text{rec}}$, remains unchanged throughout the evolution, and the merging process remains in the same reconnection regime (SP or plasmoid dominated) in which it started initially.

We can now proceed to calculate the lifetime of each generation. We note first that, by definition, the time it takes to reconnect (merge) two flux tubes is $\unicode[STIX]{x1D70F}_{\text{rec}}=\unicode[STIX]{x1D6FD}_{\text{rec}}^{-1}R/v_{A}=\unicode[STIX]{x1D6FD}_{\text{rec}}^{-1}\unicode[STIX]{x1D70F}_{\bot }$. Since $\unicode[STIX]{x1D6FD}_{\text{rec}}\ll 1$, the time scale for the $n$th stage is approximately the reconnection time at this stage

(2.6)$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D70F}_{n}\approx \unicode[STIX]{x1D70F}_{\text{rec},\text{n}}\propto \unicode[STIX]{x1D6FD}_{\text{rec}}^{-1}R_{n}/B_{\bot ,n}^{\text{tube}}. & & \displaystyle\end{eqnarray}$$

Consequently, $\unicode[STIX]{x1D70F}_{n+1}=2\unicode[STIX]{x1D70F}_{n}$. The time taken to reach the $n$th stage is thus

(2.7)$$\begin{eqnarray}\displaystyle t_{n}=\mathop{\sum }_{k=0}^{n-1}\unicode[STIX]{x1D70F}_{k}=\unicode[STIX]{x1D70F}_{0}\mathop{\sum }_{k=0}^{n-1}2^{k}\approx \unicode[STIX]{x1D70F}_{0}2^{n},\quad n\gg 1, & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D70F}_{0}=\unicode[STIX]{x1D6FD}_{\text{rec}}^{-1}R_{0}/B_{\bot ,0}^{\text{tube}}$, $R_{0}$ and $B_{\bot ,0}^{\text{tube}}$ are the radius and perpendicular magnetic field of the initial flux tubes, respectively. Therefore, we can express the index $n$ as a function of time: $n=\log _{2}t/\unicode[STIX]{x1D70F}_{0}$. The explicit formulae for the different quantities as a function of $n$ can be obtained from the recursive relations derived earlier

(2.8)$$\begin{eqnarray}\displaystyle & \displaystyle B_{\bot ,n}^{\text{tube}}=B_{\bot ,0}^{\text{tube}}\,2^{-n/2},\quad E_{M,n}=E_{M,0}2^{-n}, & \displaystyle\end{eqnarray}$$

(2.9)$$\begin{eqnarray}\displaystyle & \displaystyle k_{\bot ,n}=k_{\bot ,0}2^{-n/2},\quad R_{n}=R_{0}2^{n/2},\quad l_{n}=l_{0}2^{n}, & \displaystyle\end{eqnarray}$$

(2.10)$$\begin{eqnarray}\displaystyle & \displaystyle N_{\bot ,n}=N_{\bot ,0}2^{-n}. & \displaystyle\end{eqnarray}$$

By using $n=\log _{2}t/\unicode[STIX]{x1D70F}_{0}$, we obtain the continuous time evolution of the relevant quantities

(2.11)$$\begin{eqnarray}\displaystyle & \displaystyle B_{\bot }^{\text{tube}}=B_{\bot ,0}^{\text{tube}}\,\tilde{t}^{-1/2},\quad E_{M}=E_{M,0}\tilde{t}^{-1}, & \displaystyle\end{eqnarray}$$

(2.12)$$\begin{eqnarray}\displaystyle & \displaystyle k_{\bot }=k_{\bot ,0}\tilde{t}^{-1/2},\quad R=R_{0}\tilde{t}^{1/2},\quad l=l_{0}\tilde{t}, & \displaystyle\end{eqnarray}$$

(2.13)$$\begin{eqnarray}\displaystyle & \displaystyle N_{\bot }=N_{\bot ,0}\tilde{t}^{-1}, & \displaystyle\end{eqnarray}$$

where $\tilde{t}\equiv t/\unicode[STIX]{x1D70F}_{0}$ is time normalized to the initial reconnection time scale and $k_{\bot }\equiv 2\unicode[STIX]{x03C0}/R$ is the perpendicular wavenumber. We note that the aspect ratio of a flux tube $l/R\sim \tilde{t}^{1/2}$. Therefore, as flux tubes merge and grow, they become progressively more slender. However, they always satisfy the critical balance and remain (marginally) stable to kink-type instabilities (appendix B). We also note that $l=l_{0}\tilde{t}=l_{0}(v_{A,0}\unicode[STIX]{x1D6FD}_{\text{rec}}/R_{0})t$. Combined with the critical balance scaling of initial flux tubes $l_{0}/V_{A}\sim R_{0}/v_{A,0}$, where $V_{A}$ is the Alfvén speed based on the guide field, we obtain $l=V_{A}\unicode[STIX]{x1D6FD}_{\text{rec}}t$. We can thus consider $V_{A}\unicode[STIX]{x1D6FD}_{\text{rec}}$ as the characteristic speed for the growth of tube’s length.

To summarize our 3-D hierarchical model: in the perpendicular plane, the merging of flux tubes through reconnection results in the decay of the perpendicular magnetic field, and in the growth of the perpendicular length scale ((2.11) and (2.12)); these conclusions are identical to those obtained from two-dimensional considerations (namely, equation (3) in Z19). In the parallel direction (i.e. along the guide field), the dynamics is (causally) set by critical balance; combined with the perpendicular relations, this yields a prediction for the growth of the flux-tube coherence length in the direction along the guide field.

2.2 Self-similar properties and magnetic power spectra

The power-law time dependencies ((2.11)–(2.13)) suggest that the system evolves in a self-similar manner. We note that the perpendicular dynamics in the 3-D RMHD system is the same as that in its 2-D counterpart (Z19). Therefore, the self-similar properties and the behaviour of the perpendicular magnetic power spectrum should also be identical between the 2-D and 3-D cases. We discuss these properties in more detail in this section, for the completeness of this paper. In the perpendicular plane, the growing length scale and decreasing field strength ((2.11) and (2.12)) can be represented by a dynamical renormalization, which we show as follows. With the power-law time evolution, quantities at any given two moments of time $\tilde{t}$ and $\tilde{t^{\prime }}$ can be related through a scaling factor $\unicode[STIX]{x1D706}=(\tilde{t^{\prime }}/\tilde{t})^{1/2}$

(2.14)$$\begin{eqnarray}\displaystyle k_{\bot }(\tilde{t})=\unicode[STIX]{x1D706}k_{\bot }(\tilde{t^{\prime }}),\quad B_{\bot }^{\text{tube}}(\tilde{t})=\unicode[STIX]{x1D706}B_{\bot }^{\text{tube}}(\tilde{t^{\prime }}),\quad \tilde{t}=\unicode[STIX]{x1D706}^{-2}\tilde{t^{\prime }}. & & \displaystyle\end{eqnarray}$$

That is, the time evolution of the system is equivalent to performing the above scale transformation. In a system consisting of volume-filling flux tubes, the perpendicular magnetic field $\boldsymbol{B}_{\bot }$ follows the same scaling as $B_{\bot }^{\text{tube}}$. That is, $\boldsymbol{B}_{\bot }(\tilde{t})=\unicode[STIX]{x1D706}\boldsymbol{B}_{\bot }(\tilde{t^{\prime }})$.

Let us see how this self-similarity affects the evolution of the perpendicular magnetic power spectrum, which we define as

(2.15)$$\begin{eqnarray}M(k_{\bot },\tilde{t})\equiv \frac{1}{8\unicode[STIX]{x03C0}}\frac{2\unicode[STIX]{x03C0}k_{\bot }}{(2\unicode[STIX]{x03C0})^{2}}\int \text{d}^{2}\boldsymbol{r}\text{e}^{\text{i}\boldsymbol{k}_{\bot }\boldsymbol{\cdot }\boldsymbol{ r}}\langle \boldsymbol{B}_{\bot }(\boldsymbol{x},\tilde{t})\boldsymbol{\cdot }\boldsymbol{B}_{\bot }(\boldsymbol{x}+\boldsymbol{r},\tilde{t})\rangle ,\end{eqnarray}$$

where the angle brackets represent a 3-D average over all $\boldsymbol{x}$, and $\boldsymbol{r}$ is a vector on the perpendicular ($xy$) plane. The spectrum is also averaged over the $z$ dimension, which is not explicitly shown in equation (2.15). Using equation (2.14), the spectra at different times are related as

(2.16)$$\begin{eqnarray}\displaystyle M[k_{\bot }(t^{\prime }),\tilde{t^{\prime }}] & = & \displaystyle M(k_{\bot }/\unicode[STIX]{x1D706},\unicode[STIX]{x1D706}^{2}\tilde{t})=\frac{1}{8\unicode[STIX]{x03C0}}\frac{2\unicode[STIX]{x03C0}(k_{\bot }/\unicode[STIX]{x1D706})}{(2\unicode[STIX]{x03C0})^{2}}\int \text{d}^{2}(\unicode[STIX]{x1D706}\boldsymbol{r})\nonumber\\ \displaystyle & & \displaystyle \text{e}^{\text{i}(\boldsymbol{k}_{\bot }/\unicode[STIX]{x1D706})\cdot (\unicode[STIX]{x1D706}\boldsymbol{r})}\langle \boldsymbol{B}_{\bot }(\unicode[STIX]{x1D706}\boldsymbol{x},\unicode[STIX]{x1D706}^{2}\tilde{t})\boldsymbol{\cdot }\boldsymbol{B}_{\bot }[\unicode[STIX]{x1D706}(\boldsymbol{x}+\boldsymbol{r}),\unicode[STIX]{x1D706}^{2}\tilde{t}]\rangle \nonumber\\ \displaystyle & = & \displaystyle \unicode[STIX]{x1D706}^{-1}\frac{1}{8\unicode[STIX]{x03C0}}\frac{2\unicode[STIX]{x03C0}k_{\bot }}{(2\unicode[STIX]{x03C0})^{2}}\int \text{d}^{2}\boldsymbol{r}\text{e}^{\text{i}\boldsymbol{k}_{\bot }\boldsymbol{\cdot }\boldsymbol{r}}\langle \boldsymbol{B}_{\bot }(\boldsymbol{x},\tilde{t})\boldsymbol{\cdot }\boldsymbol{B}_{\bot }(\boldsymbol{x}+\boldsymbol{r},\tilde{t})\rangle \nonumber\\ \displaystyle & = & \displaystyle \unicode[STIX]{x1D706}^{-1}M(k_{\bot },\tilde{t}).\end{eqnarray}$$

It can be shown that the perpendicular magnetic power spectrum satisfying the relation equation (2.16) has the self-similar solution (Olesen Reference Olesen1997)

(2.17)$$\begin{eqnarray}M(k_{\bot },\tilde{t})=\tilde{t}^{-1/2}\bar{M}(k_{\bot }\tilde{t}^{1/2})=k_{\bot }F(k_{\bot }\tilde{t}^{1/2}),\end{eqnarray}$$

where $\bar{M}$ and $F$ are scaling functions of the variable $k_{\bot }\tilde{t}^{1/2}$.

The solution equation (2.17) represents the pattern of energy decay and scale growth, determined by the physics of the merger events. However, it does not contain any information about energy as a function of scale. Assuming the case of a power-law spectrum, the scaling function has to satisfy $\bar{M}(k_{\bot }\tilde{t}^{1/2})\propto (k_{\bot }\tilde{t}^{1/2})^{-\unicode[STIX]{x1D6FE}}$, where $-\unicode[STIX]{x1D6FE}$ is the index of the spectrum. Equation (2.17) can thus be further expressed as

(2.18)$$\begin{eqnarray}M(k_{\bot },\tilde{t})\propto t^{-1/2}(k_{\bot }\tilde{t}^{1/2})^{-\unicode[STIX]{x1D6FE}}\propto \tilde{t}^{-\unicode[STIX]{x1D6FC}}k_{\bot }^{-\unicode[STIX]{x1D6FE}},\end{eqnarray}$$

where $2\unicode[STIX]{x1D6FC}=\unicode[STIX]{x1D6FE}+1$.

3.1 Simulation set-up

In Viriato, the length scales on the perpendicular plane and in the parallel direction are normalized to (arbitrary) reference lengths $L_{\bot }$ and $L_{\Vert }$ respectively ($L_{\bot }\ll L_{\Vert }$ in the RMHD limit). Times are normalized to the parallel Alfvén time $\unicode[STIX]{x1D70F}_{A}=L_{\Vert }/V_{A}$. In the following description, all quantities are given in dimensionless form. The simulation domain is a box of dimensions $L_{x}\times L_{y}\times L_{z}$, where $L_{x}=L_{y}$ in all the simulations presented. We set $L_{x}=2\unicode[STIX]{x03C0}$, and show the results from simulations with $L_{z}=4L_{x}$. A more elongated simulation box allows more generations of mergers before the length $l$ of the flux tubes in the $z$ direction reaches the length of the box. In order to adequately resolve the dynamics along $z$, and considering the computational cost of the simulation, we find a box with $L_{z}=4L_{x}$ to be optimal. The box has periodic boundary conditions in all directions.

The initial condition is a straightforward extension to three dimensions of our previous 2-D study (Z19). There is no initial flow, $\unicode[STIX]{x1D719}(x,y,z,t=0)=0$, and the magnetic flux tubes are invariant in the $z$ direction: $\unicode[STIX]{x1D713}(x,y,z,t=0)=\unicode[STIX]{x1D713}_{0}\cos (k_{0}x)\cos (k_{0}y)$. It is easy to check that this initial condition satisfies the dissipationless version of equations (3.1)–(3.2) and is therefore an (ideal) equilibrium. We thus have a $2k_{0}\times 2k_{0}$ static array of guide-field-aligned magnetic flux tubes of alternating polarities. In all runs we set $k_{0}=8$ and so the initial radius of each flux tube is $R_{0}=L_{x}/4k_{0}=\unicode[STIX]{x03C0}/16$. We further fix $\unicode[STIX]{x1D713}_{0}k_{0}=1$ to set the strength of the initial in-plane magnetic field $B_{\bot ,0}\equiv \unicode[STIX]{x1D713}_{0}/R_{0}=2/\unicode[STIX]{x03C0}$. A spatially random perturbation with small (linear) amplitude is added to trigger the dynamic evolution of the system. This perturbation breaks the axial symmetry of the system.

In the RMHD limit, the magnetic helicity of the system is $H_{M}=\int \text{d}V\boldsymbol{A}\boldsymbol{\cdot }\boldsymbol{B}\approx B_{g}\int \text{d}V\unicode[STIX]{x1D713}$. We note that, while the local magnetic helicity associated with each flux tube is finite (2.1), the net magnetic helicity of the entire set-up is zero. Since $H_{M}$ is an ideal invariant, it should remain approximately zero throughout the evolution of the system. Therefore, our simulations are relevant to the study of inverse energy transfer in non-helical turbulent systems, a question briefly alluded to in the Introduction and further discussed in appendix C.

We set $\unicode[STIX]{x1D708}=\unicode[STIX]{x1D702}=10^{-4}$ in the simulation, corresponding to the Lundquist number of the initial flux tubes $S_{0}\equiv R_{0}v_{A,0}/\unicode[STIX]{x1D702}=1250$. The perpendicular plane is (spectrally) resolved with $1024^{2}$ grid points. As per our discussion of § 2.1, this value of $S_{0}$ implies that we expect the reconnection of two flux tubes to fall in the SP regime; the thickness of the SP current layers during the first generation of mergers, $\unicode[STIX]{x1D6FF}_{SP,0}\approx S_{0}^{-1/2}R_{0}$ is only marginally resolved with approximately one cell. But we note that as the system evolves and the flux tubes grow to larger scales, the thickness of the current sheets increases (the Lundquist number is expected to remain constant, as we have argued in § 2.1, but the radius of the tubes progressively increases) and the simulation becomes progressively better resolved. Naturally, it would be desirable to investigate initial conditions characterized by even higher values of $S_{0}$, such as to access the plasmoid-mediated reconnection regime. Unfortunately, the computational requirements presented by such simulations are too great: we estimate that achieving $S_{0}\gtrsim 10^{4}$ with the same number of initial flux tubes ($k_{0}=8$) requires a numerical resolution of at least 3000 grid points in each perpendicular dimension. This is too demanding even if the parallel resolution requirements were to remain unchanged.

The parallel length of the box is resolved with 512 grid points. A seventh-order upwind scheme is used to advect the fields in this direction. To check the convergence with respect to parallel resolution, we performed a simulation with the same physical parameters but with $1024^{3}$ grid points instead (thus doubling the resolution in the $z$-direction), and evolved the system for the first few generations of mergers (during which the simulation is least resolved). The results are indistinguishable from the fiducial run with $1024^{2}\times 512$ grid points.

Finally, we also perform a $512^{3}$ run in which the resistivity and viscosity are replaced with hyper-dissipation. Specifically, the dissipative terms in equation (3.1) and equation (3.2) become $\unicode[STIX]{x1D702}_{H}\unicode[STIX]{x1D6FB}_{\bot }^{6}\unicode[STIX]{x1D713}$ and $\unicode[STIX]{x1D708}_{H}\unicode[STIX]{x1D6FB}_{\bot }^{12}\unicode[STIX]{x1D719}$, respectively. Since the time step with which the equations are integrated is constantly recomputed by Viriato to maximize the efficiency of the calculation, it follows that the values of $\unicode[STIX]{x1D702}_{H}$ and $\unicode[STIX]{x1D708}_{H}$ are also continuously adjusted in the code to ensure that the hyper-dissipation is reduced to the minimum possible for a given spatial resolution (Loureiro et al. Reference Loureiro, Dorland, Fazendeiro, Kanekar, Mallet, Vilelas and Zocco2016). However, in the simulation reported here, we observe that the temporal variation of these coefficients is so small that it is justified to consider $\unicode[STIX]{x1D702}_{H}$ and $\unicode[STIX]{x1D708}_{H}$ as constants. The corresponding effective hyper-Lundquist number is $S_{H,0}\equiv R_{0}^{5}v_{A,0}/\unicode[STIX]{x1D702}_{H}\approx 5\times 10^{8}$.

3.2 Visualizations

We begin the analysis of the results of the numerical simulations by showing visualizations at different times of the run with regular dissipation, $S_{0}=1250$. The hyper-dissipative run looks qualitatively similar, except that structures are, of course, sharper. In figure 2 we show the vorticity $\unicode[STIX]{x1D714}_{z}\equiv \unicode[STIX]{x1D6FB}_{\bot }^{2}\unicode[STIX]{x1D719}$ and current density $J_{z}\equiv \unicode[STIX]{x1D6FB}_{\bot }^{2}\unicode[STIX]{x1D713}$ at various times. The top panel shows $\unicode[STIX]{x1D714}_{z}$ on the $xy$ plane at $z=0$ (i.e. $\unicode[STIX]{x1D714}_{z}(x,y,z=0)$), and the middle panel shows $J_{z}$ on that same plane (i.e. $J_{z}(x,y,z=0)$). The bottom panel shows $J_{z}$ on the $xz$ plane at $y=0$ (i.e. $J_{z}(x,y=0,z)$), plotted at the same times as the middle panels.

Figure 2.

Visualizations for the run with $S_{0}=1250$: vorticity density $\unicode[STIX]{x1D714}_{z}$ (normalized by its root mean square value) at various times on the $xy$ planes at $z=0$ (a); current density $J_{z}$ (normalized by its root mean square value) at various times on the $xy$ planes at $z=0$ (b) and on the $xz$ planes at $y=0$ (c).

Figure 3.

Current density contours selected according to $|J_{z}|>3J_{\text{ rms}}$ at various times for the run with $S_{0}=1250$.

These plots suggest a two-stage evolution of the system. The length of the initial flux tubes ($8\unicode[STIX]{x03C0}$) is much longer than the length that critical balance would dictate (approximately $0.3$) and, therefore, the flux tubes are not causally connected from one end of the box to the other. When reconnection events start at different positions along $z$-axis (i.e. on different $xy$ planes), the flux tubes break in the $z$-direction (kink-type instabilities are not expected in our RMHD system, as we discuss in appendix B). Due to this unstable set-up of elongated flux tubes, the system first enters a transient stage of developing turbulence. At this stage, finer structures appear and the system becomes fully turbulent (see the snapshot at $t=11\unicode[STIX]{x1D70F}_{A}$). After the transient stage, the system relaxes to a more natural state, where the initial set-up is ‘washed out’ by the developing turbulence. The system then enters a prolonged stage of decaying turbulence, during which progressively larger flux tubes emerge (both in the $xy$ and in $xz$ planes). The length of the flux tubes grows faster than the radius, consistent with our prediction in (2.12) (this can be seen in figure 2, and a quantitative demonstration is shown in § 3.6).

In figure 3 we show current sheets, defined as the regions where $|J_{z}(x,y,z)|>3J_{\text{rms}}$ and $J_{\text{rms}}$ is the root mean square of $J_{z}$. The increasing complexity at the transient stage is manifested by this plot: the breaking of flux tubes leads to increasing number of current sheets, indicating more structures and reconnection sites. This phenomenon is also shown quantitatively by the probability distribution function of current density square $J_{z}^{2}$ in figure 7, which we will describe in § 3.4.

3.3 Magnetic topology

In this subsection we introduce a method to study the topology of the magnetic field; specifically, we aim to quantify the number of flux tubes on the perpendicular planes so as to make contact with our predictions of § 2.

Any flux tube crossing a given $xy$ plane will leave a local maximum or minimum point of the flux function on the plane (called the O-point). That is, the cross-section of the flux tube on the plane has a structure of a magnetic island. Similarly, the current layers between two merging flux tubes will leave a saddle point on the flux function of the plane (called the X-point). Both the O-points and the X-points are the critical (or stationary) points of a scalar field – in this case, the flux function $\unicode[STIX]{x1D713}(x,y)$ on the plane. In principle, on each $xy$ plane, the number of O-points should be the same as that of X-points. They help characterize in a quantitative way the topology of the magnetic field and can be identified using the technique which we describe next.

We identify the X/O points using the method of Hessian matrix, following previous work (Servidio et al. Reference Servidio, Matthaeus, Shay, Cassak and Dmitruk2009, Reference Servidio, Matthaeus, Shay, Dmitruk, Cassak and Wan2010). On a given $xy$ plane, we first identify all the critical points of the magnetic potential, defined by $\unicode[STIX]{x1D735}_{\bot }\unicode[STIX]{x1D713}=0$, where the perpendicular magnetic field $\boldsymbol{B}_{\bot }=0$. Then we numerically calculate the second-order partial derivatives of the flux function $\unicode[STIX]{x1D713}$ at the critical points, constructing the Hessian matrix

(3.4)$$\begin{eqnarray}H=\left(\begin{array}{@{}cc@{}}\unicode[STIX]{x2202}_{xx}\unicode[STIX]{x1D713} & \unicode[STIX]{x2202}_{xy}\unicode[STIX]{x1D713}\\ \unicode[STIX]{x2202}_{yx}\unicode[STIX]{x1D713} & \unicode[STIX]{x2202}_{yy}\unicode[STIX]{x1D713}\end{array}\right)\end{eqnarray}$$

and computing its determinant and eigenvalues. If $H$ is positive definite at point $(x,y)$, then the point is an O-point with local minimum; and vice versa: a point with negative definite $H$ is an O-point with local maximum. If $H$ has one positive and one negative eigenvalues, such point is a saddle point, called X-point. Rarely, the determinant of $H$ is zero; the locations where this happens are called degenerate critical points.

We apply this diagnostic to time slices throughout the simulation, thus obtaining the time evolution of the number of X/O-points, $N_{\bot }(t)$. For each time slice, we examine eight $xy$ planes spread evenly along the $z$-direction. An example is shown in figure 4. Lastly, the number of X/O-points is averaged over the eight planes (all $xy$ planes should be statistically equivalent).

Figure 4.

Current density $J_{z}$ (colours) and magnetic flux $\unicode[STIX]{x1D713}$ (contours) on the $xy$ plane at $z=L_{z}$, at $t=100\unicode[STIX]{x1D70F}_{A}$ for the run with $S_{0}=1250$. Green circles (‘○’ symbols) identify the positions of O-points; green crosses (‘$\times$’ symbols) mark the positions of X-points.

Figure 5.

Time evolution of the number of X and O points, $N_{\bot }$, for $S_{0}=1250$. The data points show the values of the number of X/O-points averaged over eight $xy$ planes for each timeslice. The shadowed area shows the variation between different $xy$ planes. A $t^{-1}$ power law, predicted in § 2, is shown for reference.

Figure 5 shows the time evolution of $N_{\bot }$ for the run with $S_{0}=1250$. The data points are the number of X/O-points averaged over eight $xy$ planes and the shadowed area shows the variation between different $xy$ planes. Initially, there are $16\times 16=256$ X/O points, as specified by our initial condition. During the transient stage (between ${\sim}6\unicode[STIX]{x1D70F}_{A}-30\unicode[STIX]{x1D70F}_{A}$), the number of X/O points first increases by approximately a factor of 2, and then rapidly decreases. During this stage, the increase in the number of X/O-points is correlated to the deformation of flux tubes and the formation of current sheets (figure 2 and figure 3). This indicates the formation of reconnection sites at multiple places and leads to rapid magnetic energy dissipation, as we will discuss in § 3.4.

After the transient stage, the decay of $N_{\bot }$ slows down and scales with time as $N_{\bot }\sim t^{-1}$, as it predicted by our hierarchical model (2.13). This suggests that the late-time evolution of the system is indeed governed by the merger of magnetic flux tubes.

3.4 Energy decay

From the visualizations and the number of magnetic structures, we can see that the system first goes through a transient stage of developing turbulence. This is followed by a prolonged stage of decaying turbulence, with the formation of progressively larger magnetic structures. The behaviour of magnetic energy evolution, as we will now describe, is consistent with this physical picture.

Figure 6 shows the time evolution of magnetic and kinetic energies for the run with $S_{0}=1250$. Rapid magnetic energy decay occurs during the transient stage (around 6$\unicode[STIX]{x1D70F}_{A}$ to $30\unicode[STIX]{x1D70F}_{A}$), corresponding to the increase, followed by rapid decrease, of the number of X/O points (figure 5). Also, it corresponds to the stage, as we observe from the visualizations, of increasing complexity with the breaking of flux tubes and the emergence of finer structures. Therefore, we attribute the rapid magnetic energy dissipation during this stage to the formation of a large amount of reconnection sites with localized current sheets among the twisted and deformed structures in this complex system. After this stage ($t\gtrsim 30\unicode[STIX]{x1D70F}_{A}$), the magnetic energy decay clearly follows a $t^{-1}$ scaling, as predicted by our hierarchical model, equation (2.11). It coincides with the stage when the number of X/O points exhibits the same scaling, $N_{\bot }\sim t^{-1}$ (figure 5), also in accordance with our predictions, equation (2.13). These numerical confirmations of our predictions suggest that indeed the key process underlying the evolution of our system at this stage is the merger of flux tubes, which appears to proceed as we envisioned in § 2.

Further evidence to support the claims on enhanced energy dissipation during the transient stage is provided in figure 7, which shows the compensated probability distribution function (PDF) of current density square, $J_{z}^{2}$, normalized to its mean square, $J_{\text{ rms}}^{2}$, at different moments of time (in resistive RMHD with constant resistivity, the magnetic energy dissipation (ohmic heating) rate per unit volume is proportional to $J_{z}^{2}$). We arbitrarily choose one snapshot at the transient stage ($t=11\unicode[STIX]{x1D70F}_{A}$) and two from the self-similar stage ($t=65\unicode[STIX]{x1D70F}_{A}$ and $t=214\unicode[STIX]{x1D70F}_{A}$). The PDFs of the two later snapshots (during the self-similar stage) are similar, while the PDF at $t=11\unicode[STIX]{x1D70F}_{A}$ has a clear tail at the high current density end. This suggests that, during the transient stage, there are more localized current sheets facilitating the rapid dissipation of magnetic energy.

The kinetic energy, starting from zero, rapidly increases as the flux tubes start to interact. It then traces the magnetic energy and undergoes the transient rapid-decay stage and $t^{-1}$ self-similar stage. Comparing the vorticity density $\unicode[STIX]{x1D714}_{z}$ and the current density $J_{z}$ on the same $xy$ planes at the same moments of time (figure 2), we find that the flow in the system is dominated by the reconnection outflows with the Alfvén speed based on the perpendicular magnetic field. Therefore, when the magnetic energy decays through the merger of flux tubes, the kinetic energy follows magnetic energy and decays as $t^{-1}$. During the self-similar stage, a small oscillating exchange between magnetic energy and kinetic energy is observed, which we will discuss in § 3.5.

Figure 6.

Time evolution of total magnetic energy ($E_{M}$) and kinetic energy ($E_{K}$) for $S_{0}=1250$. A $t^{-1}$ power law is shown for reference, in accordance with our prediction, equation (2.11).

Figure 7.

Compensated probability distribution function of $J_{z}^{2}/J_{\text{rms}}^{2}$ for $S_{0}=1250$ at three moments of time. The tail at large $J_{z}^{2}/J_{\text{rms}}^{2}$ for the PDF at $t=11\unicode[STIX]{x1D70F}_{A}$ indicates a greater number of localized current sheets in the transient stage.

The rapid magnetic energy decay that we observe in the transient stage bears similarities to that observed in the numerical study of magnetized jets from active galaxies (O’Neill, Beckwith & Begelman Reference O’Neill, Beckwith and Begelman2012; Alves, Zrake & Fiuza Reference Alves, Zrake and Fiuza2018). A jet can be viewed as a screw-pinch/single flux tube that is unstable to current-driven instabilities, such as the kink instability. Such instabilities distort the flux tube and create localized transient current sheets, enhancing the dissipation of magnetic energy. In our system with a strong guide field, we do not expect kink-type instabilities (see discussion in appendix B). However, the configuration of a large ensemble of flux tubes has more degrees of freedom to form reconnection sites via their twisting and deformation.

3.5 Spectra and related quantities

We now proceed to investigate the time evolution of the magnetic and kinetic energy spectra. In order to obtain the clearest spectra possible, we analyse the run with hyper-dissipation instead of normal (Laplacian) diffusivity. For reference, the time evolution of energies for this run is shown in figure 8. Comparison with figure 6 shows qualitatively similar behaviour between the two runs, although two differences stand out: (i) in the hyper-dissipative run the kinetic energy remains approximately a factor of 2 below the magnetic energy, whereas in the $S_{0}=1250$ resistive run they have comparable magnitudes in the self-similar part of the evolution (i.e. for $t/\unicode[STIX]{x1D70F}_{A}\gtrsim 40$); and (ii) the amplitude of the oscillations in the late evolution of the energies is significantly reduced in the hyper-dissipative case. We will first discuss the energy spectra and then turn to the reasons for these differences at the end of this subsection.

Figure 8.

Time evolution of total magnetic ($E_{M}$) and kinetic ($E_{K}$) energies for the run with hyper-dissipation.

Generally speaking, the magnetic energy in the system is contained not only in magnetic flux tubes but also in Alfvén wave packets launched by the convective motion of flux tubes in the perpendicular plane at the scale of flux tubes. Magnetic energy contained in flux tubes is transferred to larger scales through the flux-tube mergers, while that contained in Alfvén wave packets cascades to smaller scales through nonlinear interactions. In the Alfvén wave packets, the kinetic energy is in equipartition with the magnetic energy, cascading to smaller scales. A $k^{-3/2}$ spectrum is expected for both the magnetic and the kinetic energy as a result of the turbulent direct cascade (Boldyrev Reference Boldyrev2006). There is another ingredient in the system that can potentially contribute significantly to the magnetic and kinetic energy spectrum – the thin current sheets between merging flux tubes. The sharp magnetic field reversals at the current sheets are expected to yield a $k^{-2}$ magnetic energy spectrum (Burgers Reference Burgers1948). The Alfvénic outflows from reconnection layers, whose spatial profile is derived in Loureiro, Schekochihin & Uzdensky (Reference Loureiro, Schekochihin and Uzdensky2013), are expected to yield a flat kinetic energy spectrum. That is, the thin current sheets in the system can potentially steepen the turbulent magnetic energy spectrum, and flatten the turbulent kinetic energy spectrum.

With these considerations borne in mind, let us analyse the numerical results. The perpendicular magnetic energy spectrum $M(k_{\bot },t)$ at different moments of time is shown in figure 9(a). The initial spectrum peaks at the scale corresponding to the cross-section of the initial flux tubes. As the flux tubes merge, the peak of the spectrum shifts to larger scales, until it becomes limited by the box size. An inertial range appears to the right of the peak, with a power-law exponent $\unicode[STIX]{x1D6FE}\approx 2$ at earlier times, and $\unicode[STIX]{x1D6FE}\approx 3/2$ at later times. The kinetic spectrum, shown in figure 10, does not have a steady power law. At later times, while the magnetic spectrum peaks at $k_{\bot }=4$, corresponding to the size of islands, the kinetic spectrum peaks at $k_{\bot }=2$, corresponding to the large-scale motion of the islands. Such motion injects kinetic energy into the system, potentially yielding a power-law inertial range at smaller scales, which we tentatively fit with $k_{\bot }^{-3/2}$. However, at earlier times, the kinetic spectrum follows a shallower power law ${\sim}k_{\bot }^{-1}$. While we do not have an explanation for the $-1$ index, we think it is consistent with the observation that at earlier times the system has many localized thin current sheets, where the kinetic energy is dominated by the (Alfvénic) outflows. The outflows have a spatial profile yielding a flat spectrum as we described above, and act as a source of kinetic energy at all scales. As the flux tubes merge and grow to larger scales, the thickness of the current sheets also grows and the magnetic reversal becomes less sharp, resulting in a weaker flattening effect on the spectrum. The kinetic spectrum thus becomes steeper at later times. At early times (before the turbulence is fully developed), the magnetic spectrum follows a $k_{\bot }^{-2}$ spectrum resulting from the sharp field reversals (Burgers Reference Burgers1948). At later times, as the turbulent cascade becomes stronger, the $k_{\bot }^{-2}$ spectrum is replaced by a shallower $k_{\bot }^{-3/2}$ turbulent spectrum. Interestingly, in our 2-D study, a $k_{\bot }^{-2}$ magnetic energy spectrum was maintained throughout the evolution (figure 4 in Z19). The reason for this difference is that, in the 2-D system, the magnetic energy significantly dominates the kinetic energy, making the system less turbulent than its 3-D counterpart. The direct turbulent cascade is too weak to form a $k^{-3/2}$ spectrum, and therefore, the total magnetic energy spectrum in our 2-D system is mostly determined by the current sheets, yielding a $k_{\bot }^{-2}$ spectrum.

Figure 9.

Raw (a) and normalized (b) perpendicular magnetic power spectra, for simulation with hyper-dissipation. Dotted and dashed lines indicate reference slopes of $k_{\bot }^{-2}$ and $k_{\bot }^{-3/2}$, respectively.

Figure 10.

Perpendicular kinetic power spectra for the simulation with hyper-dissipation. $k_{\bot }^{-1}$ and $k_{\bot }^{-3/2}$ slopes are shown for reference.

Our analytical model of § 2 leads to the prediction that the long-time evolution of the magnetic energy spectrum should be self-similar, equation (2.18). To test this prediction, we show in figure 9(b) the spectra normalized to their instantaneous maximum value, $M_{\max }(t)$, as a function of wavenumber normalized to the values $k_{\bot ,\max }(t)$ at which $M_{\max }(t)$ are achieved. We observe that, except for the transient-stage spectrum at $t=10\unicode[STIX]{x1D70F}_{A}$, the spectra thus normalized approximately collapse onto the same distribution, and that this collapse becomes progressively better as time increases. Mathematically, this behaviour is expressed as $M(k_{\bot },t)=M_{\max }(t)\bar{M}(k_{\bot }/k_{\bot ,\max })$, where $\bar{M}(k_{\bot }/k_{\bot ,\max })$ is the static universal distribution to which all the curves collapse after normalization.

We also observe that, after the transient stage, the time evolutions of $k_{\bot ,\max }$ and $M_{\max }$ approximately follow the power-law-in-time behaviour: $k_{\bot ,\max }\propto t^{-1/2}$, and $M_{\max }\propto t^{-1/2}$ – see the second and third panels in figure 11. Thus, the magnetic energy spectrum can be expressed as $M(k_{\bot },t)\propto t^{-1/2}\bar{M}(k_{\bot }t^{1/2})$, in agreement with the prediction of (2.18). Furthermore, since the scaling function is itself a power law in the inertial range, $\bar{M}(k_{\bot }t^{1/2})\propto k_{\bot }^{-3/2}$ ($\unicode[STIX]{x1D6FE}=3/2$), according to our prediction, the time evolution of spectral magnetic energy density at any fixed wavenumber should behave as $M_{k_{\bot }}(t)\propto t^{-\unicode[STIX]{x1D6FC}}$, where $\unicode[STIX]{x1D6FC}=(\unicode[STIX]{x1D6FE}+1)/2=5/4$. This is indeed observed in the bottom panel of figure 11, where after the transient stage, $M_{k_{\bot }}(t)$ (for moderate wavenumbers) approaches the $t^{-5/4}$ power law. At the largest scale the energy keeps building up, explicitly showing the inverse transfer of magnetic energy. From the measured scaling laws described above, we conclude that the self-similar solution of the magnetic power spectrum, based on our hierarchical analytical model, is confirmed by our numerical results.

Figure 11.

Time evolution of $\unicode[STIX]{x1D709}_{M}$ (a), $k_{\bot ,\max }$ (b), $M_{\max }$ (c) and $M_{k_{\bot }}$ (for selected values of $k$; d), for the simulation with hyper-dissipation.

Complementary to the time evolution of $k_{\bot ,\max }$, the growth of the perpendicular coherence length of our system is illustrated by computing the magnetic-energy-containing scale, defined as

(3.5)$$\begin{eqnarray}\unicode[STIX]{x1D709}_{M}(t)\equiv \frac{\displaystyle \int 2\unicode[STIX]{x03C0}k_{\bot }^{-1}M(k_{\bot },t)\,\text{d}k_{\bot }}{\displaystyle \int M(k_{\bot },t)\,\text{d}k_{\bot }}.\end{eqnarray}$$

As shown in the top panel of figure 11, after the transient stage, $\unicode[STIX]{x1D709}_{M}$ increases with time roughly following the $t^{1/2}$ scaling, consistent with our model, until the later stage ($t\gtrsim 100\unicode[STIX]{x1D70F}_{A}$) during which the parallel length of flux tubes becomes limited by the length of the box.

To characterize the system’s dynamics along the guide field, we show in figure 12 the magnetic power spectrum as a function of $k_{z}$. A broad energy distribution develops once the flux tubes break in the $z$ direction. As flux tubes merge, magnetic energy at high $k_{z}$ rapidly decreases, indicating the increasing coherence length in $z$ direction. This is consistent with the measurement of magnetic structure functions (figure 13), showing the rapid growth of the characteristic flux tube, as we explain in § 3.6.

Finally, let us address the issue of the oscillatory energy exchange between the magnetic field and the flow. From the time evolution of magnetic energy density at different wavenumbers $k_{\bot }$ (figure 11d), we see that energy oscillation mostly happens at the smallest $k_{\bot }$, corresponding to system-size magnetic flux tubes. These tubes cannot merge and grow into larger scales due to the system-size limit (in a sense, they have reached a final ‘equilibrium’ consisting of two flux tubes of opposite polarities). Instead, they exchange energy with the flow – possibly both in the form of Alfvén waves, and bouncing motions: background flows can push these flux tubes against each other, leading to their (incompressible) deformation. This increases magnetic field line tension and results in a restitution force leading to an oscillation whose amplitude, therefore, depends on the ratio of the kinetic-to-magnetic energy. Since the magnetic energy is larger in the hyper-dissipative case than in the Laplacian resistivity case (simply because hyper-resistivity leads to better conservation of magnetic energy), the hyper-dissipation case has smaller amplitude oscillations.

Figure 12.

Parallel (to the guide-field $\boldsymbol{B}_{g}$) magnetic power spectra at different moments of time for the simulation with hyper-dissipation.

3.6 Magnetic structure functions

Figure 13.

Contours of normalized (to its maximum value) magnetic structure function $S_{B}^{(2)}(\unicode[STIX]{x1D6FF}x_{\bot },\unicode[STIX]{x1D6FF}x_{\Vert })$ at various times for the run with hyper-dissipation.

In order to obtain statistical information on the growth of flux tubes, and to test the assumption of critical balance, we calculate the second-order structure function of the magnetic field $\boldsymbol{B}(\boldsymbol{x},t)$, defined as

(3.6)$$\begin{eqnarray}S_{B}^{(2)}(\unicode[STIX]{x1D6FF}\boldsymbol{x},t)\equiv \langle |\boldsymbol{B}(\boldsymbol{x}_{2},t)-\boldsymbol{B}(\boldsymbol{x}_{1},t)|^{2}\rangle _{x_{1}},\end{eqnarray}$$

where $\unicode[STIX]{x1D6FF}\boldsymbol{x}\equiv \boldsymbol{x}_{2}-\boldsymbol{x}_{1}$ and the angle brackets indicate the average over the position $\boldsymbol{x}_{1}$ in the simulation domain. We define the local mean field as $\boldsymbol{B}_{m}\equiv [\boldsymbol{B}(\boldsymbol{x}_{1})+\boldsymbol{B}(\boldsymbol{x}_{2})]/2$ and the local unit magnetic field vector $\hat{\boldsymbol{b}}_{m}\equiv \boldsymbol{B}_{m}/|\boldsymbol{B}_{m}|$ represents its direction. These definitions ensure that the structure function is computed with respect to the local mean magnetic field (as it should, if it is to be used to gauge whether the turbulence we observe is critically balanced (Cho & Vishniac Reference Cho and Vishniac2000; Mallet et al. Reference Mallet, Schekochihin, Chandran, Chen, Horbury, Wicks and Greenan2016)), rather than the global guide field. The local coordinates ($\unicode[STIX]{x1D6FF}x_{\bot },\unicode[STIX]{x1D6FF}x_{\Vert }$) can thus be specified by

(3.7)$$\begin{eqnarray}\unicode[STIX]{x1D6FF}x_{\bot }\equiv |\hat{\boldsymbol{b}}_{m}\times \unicode[STIX]{x1D6FF}\boldsymbol{x}|,\quad \unicode[STIX]{x1D6FF}x_{\Vert }\equiv \hat{\boldsymbol{b}}_{m}\boldsymbol{\cdot }\unicode[STIX]{x1D6FF}\boldsymbol{x}.\end{eqnarray}$$

The correlation lengths of the magnetic field parallel and perpendicular to the local mean field, which we interpret as the statistical length and radius of the flux tubes respectively, are thus expressed by the contours of the structure function $S_{B}^{(2)}(\unicode[STIX]{x1D6FF}x_{\bot },\unicode[STIX]{x1D6FF}x_{\Vert })$ on the ($\unicode[STIX]{x1D6FF}x_{\bot },\unicode[STIX]{x1D6FF}x_{\Vert }$) plane. This is shown in figure 13 (normalized to the instantaneous maximum value of the structure function). We then define the characteristic length $l$ and the radius $R$ of the flux tubes using the intersection of the contour line $S_{B}^{(2)}/S_{\text{B},\max }^{(2)}=0.7$ (roughly the largest value at which the contour lines are not affected by noise) with the $\unicode[STIX]{x1D6FF}x_{\Vert }$ axis and the $\unicode[STIX]{x1D6FF}x_{\bot }$ axis, respectively. Note that because of periodic boundary conditions, the largest values of $\unicode[STIX]{x1D6FF}x_{\bot }$ and $\unicode[STIX]{x1D6FF}x_{\Vert }$, and hence maximum radius and length of a flux tube that can be attained in our simulations, are roughly $L_{x}/\sqrt{2}=\sqrt{2}\unicode[STIX]{x03C0}$ and $L_{z}/2=4\unicode[STIX]{x03C0}$, respectively. The growth of the typical flux tube can be clearly observed as the system evolves with time. The growth of parallel length is faster than that in the perpendicular direction, which agrees with our discussion in § 2 that the aspect ratio of a flux tube increases with time ($l/R\sim \tilde{t}^{1/2}$).

Figure 14.

Time evolution of $l$ (a, dashed line represents a $t^{1}$ power law, equation (2.5)) and scaling of $l$ versus $R/\bar{B}_{\bot }$ (b, dashed line represents the critical balance scaling, equation (2.4)) for the run with hyper-dissipation, where $l$ and $R$ are the length and radius of the characteristic flux tube respectively.

In order to characterize the time evolution of magnetic flux tubes, we take eight snapshots between $t=10\unicode[STIX]{x1D70F}_{A}$ and $t=80\unicode[STIX]{x1D70F}_{A}$ (roughly every $10\unicode[STIX]{x1D70F}_{A}$). For each snapshot, we compute the values of $l$ and $R$. The time evolution of $R$ (not shown) roughly follows a $t^{1/2}$ scaling (similar to the time evolution of $\unicode[STIX]{x1D709}_{M}$ shown in figure 11), whereas the time evolution of $l$ scales linearly with time (figure 14a), in accordance with the prediction of equation (2.12). We also compute the root-mean-square value of the perpendicular magnetic field, $\bar{B}_{\bot }$, and test if it satisfies the critical balance scaling, $l/\hat{B}_{g}\sim R/\bar{B}_{\bot }$, where $\hat{B}_{g}=B_{g}L_{\bot }/L_{\Vert }=1$ is the guide field normalized in the same way as used to compute the structure function. We plot $l$ versus $R/\bar{B}_{\bot }$ for these eight snapshots in figure 14(b). As seen, the critical balance condition is indeed verified (except for the last snapshot, whereupon the parallel length of flux tubes has reached the length of the simulation box, and the critical balance scaling is thus expected to break). This confirms our assumption made in § 2 that in our system of merging flux tubes, the aspect ratio of the flux tubes, $l/R$, is determined by critical balance.

3.7 Energy transfer function

In this subsection we use shell-to-shell energy transfer functions to quantify the transfer between magnetic and kinetic energy at different scales. This enables us to identify direct and inverse energy transfer in our simulations. This technique has been used in studies of MHD and kinetic turbulence (e.g. Alexakis, Mininni & Pouquet (Reference Alexakis, Mininni and Pouquet2005), Told et al. (Reference Told, Jenko, TenBarge, Howes and Hammett2015)), and we made modifications for its use in the RMHD system.

We define the shell-filtered variables as

(3.8)$$\begin{eqnarray}\boldsymbol{B}_{\bot }^{K}(\boldsymbol{x}_{\bot },z)=\mathop{\sum }_{K<k_{\bot }\leqslant K+1}\boldsymbol{B}_{\bot }(\boldsymbol{k}_{\bot },z)\text{e}^{\text{i}\boldsymbol{k}_{\bot }\boldsymbol{\cdot }\boldsymbol{ x}_{\bot }}\end{eqnarray}$$

and

(3.9)$$\begin{eqnarray}\boldsymbol{u}_{\bot }^{K}(\boldsymbol{x}_{\bot },z)=\mathop{\sum }_{K<k_{\bot }\leqslant K+1}\boldsymbol{u}_{\bot }(\boldsymbol{k}_{\bot },z)\text{e}^{\text{i}\boldsymbol{k}_{\bot }\boldsymbol{\cdot }\boldsymbol{ x}_{\bot }}.\end{eqnarray}$$

Here, $\boldsymbol{x}_{\bot }$ is the projection of the position vector onto the $xy$ plane and the Fourier and inverse Fourier transforms are only performed in the $x$- and $y$- directions. The shell-filtered variables contain only the magnetic and velocity fields within a certain range of perpendicular scales. The shell transfer functions are obtained by calculating the evolution equation of energy in each shell averaged over the $z$ direction. We are interested in the time evolution of magnetic energy in each shell $K$

(3.10)$$\begin{eqnarray}\displaystyle & & \displaystyle \unicode[STIX]{x2202}_{t}\left\langle \int \text{d}^{2}\boldsymbol{x}_{\bot }\frac{1}{2}(\boldsymbol{B}_{\bot }^{K})^{2}\right\rangle _{z}\nonumber\\ \displaystyle & & \displaystyle \quad =\left\langle \int \text{d}^{2}\boldsymbol{x}_{\bot }\mathop{\sum }_{Q}[\boldsymbol{B}_{\bot }^{K}\boldsymbol{\cdot }(\boldsymbol{B}_{\bot }\boldsymbol{\cdot }\unicode[STIX]{x1D735}_{\bot })\boldsymbol{u}_{\bot }^{Q}-\boldsymbol{B}_{\bot }^{K}\boldsymbol{\cdot }(\boldsymbol{u}_{\bot }\boldsymbol{\cdot }\unicode[STIX]{x1D735}_{\bot })\boldsymbol{B}_{\bot }^{Q}]-\unicode[STIX]{x1D702}|\unicode[STIX]{x1D6FB}_{\bot }\boldsymbol{B}^{K}|^{2}\right\rangle _{z}.\qquad\end{eqnarray}$$

The transfer of magnetic energy from shell $Q$ to shell $K$ is

(3.11)$$\begin{eqnarray}T_{bb}(Q,K)=-\left\langle \int \text{d}^{2}\boldsymbol{x}_{\bot }\boldsymbol{B}_{\bot }^{K}\boldsymbol{\cdot }(\boldsymbol{u}_{\bot }\boldsymbol{\cdot }\unicode[STIX]{x1D735}_{\bot })\boldsymbol{B}_{\bot }^{Q}\right\rangle _{z},\end{eqnarray}$$

and the transfer of kinetic energy in shell $Q$ to magnetic energy in shell $K$ is

(3.12)$$\begin{eqnarray}T_{ub}(Q,K)=\left\langle \int \text{d}^{2}\boldsymbol{x}_{\bot }\boldsymbol{B}_{\bot }^{K}\boldsymbol{\cdot }(\boldsymbol{B}_{\bot }\boldsymbol{\cdot }\unicode[STIX]{x1D735}_{\bot })\boldsymbol{u}_{\bot }^{Q}\right\rangle _{z}.\end{eqnarray}$$

Figure 15.

Averaged energy transfer from shell $Q$ to shell $K$ at $t=10\unicode[STIX]{x1D70F}_{A}$ (first stage) and $t=60\unicode[STIX]{x1D70F}_{A}$ (second stage) for the run with hyper-dissipation. (a) magnetic-to-magnetic transfer $T_{bb}(Q,K)$. (b) kinetic-to-magnetic transfer $T_{ub}(Q,K)$.

We perform the energy transfer function diagnostic on the results of the run with hyper-dissipation. We arbitrarily choose one snapshot during the first stage of developing turbulence ($t=10\unicode[STIX]{x1D70F}_{A}$) and another during the second stage of decaying turbulence ($t=60\unicode[STIX]{x1D70F}_{A}$) to illustrate the behaviour of shell-to-shell energy transfer, as shown in figure 15. To guide the reader, we note that red colour below (above) the diagonal in these plots denotes an inverse (direct) energy transfer; and vice versa for the blue colour. We first focus on the contour plot of $T_{bb}(Q,K)$ at $t=10\unicode[STIX]{x1D70F}_{A}$, shown in the top-left panel. At this moment, the magnetic energy spectrum shown in figure 9 peaks at $k_{\bot }=8$. There is a strong direct energy transfer from shells $Q=7,8$ to shells $K>Q$, indicating the non-local direct transfer characterizing the formation of thin current sheets (consistent with the $k_{\bot }^{-2}$ magnetic energy spectrum as we explained earlier, shown in figure 9a). There is also a weak but statistically significant inverse transfer of magnetic energy from shells $Q=7,8$ to shells $K=5-7$ (corresponding to the peak of magnetic spectrum). This suggests that the inverse transfer happens at the scale of flux tubes and is caused by flux-tube merger. Due to the development of turbulence, $T_{bb}$ is dominated by the direct transfer from small-$Q$ shells to large-$K$ shells overall. At this stage, there is also a transfer of magnetic energy to kinetic energy close to the diagonal, as shown in the contour plot of $T_{ub}$ in the bottom-left panel, consistent with the formation of plasma flow.

At $t=60\unicode[STIX]{x1D70F}_{A}$, flux tubes have reached larger scales (with a peak at around $k_{\bot }=4$). Shells $K,Q>3$ are dominated by the strong direct transfer between a wide range of scales. However, although relatively weak in strength, there is also a finite inverse magnetic energy transfer from shells $Q>3$ to shells $K=1-2$ shown by the red horizontal band along the lower boundary of the contour plot of $T_{bb}$ in the top-right panel. From the $T_{ub}$ plot in the bottom-right panel, we observe a mild transfer from magnetic energy in shells $K=1,2$ to kinetic energy in shells $Q>4$. For $K>3$, kinetic energy in shells $Q$ is transferred to magnetic energy in shells $K$ if $Q<K$, while for $Q>K$, the transfer between kinetic and magnetic energy happens in both directions.

In summary, the energy transfer diagnostics paint a complex picture that is consistent with the analytical model and numerical results discussed earlier. Both inverse and direct energy transfer occurs: the former driven by flux-tube coalescence; the latter by the turbulent dynamics that arises. While the direct energy transfer is dominant, it is the inverse transfer that progressively delivers energy to larger and larger scales until reaching the simulation domain size.