2.1. Expression for $\mathcal {C}_\textrm {{Snk}}(\psi,A,t)$

We assume that any single coalescence event takes place between two islands only, i.e. we effectively treat island–island interaction as binary. This assumption imposes an exclusion principle on islands undergoing merging. Once an island starts to merge with another island, it is unavailable for other merging events until the current coalescence process finishes.Footnote ³ When one island with area $A$ and flux $\psi$ merges with another island with area $A'$ and flux $\psi '$, the resulting island will possess an area $A+A'$ and a flux with the larger value between $\psi$ and $\psi '$ (Fermo et al. Reference Fermo, Drake and Swisdak2010).

The rate of change of $f$ can be calculated using the 2-D hard-sphere scattering model. According to our normalisation, the number density of islands (in $A-\psi$ space) with flux $\psi$ and area $A$ is $f(\psi,A,t)$. The number of islands with arbitrary $A'$ and $\psi '$ is given by $\int _0^\infty \,\textrm {d}A' \int _0^\infty \,\textrm {d}\psi ' f(\psi ',A',t)$. The probability rate for two islands to meet and merge in real 2-D space is given by $\sigma \delta v /L^2$, where the cross-section for an island with area $A$ interacting with an island with area $A'$ is $\sigma (A,A^\prime ) = \sqrt {A}+\sqrt {A'}$, multiplied by a geometric factor of order unity, which we take to be one for simplicity.Footnote ⁴ The relative velocity between two islands, $\delta v$, is approximated with the larger Alfvén velocity of each of them, $\delta v = \max \{\psi /\sqrt {A},\psi '/\sqrt {A'}\}$.Footnote ⁵ With such an approximation, the maximum possible error in the relative velocity is $\delta v$.

Using the model described above, the sink term $\mathcal {C}_{\text {Snk}}(\psi,A,t)$, which equals the decreasing rate of $f(\psi,A)$, can be expressed as

(2.3)\begin{equation} \mathcal{C}_{\text{Snk}}(\psi,A,t) ={-}\frac{1}{L^2} f(\psi,A,t)\int_0^\infty\,\textrm{d}A' \int_0^\infty\,\textrm{d}\psi'\ f(\psi',A',t) \sigma(A,A^\prime) \delta v. \end{equation}

We perform a simple dimensional analysis to show that our IKE (dimensionally $\partial _t f \sim \mathcal {C}_\textrm {snk}$) has the standard form $\partial _t f \sim \nu f$, where $\nu$ is a characteristic changing rate of $f$. The mean-free path of the system is $\lambda _\textrm {mfp}\sim 1/(n \sigma )$, where $n=N/L^2=(\int f\,\textrm {d}A\,\textrm {d}\psi )/L^2$ is the island density in configuration space. Plugging the scaling $\int f\,\textrm {d}A\,\textrm {d}\psi \sim L^2/(\lambda _\textrm {mfp} \sigma )$ into (2.3), we obtain $\mathcal {C}_\textrm {snk}(A,\psi,t) \sim (\delta v/\lambda _\textrm {mfp}) f(\psi,A,t)$ and hence the characteristic decreasing rate of $f$ is as expected: $\nu \sim \delta v/\lambda _\textrm {mfp}$. The mean-free path is determined both by the density and typical cross-section of islands. In the case of volume-filling islands, the density of islands is $n=N/L^2 \sim 1/ \langle{A}\rangle $, where $ \langle{A}\rangle $ is the typical area of the islands. The typical cross-section of the islands is thus $\sqrt {\langle {A}\rangle}$, which gives rise to $\lambda _\textrm {mfp} \sim 1/(n \sigma ) \sim \sqrt {\langle {A}\rangle}$. This is consistent with the intuition that when the system is packed with islands, each island will interact with an adjacent one, and $\lambda _\textrm {mfp}$ should be comparable to the size of islands. In this case, the convection time of islands ($\sim \lambda _\textrm {mfp}/v_A$) is the same as their local Alfvén time ($\sim \sqrt {\langle {A}\rangle}$), where $v_A$ is the local Alfvén velocity determined by the magnetic field of islands. In our IKE, $\lambda _\textrm {mfp}$ can be adjusted by tuning the size and density of islands, i.e. changing the ‘volume-filling fraction’ of islands in the system. We will define this quantity in § 2.5.

2.2. Time delay between $\mathcal {C}_\textrm {{Snk}}(\psi,A,t)$ and $\mathcal {C}_\textrm {{Src}}(\psi,A,t)$

The source term $\mathcal {C}_{\text {Src}}(\psi,A,t)$ can be obtained in a similar way, except that there is an additional feature that we now explain. Any given merging event is described by both a sink term and a source term in the collision operator. In traditional collision theory, collisions (between particles) are regarded as instantaneous, because the particle interaction time is short compared with the particle flight time between collisions. In our island coalescence problem, in contrast, the coalescing process is slow, taking a long time compared with the advection of islands. The coalescence time of islands thus becomes an important dynamical time scale in our system. Therefore, it is necessary to consider the time difference between the sink term and source term for any given merger, where the former corresponds to the start of a merger and the latter corresponds to the end. The time delay for a given merging event is denoted by $\tau _\textrm {rec}$, which depends on the areas and fluxes of both islands. The expression for $\tau _\textrm {rec}$ will be given in § 2.3 and the way we treat the islands during the merger will be described in § 2.4.

We assume that the decrease in the number of islands happens as soon as the two islands meet in real space. However, the increase in the number of islands happens when the ‘parent’ islands finish the merging process. That is, for a merging event between two islands characterised by ($\psi,A$) and ($\psi ',A'$) (with $\psi '<\psi$), starting at time $t$, the sink term acts to decrease $f(\psi,A,t)$ and $f(\psi ',A',t)$ immediately, while the source term acts to increase $f(\psi, A+A', t+\tau _\textrm {rec})$ when islands finish merging at $t+\tau _\textrm {rec}(\psi,A,\psi ',A')$. Equivalently, for a source term increasing $f(\psi,A,t)$ at time $t$, the corresponding sink term resulting from the same merging event decreases $f(\psi,A',t-\tau _\textrm {rec})$ and $f(\psi ',A-A',t-\tau _\textrm {rec})$ when the islands with $(\psi,A')$ and $(\psi ',A-A')$ start to merge at $t-\tau _\textrm {rec}(\psi,A',\psi ',A-A')$.

The way we implement this time delay $\tau _\textrm {rec}$ in the source term $\mathcal {C}_{\text {Src}}(\psi,A,t)$ is as follows. At any given time $t$ and given point $(\psi,A)$ in the phase space, we trace back the history of island distributions at all previous times $t-\tau$ ($\tau \in [0,t]$) and consider all possible pairs of islands ($\psi,A'$) and ($\psi ',A-A'$) from which the new islands can be formed. The merging times $\tau _\textrm {rec}(\psi, A', \psi ', A-A')$ of these islands are calculated. Those islands whose merging time $\tau _\textrm {rec}$ matches $\tau$ contribute to the source term.

Following the above description, the source term can thus be formally expressed as

(2.4)\begin{gather} \mathcal{C}_{\text{Src}}(\psi,A,t) = \frac{1}{L^2} \int_0^t \,\textrm{d}\tau \int_0^A\,\textrm{d}A'\ f(\psi,A',t-\tau)\int_0^\psi\,\textrm{d}\psi' \nonumber\\ f(\psi',A-A',t-\tau)\, \sigma(A-A^\prime, A^\prime)\, \delta v\, \delta[\tau-\tau_\textrm{rec}(\psi, A', \psi', A-A')]. \end{gather}

We note that introducing the delay of islands’ merger time is important to obtain the correct dynamical time scale of the evolution of the system. With an instantaneous collision operator, the system would evolve on the Alfvénic time scale, whereas with the consideration of merger time delay, the system evolves on the reconnection time scale (which has indeed been shown to be the correct evolution time via direct numerical simulations of this problem Zhou et al. Reference Zhou, Bhat, Loureiro and Uzdensky2019, Reference Zhou, Loureiro and Uzdensky2020; Bhat et al. Reference Bhat, Zhou and Loureiro2021; Hosking & Schekochihin Reference Hosking and Schekochihin2021).

2.3. Reconnection time $\tau _\textrm {rec}$

The time for two islands to merge, $\tau _\textrm {rec}$, is determined by the physics of the reconnection process. In this subsection, we first consider the collisional, low Lundquist number regime $S\lesssim 10 ^4$ and calculate the reconnection time $\tau _\textrm {rec}$ based on the Sweet–Parker (SP) reconnection model (Parker Reference Parker1957; Sweet Reference Sweet1958). We note that different reconnection models can be adopted to calculate $\tau _\textrm {rec}$ and our IKE model can be modified to study systems in different reconnection regimes.

In the SP regime, the merging of two islands with different sizes ($A_1,A_2$) and fluxes ($\psi _1,\psi _2$) is a process of asymmetric reconnection (the symmetric case is a particular scenario in this description), where the outflow Alfvénic velocity $v_A$ and the pertinent Lundquist number $S$ are calculated using geometric averages (Cassak & Shay Reference Cassak and Shay2007), as follows:

(2.5)\begin{gather} v_A = \left(B_1 B_2\right)^{1/2} = \left(\frac{\psi_1}{\sqrt{A_1/{\rm \pi}}}\frac{\psi_2}{\sqrt{A_2/{\rm \pi}}}\right)^{1/2} \simeq \frac{\left(\psi_1 \psi_2\right)^{1/2}}{\left(A_1A_2\right)^{1/4}}, \end{gather}

(2.6)\begin{gather} S = \frac{\sqrt{\psi_1\psi_2}}{\eta}. \end{gather}

The merging time between an island with ($\psi _1,A_1$) and island with ($\psi _2,A_2$) is calculated as

(2.7)\begin{equation} \tau_\textrm{rec}(\psi_1, A_1, \psi_2, A_2) \simeq \frac{\min \left( \sqrt{A_1},\sqrt{A_2}\right) } {v_\textrm{rec}}, \end{equation}

where $v_\textrm {rec}$ is the merging velocity of the two islands. The $v_\textrm {rec}$ and outflow Alfvénic velocity $v_{A}$ are related by the dimensionless reconnection rate $\beta _\textrm {rec}=v_\textrm {rec}/v_A$. In the SP model, $\beta _\textrm {rec}$ is related to the local Lundquist number of the merging islands $S$ by $\beta _\textrm {rec} = S^{-1/2}$. Therefore,

(2.8)\begin{equation} v_\textrm{rec} = \beta_\textrm{rec} v_A = S^{{-}1/2} v_A \simeq \left(\frac{\sqrt{\psi_1\psi_2}}{\eta}\right)^{{-}1/2} \frac{\left(\psi_1 \psi_2\right)^{1/2}}{\left(A_1A_2\right)^{1/4}} = \eta^{1/2}\left(\frac{\psi_1\psi_2}{A_1A_2}\right)^{1/4}, \end{equation}

and the reconnection time $\tau _\textrm {rec}$ can be calculated.

We emphasise that the elements of the reconnection process that affect the IKE are the conservation of magnetic flux and area (both of which are quite generally valid), and the reconnection rate. Beyond this, the actual detailed reconnection physics – such as, for example, the specific mechanisms that set the reconnection rate – do not enter our model. Therefore, our IKE can be straightforwardly extended to other reconnection regimes where the reconnection rate is a constant. This includes collisionless plasmas or, indeed, any regime where Hall physics is important; in either case, the reconnection rate is $\beta _\textrm {rec} \simeq 0.1$ (see, e.g. Biskamp, Schwarz & Drake Reference Biskamp, Schwarz and Drake1995; Birn et al. Reference Birn, Drake, Shay, Rogers, Denton, Hesse, Kuznetsova, Ma, Bhattacharjee and Otto2001; Shay et al. Reference Shay, Drake, Rogers and Denton2001; Cassak, Liu & Shay Reference Cassak, Liu and Shay2017). For convenience, in the rest of the paper, we refer to this regime with $\beta _\textrm {rec} \simeq 0.1$ as the collisionless regime but, we stress, it is in fact more general than that.

2.4. Accounting for islands undergoing merging processes

Owing to the time delay between the sink term and the source term of the collision operator, the contributions of currently merging islands to the distribution function $f(\psi,A,t)$ have already been subtracted by the sink term, while the components for the resulting new islands have not yet been added by the source term. This would lead to the problem that, when calculating global physical quantities using weighted averages with the distribution function, the islands in the merging process would not be taken into account.

To fix this problem, we add the components of the distribution function that have been sunk but not yet re-added by the delayed source term. We denote by $\tilde {f}(\psi,A,t)$ the distribution function of the islands that are undergoing the merging process; as stated above, under the assumption of binary merging, these islands do not participate in the interaction with additional islands. The total distribution function is then given by $f_\textrm {tot}(\psi,A,t) = f(\psi,A,t) + \tilde {f}(\psi,A,t)$; all the physical quantities should be calculated using this total distribution function.

A kinetic equation for the time evolution of $\tilde {f}(\psi,A,t)$ may also be calculated from the merging of islands and can be decomposed into the source term $\mathcal {\tilde {C}}_\textrm {Src}(\psi,A,t)$ and the sink term $\mathcal {\tilde {C}}_\textrm {Snk}(\psi,A,t)$, as follows:

(2.9)\begin{equation} \partial_t \tilde{f}(\psi,A,t) = \mathcal{\tilde{C}}_\textrm{Src}(\psi,A,t) + \mathcal{\tilde{C}}_\textrm{Snk}(\psi,A,t). \end{equation}

We note that the advective terms $\partial _\psi (\dot {\psi }\tilde {f})$ and $\partial _A(\dot {A}\tilde {f})$ are neglected in (2.9). That is, we do not consider the continuous evolution of the merging island parameters during the reconnection process. Instead, we assume that when islands start to merge, their parameters remain the same until the moment when islands finish merging. At the end of mergers, the islands’ parameters jump from values of previous islands to values of the new islands (in particular, one of the islands disappears). This is implemented by the collisional integrals $\mathcal {\tilde {C}}_\textrm {Src}$ and $\mathcal {\tilde {C}}_\textrm {Snk}$. On a time scale much longer than the life time of one island generation, the replacement of the continuous change of islands parameters during mergers with a discontinuous change at the end of mergers will not significantly change the system evolution.

The merging islands distribution $\tilde {f}(\psi,A,t)$ increases when islands start to merge; therefore, the source term $\mathcal {\tilde {C}}_\textrm {Src}(\psi,A,t)$ should act immediately on ${\tilde {f}}(\psi,A,t)$. Obviously $\mathcal {\tilde {C}}_\textrm {Src}(\psi,A,t)=-\mathcal {C}_{\text {Snk}}(\psi,A,t)$, because the components of islands starting to merge should be subtracted from $f(\psi,A,t)$ and added to ${\tilde {f}}(\psi,A,t)$ instantaneously.

The sink term is more complex. We want to subtract from $\tilde {f}(\psi,A,t)$ once the islands finish merging. That is, at any given time $t$ when islands finish merging, we want to subtract the components that we added at $t-\tau _\textrm {rec}$ when they started to merge. This can be implemented by introducing the time delay to the source term $\mathcal {\tilde {C}}_\textrm {Src}(\psi,A,t)$:

(2.10)\begin{align} \mathcal{\tilde{C}}_\textrm{Snk}(\psi,A,t) & ={-}\int_0^t \,\textrm{d}\tau\, \mathcal{\tilde{C}}_\textrm{Src}(\psi,A,t-\tau)\,\delta[\tau-\tau_\textrm{rec}(\psi, A', \psi', A)]\nonumber\\ & ={-}\frac{1}{L^2} \int_0^t \,\textrm{d}\tau f(\psi,A,t-\tau) \int_0^\infty\,\textrm{d}A' \int_0^\infty\,\textrm{d}\psi' f(\psi',A',t-\tau)\, \sigma(A,A^\prime)\, \delta v \notag\\ & \quad \times \delta[\tau-\tau_\textrm{rec}(\psi, A', \psi', A)]. \end{align}

One of our merging rules says that the area of the resulting islands should be the sum of the areas of the two merging islands. Therefore, in this model, the total area of all the islands should be conserved. The total area can be calculated by

(2.11)\begin{equation} A_\textrm{tot}=\int\,\textrm{d}A\,\textrm{d}\psi\ f_\textrm{tot}(\psi,A,t)A. \end{equation}

In the numerical implementation of our IKE that we discuss in §§ 3 and 4, the conservation of $A_\textrm {tot}$ has been tested and confirmed numerically as a benchmark of the equation solver.

2.5. Macroscopic quantities and magnetic spectrum

We first introduce the parameter $V_\textrm {fill} \equiv A_\textrm {tot}/L^2 \in [0,1]$ to represent the volume-filling fraction of islands in the system. This quantity is determined by both the density and sizes of islands. As we discussed at the end of § 2.1, larger $V_\textrm {fill}$ implies a smaller mean-free path $\lambda _\textrm {mfp}$ for the islands. When the islands are volume filling, i.e. $V_\textrm {fill} \approx 1$, the mean-free path is close to the typical length scale of islands $\lambda _\textrm {mfp}\approx \sqrt {A_\textrm {tot}/N}$. This is the case that we focus on in this paper, but we also show in § 4.1 a scan in $V_\textrm {fill}$, which indicates that this parameter does not control the evolution of the main quantities of interest in our study.

The other main macroscopic quantities of interest are the total number of islands $N$, the average area of islands $ \langle {A}\rangle $ and the total magnetic energy $\mathcal {E}$. Their expressions are

(2.12)\begin{gather} N(t) = \int\,\textrm{d}A\,\textrm{d}\psi\ f_\textrm{tot}(\psi,A,t); \end{gather}

(2.13)\begin{gather} \langle{A}\rangle (t) = \frac{A_\textrm{tot}(t)}{N(t)}=\frac{1}{N(t)}\int\,\textrm{d}A\,\textrm{d}\psi\ f_\textrm{tot}(\psi,A,t)A; \end{gather}

(2.14)\begin{gather}\mathcal{E}(t) = \int\,\textrm{d}A\,\textrm{d}\psi\ f_\textrm{tot}(\psi,A,t) \psi^2. \end{gather}

The scalings for the time evolution of these three quantities are the most important predictions from our hierarchical model (Zhou et al. Reference Zhou, Bhat, Loureiro and Uzdensky2019), which we repeat here:

(2.15(a–c))\begin{equation} N(t)\sim t^{{-}1}; \quad \langle{A}\rangle (t) \sim t; \quad \mathcal{E}(t) \sim t^{{-}1}. \end{equation}

One of the main goals of this study is to test if the above scaling laws remain valid for a system consisting of a non-trivial distribution of islands, which we describe in detail in § 4.

Apart from the evolution of macroscopic quantities, we are also interested in the magnetic energy distribution over different length scales. This can be first represented by the magnetic energy density associated with different areas $A$ of islands, denoted as $U(A,t)$, which can be calculated from the distribution function $f(\psi,A,t)$ as

(2.16)\begin{equation} U(A,t) = \int\,\textrm{d}\psi\ f_\textrm{tot}(\psi,A,t) \psi^2, \end{equation}

and is related to the magnetic energy as $\mathcal {E}(t)=\int U(A,t)\,\textrm {d}A$. The conventional definition of magnetic energy spectrum in the wavenumber ($k$) space, denoted as $U(k,t)$, is related to $U(A,t)$. The wavenumber $k$, corresponding to a length scale, is related to the area of islands $A$ as

(2.17)\begin{equation} k \equiv \frac{2{\rm \pi}}{4\sqrt{A/{\rm \pi}}} \simeq A^{{-}1/2}. \end{equation}

The $k$ space can therefore be constructed from the $A$ space. The magnetic energy associated with each wavenumber and that associated with each value of area are related as $U(k)\,\textrm {d}k = U(A)\,\textrm {d}A$, where $A \sim k^{-2}$. Therefore, the magnetic energy spectrum $U(k,t)$ can be calculated from

(2.18)\begin{equation} U(k,t) = U(A,t)\frac{\textrm{d}A}{\textrm{d}k} \sim k^{{-}3} U(A,t). \end{equation}

Given an $A$-spectrum $U(A) \sim A^{-\alpha }$, we can obtain the corresponding $k$-spectrum $U(k) \sim k^{-\gamma }$, where $\gamma =3-2\alpha$, according to (2.18). We note that in some circumstances, the magnetic energy spectrum calculated using the definition (2.18) can provide more accurate and meaningful information on the energy distribution over scales than the standard one based on Fourier transform of real-space configurations of magnetic fields. For example, in Zhou et al. (Reference Zhou, Bhat, Loureiro and Uzdensky2019), a $k^{-2}$ magnetic energy spectrum was observed in 2-D RMHD simulations of island mergers, but the origin of this spectrum was, in fact, traced simply to the sharp magnetic field reversals at the thin current sheets between merging islands (Burgers Reference Burgers1948).

Finally, to visualise the island distribution function more easily, we calculate the 1-D distribution functions $F(A,t)$ (the distribution function over $A$) and $F(\psi,t)$ (the distribution function over $\psi$):

(2.19)\begin{equation} \left. \begin{aligned} F_A(A,t) & = \int\,\textrm{d}\psi f_\textrm{tot}(\psi,A,t),\\ F_\psi(\psi,t) & = \int\,\textrm{d}A f_\textrm{tot}(\psi,A,t). \end{aligned} \right\} \end{equation}

4.1. Delta-distribution initial condition

We first study the system with a delta-distribution initial condition, expressed as

(4.2)\begin{equation} f(\psi,A,t=0)=f_0\,\delta(\psi-\psi_\textrm{ini})\,\delta(A-A_\textrm{ini}). \end{equation}

That is, the system starts with an ensemble of identical magnetic islands with flux $\psi _\textrm {ini}$ and area $A_\textrm {ini}$ (the same initial condition as used in Zhou et al. Reference Zhou, Bhat, Loureiro and Uzdensky2019). In the following calculations, we set $\psi _\textrm {ini} = 5 B_0 L/6$ and area $A_\textrm {ini} = L^2/N_A$ (the smallest area we can resolve). We use $N_A = 256$ and $N_\psi = 4$ to resolve the system. In this and subsequent studies, we vary the Lundquist number $S_L$ between $10^3$ and $10^4$ and study the effect of non-local interactions by varying $R$ from $0$ to $1$. For this initial condition, we also vary the value of $f_0$ to set the island volume-filling fraction $V_\textrm {fill} \in \{0.125,0.25,0.5,1\}$. We focus mostly on the case $V_\textrm {fill}=1$ and study the effects of $S_L$ and $R$ with fixed $V_\textrm {fill}=1$.

This is the ideal set-up to test the scalings from the hierarchical model: total magnetic energy $\mathcal {E} \sim t^{-1}$, number of islands $N \sim t^{-1}$ and average area of islands $\langle A \rangle \sim t$ (Zhou et al. Reference Zhou, Bhat, Loureiro and Uzdensky2019). The evolution of these quantities is shown in figure 2. Figure 2(a,c,e) show the $R=0$ case (i.e. only identical islands interact) for $S_L = 10^3$, $10^4$. In this and subsequent figures, the unit of time is the global Alfvén time $L/v_{A0}$. Both $\mathcal {E}$ and $N$ show a $t^{-1}$ scaling to a good approximation, while $\langle A \rangle$ scales as $t$. This shows agreement with the hierarchical model described by Zhou et al. (Reference Zhou, Bhat, Loureiro and Uzdensky2019). There are discrete ‘steps’ that can be seen in the time evolution. The horizontal segments arise from the finite reconnection time, as the macroscopic quantities only change after merging is complete, while the steepness of the mostly vertical segments is controlled by the collision integral.

Figure 2. Delta-distribution initial conditions. (a,c,e) The evolution of macroscopic quantities ($\mathcal {E}$, $N$ and $ \langle{A}\rangle $) allowing only local interactions ($R=0$) between volume-filling islands ($V_\textrm {fill}=1$) with varying $S_L$. The corresponding scaling laws from the hierarchical model are shown for reference (dashed lines). (b,d,f) The evolution of macroscopic quantities with $S_L=10^4$, $V_\textrm {fill}=1$ and varying $R$. Allowing interactions between non-identical islands does not significantly change the evolution of macroscopic quantities (all the curves overlap each other such that the black dots are not visible, with only the $R = 1$ case showing minor deviations for $t \gtrsim 10$). (g,h,i) The evolution of macroscopic quantities with fixed $S_L=10^4$, $R=1$ and varying the volume-filling fraction of islands, $V_\textrm {fill}$. The evolution of the system does not depend strongly on $V_\textrm {fill}$.

The effect of allowing non-identical islands to interact is shown in figure 2(b,d,f) that compare runs with $R \in \{0,0.1,1\}$ and fixed $S_L=10^4$. There is little effect, as can be seen by the overlapping traces as $R$ is varied. This is because merging in this system is dominated by identical islands. This result shows that the assumption of only considering merging between identical islands used in Zhou et al. (Reference Zhou, Bhat, Loureiro and Uzdensky2019) is valid for the system it analyses.

The effect of the volume-filling fraction of islands is shown in figure 2(g,h,i), where $V_\textrm {fill}$ is varied from $0.125$ to $1$. We see that the overall evolution does not strongly depend on how compactly the islands are distributed in the system: all curves exhibit approximately the same slope. This is consistent with the fact that reconnection is the dominant dynamical process governing the evolution of the system. The islands’ mean-free path and volume-filling fraction only affect the convection time (i.e. the efficiency with which they pair up), while the merging process remains unchanged. Therefore, $V_\textrm {fill}$ should not be a critical parameter for system evolution within the regime where the convection time, $\lambda _\textrm {mfp}/v_A$, is much shorter than the merger time, $\tau _\textrm {rec} \sim \beta _\textrm {rec}^{-1} \sqrt {\langle{A}\rangle}/v_A$, which sets the condition $\beta _{\textrm {rec}} (\lambda _\textrm {mfp}/\sqrt {\langle{A}\rangle}) \ll 1$. Nevertheless, some small differences still occur between cases with different $V_\textrm {fill}$. With decreasing $V_\textrm {fill}$, the evolution of macroscopic quantities is slightly slower, which can be justified by the larger convection time (i.e. smaller pairing-up rate) of islands. The curve also becomes smoother with smaller $V_\textrm {fill}$, because islands have to travel to pair with one another, as opposed to just merging with immediately adjacent ones in the volume-filling case. This leads to different merger starting times for different pairs of islands, thus smoothing out the ‘step’ feature in the evolution curves.

The evolution of the distribution functions $F_A$ and $F_\psi$ for $S_L=10^4$ and $R=1$ is shown in figure 3. As the merging progresses, the peak of the area distribution moves from small to large scale (figure 3a). The peak value of $F_A$, which is approximately the number of islands $N$, is proportional to $\langle A\rangle ^{-1}$ as the total area is conserved. Because merging of different-sized islands is allowed, the distribution becomes wider with time, though the value at the peak remains 1–2 orders of magnitude larger than in the rest of the distribution, which confirms that the overall merging process is dominated by interactions between identical islands. Figure 3(b) shows the evolution of $F_\psi$. As there are only islands with one value of $\psi$ initially, the only change during evolution is the magnitude of the peak, which reflects the reduction in the number of islands. The spectra (not shown) are qualitatively similar to the distribution functions, which have one dominant peak.

Figure 3. The evolution of the 1-D distribution function in area, $F_A$ (a), and in flux, $F_\psi$ (b), for the run with $\delta$-function initial distribution, with $S_L=10^4$ and unconstrained interactions ($R=1$).

The results with the delta-distribution initial conditions provide a benchmark of the IKE model. We show that when the initial conditions and assumptions about merging are restricted to those of the model in Zhou et al. (Reference Zhou, Bhat, Loureiro and Uzdensky2019), the evolution of macroscopic quantities such as $\mathcal {E}$, $N$ and $\langle A \rangle$ agrees with the analytic predictions and with the direct numerical solution of the RMHD equations reported in that reference.

4.2. Gaussian-distribution initial condition

In a realistic system, magnetic islands will not be identical. For example, the population of large flux ropes (magnetic clouds) in the solar wind has a Gaussian distribution (Janvier, Démoulin & Dasso Reference Janvier, Démoulin and Dasso2014). We thus study an initial Gaussian distribution with a spread around its mean area and flux, expressed as

(4.3)\begin{equation} f(\psi,A,t=0)=f_0\exp({-(A-\bar{A})^2/(2\sigma_A^2)})\exp({-(\psi-\bar{\psi})^2/(2\sigma_{\psi}^2)}). \end{equation}

Our fiducial run has $\bar {A}= L^2/40 \approx 0.99$ and $\bar {\psi }=B_0L/2={\rm \pi}$, and we use $N_A = 256$ and $N_\psi = 4$ to resolve the system (a numerical convergence study is included in Appendix B).

The standard deviations are set as $\sigma _{A}=3 L^2/N_{A} \approx 0.3\bar {A} \approx 0.47$ and $\sigma _\psi =3B_0L/N_\psi \approx 3/2\bar {\psi } \approx 4.71$. Other runs (not shown in this paper) with different values of $\sigma _A$ and $\sigma _\psi$ show qualitatively similar results; in particular, when the values of $\sigma _A$ and $\sigma _\psi$ are small enough, the evolution of the system is similar to that with the delta-distribution initial condition reported in § 4.1. To study the difference between collisional reconnection, where $\beta _\textrm {rec} \simeq S^{-1/2}$, and the collisionless case, where $\beta _\textrm {rec} \simeq 0.1$, we perform an additional calculation in the collisionless regime with this same initial condition.

Figure 4 shows the time evolution of the macroscopic quantities $\mathcal {E}, N$ and $\langle A\rangle$ for different values of $S_L$ (panels a,c,e; the collisionless case is also included for comparison) and of $R$ (panels b,d,f). The time evolution of these quantities retains a power-law behaviour where the scaling of $\mathcal {E}$ and $N$ is $t^{-a}$ with $0.7 < a < 0.8$ and $\langle A\rangle$ scales as $t^a$. That is, an initial Gaussian distribution of islands leads to a somewhat slower evolution than the delta-distribution case. The ‘steps’ are still present during the time evolution in this case, as the distribution is narrow enough that the difference in merging time between the smallest and largest islands is not too large, and a large fraction of the islands complete the same number of mergers. The effect of changing $S_L$ or using the collisionless reconnection rate is shown in figure 4(a,c,e) (at fixed $R=1$). The reconnection rate, and hence the merging rate, increases going from $S_L = 10^4$ to $10^3$ to the collisionless regime. Merging starts earlier when the reconnection rate is higher. However, the slope of the power-law evolution is not strongly affected as the reconnection time scale is much longer than the Alfvénic time scale both in the SP and in the collisionless regimes.

Figure 4. Gaussian-distribution initial condition: (a,c,e) time evolution of macroscopic quantities ($\mathcal {E}$, $N$ and $ \langle{A}\rangle $) with $R=1$ and $S_L=10^3$ (black curve), $S_L=10^4$ (orange curve) and the collisionless case (green curve); (b,d,f) time evolution of macroscopic quantities at $S_L=10^4$ and varying $R$.

Figure 5. Gaussian initial distribution for $R=1$ and $S_L=10^4$: (a) time evolution of $F_A$; (b,c) time evolution of the energy distribution $U(A,t)$ and $U(k,t)$, respectively. Vertical lines indicate mean values of $A$ and $k$ at the corresponding moments of time. The solid (dotted) black curves in each panel indicate fittings of Gaussian distribution using $\sigma _A=0.47$ ($\sigma _A=0.6$) for $t=0$ ($t=10$).

Figure 6. Gaussian initial distribution with $R=1$ and the collisionless reconnection model: (a) time evolution of $F_A$; (b,c) time evolution of the energy distribution $U(A,t)$ and $U(k,t)$, respectively. Vertical lines indicate mean values of $A$ and $k$ at the corresponding moments of time.

Figure 4(b,d,f) show the effect of allowing for interactions between non-identical islands (i.e. changing $R$ for fixed $S=10^4$).

The evolution of macroscopic quantities is found to be (marginally) faster for systems with finite $R$ than for that with $R=0$, i.e. allowing for only identical-island interactions. The merging time of islands is proportional to $\sqrt {A_s}(A_s A_l)^{1/4} = (A_s^3 A_l)^{1/4}$ ((2.7) and (2.8)), where $s$ and $l$ refer to the smaller and larger islands in an interaction. Therefore, the merging time depends more strongly on the smaller island (with $A_s$). In the case of the Gaussian distribution, islands with sizes close to the mean area have the largest population and thus have a large effect on the overall dynamics. Compared with the merging time between two mean-sized islands, the merging time between a mean-sized island and a smaller island is shorter and dominated by the smaller island, while the merging time between a mean-sized island and a larger island is longer but only has a weak dependence on the area of the larger island. Therefore, allowing the interaction between non-identical islands essentially reduces the merging time. Additionally, more islands are allowed to merge at any given time. The combination of these effects makes the overall dynamics faster when $R$ is finite.

We also note that for the evolution of $N$ and $\langle A \rangle$, the curves for $R=0.01,0.03,0.1,1$ are essentially identical and are different from the $R=0$ curve. This suggests that the value of $R$ does not have a significant effect on the overall dynamics, as long as it is not zero. We conclude that highly non-local interactions (corresponding to mergers between islands differing a lot in size) are not dynamically important. This observation has implications for the discussion of decaying turbulence found in § 5.

The evolution of the distribution function $F_A$ and magnetic energy spectra are shown in figure 5 for the case $S_L = 10^4, R=1$ (given the very limited resolution in the $\psi$ coordinate, as mentioned earlier, there is not much value in studying $F_\psi$ and therefore it is not shown). At $t=10$, the initial distribution has decayed and a new peak grows at $A \approx 2 \bar {A}$ (panel a), as a result of the merging of islands with areas close to the first peak at $\bar {A}$ merging. The second peak of the distribution remains Gaussian, and has a somewhat larger width (comparing to the initial one), which we fit with $\sigma _A = 0.6$. After the merging of the initial islands, the peak of distribution continues to shift to larger area with a greater width. Owing to the longer merging time, the old peaks have not always decayed when the new peaks form (e.g. at $t=30, t=60$). We note that, at late times ($40< t<60$), the tail of $F_A$ at small $A$ remains (almost) unchanged and the mergers happen mainly between islands at the peak of $F_A$. This is because (i) the relatively fast mergers between (identical or similar) small islands are less frequent with larger values of $R$; and (ii) the probability for mergers between islands at the peak of $F_A$ is much higher than for mergers involving small islands because of their larger number density and cross sections.

The spectra evolution (panels b and c) show similar behaviour as the distributions.Footnote ⁶ We remark again that our definition of magnetic energy spectrum ((2.16) and (2.18)) is different from the conventional definition based on the Fourier transform of the magnetic field as a function of position in configuration space; the former straightforwardly presents magnetic energy distribution over islands with different scales, while the latter can be dominated by local features such as magnetic reversals at current sheets. The difference between the multi-peak spectra presented here and the $\sim k^{-2}$ spectrum in Zhou et al. (Reference Zhou, Bhat, Loureiro and Uzdensky2019) is caused by the different definitions of spectrum, without incompatibility. The collisionless case is shown in figure 6. The multiple peaks are still visible in the system both for $F_A$ and the spectra, but there is no clear transition from one Gaussian to another. Instead, multiple peaks are present because of the smaller merging time. This means that the first islands in a given generation to complete merging can start their second merger while other islands of that original generation are still merging, which gives rise to the multiple peaks.

4.3. Power-law initial condition

Another situation of interest is that of an initial power-law distribution of islands, which may be relevant to various astrophysical systems such as small-scale ($< 0.1$ AU) flux ropes in the solar wind (Janvier et al. Reference Janvier, Démoulin and Dasso2014) and magnetic structures in the downstream of a collisionless shock as a result of the Weibel instability (Katz, Keshet & Waxman Reference Katz, Keshet and Waxman2007). To reduce computational expense, we set up the initial islands with a power-law distribution over the size, but identical magnetic flux:

(4.4)\begin{equation} f(\psi,A,t=0)=f_0A^{-\alpha}\delta(\psi-\psi_\textrm{ini}). \end{equation}

The power-law exponent, $-\alpha$, is the same as that of its corresponding $A$-spectrum: $U(A) = \int \,\textrm {d}\psi f_0A^{-\alpha }\delta (\psi -\psi _\textrm {ini}) \psi ^2 =f_0 \psi _\textrm {ini}^2 A^{-\alpha }$. In our calculation, the power-law distribution is set up in the range $A_\textrm {min} < A < A_\textrm {max}$, where $A_\textrm {min} = L^2/N_A$ and $A_\textrm {max} = 5 L^2/N_A$. Initial conditions with $\alpha \in \{0.5,1,2,4\}$ are studied. We set $\psi _\textrm {ini} = (5/6)LB_0$ and resolve the system with $N_A=256$.

The evolution of various macroscopic quantities with $\alpha = 1$ in the initial condition is shown in figure 7. Panels (a,c,e) show cases with $S_L=10^3, 10^4$ and the collisionless case, for $R=1$. The decay of $\mathcal {E}$ and $N$ shows a $t^{-a}$ power-law scaling, and $\langle A \rangle$ increases as $t^{a}$, where $a \approx 0.9$ and does not strongly depend on the reconnection rate. This is only slightly different from our results for the delta and Gaussian initial distributions (reported in §§ 4.1 and 4.2). We note that in this case, the curves of evolution of macroscopic quantities are smooth: the discrete ‘steps’ that occur for the delta (figure 2) and Gaussian (figure 4) distributions are absent. This is because no dominant island size or merging time scale exists in the system owing to the wide spread of island sizes pertaining to a power-law distribution. Islands then merge over a wide range of time scales, which causes the evolution to appear smoother. The effect of changing $R$ is shown in panels (b,d,f) for $S_L=10^4$. As $R$ increases, the decay of $\mathcal {E}$ and $N$ becomes slightly steeper, as does the corresponding increase of $\langle A\rangle$, which indicates a faster overall merging process. The explanation is the same as for the Gaussian case: allowing non-local interactions in the collision integral increases the number of islands that can merge and increases the average merging rate.

Figure 7. Power-law initial distribution: (a,c,e) the evolution of macroscopic quantities ($\mathcal {E}$, $N$ and $ \langle {A}\rangle $) for unconstrained interaction ($R=1$) with $\alpha =1$ and $S_L=10^3,~10^4$ and the collisionless case; (b,d,f) the evolution of macroscopic quantities for $S_L=10^4$ and $\alpha =1$, with varying $R$.

The effect of varying the index of the initial power law from $\alpha = 0.5$ to $4$ on the evolution of the macroscopic quantities is shown in figure 8. The difference between the $\alpha =0.5$ and $\alpha =1.0$ runs is minor, which indicates that for shallow initial power-law distributions, the evolution of macroscopic quantities depends only weakly on the initial power-law slopes. As $\alpha$ increases, merging becomes faster, shown by somewhat steeper traces of $N$, $\mathcal {E}$ and $\langle A \rangle$. This can be understood by noting that the distribution function becomes progressively more similar to a delta-function distribution with an increasing $\alpha$, and so this steepening is consistent with the results of § 4.1 and should tend to ${\sim }t^{-1}$.

Figure 8. Power-law initial distribution. Time evolution of macroscopic quantities ($\mathcal {E}$, $N$ and $ \langle {A}\rangle $) for unconstrained interaction ($R=1$) and $S_L=10^4$, with varying $\alpha$.

The evolution of the energy spectra and area distribution function is shown in figure 9, for $\alpha =1$ (corresponding to an initial $k^{-\gamma }$ spectrum where $\gamma = 3-2 \alpha =1$), $R=1$ and $S_L=10^4$. The distribution function $F_\psi$ is not shown because it starts as a delta distribution and remains so, according to our merging rules. Interestingly, the area distribution $F_A$ (panel a) spontaneously forms a range of $A$ where a power-law distribution is maintained with $F_A \sim A^{-2}$, though the $-2$ index differs from the initial index. The formation of this $A^{-2}$ distribution is still present when the initial distribution has $\alpha = 0.5$ and $\alpha =2$,Footnote ⁷ which implies that this $A^{-2}$ distribution is an ‘attractor’ that the system tends to evolve to, independent of the initial conditions. Moreover, there is an extended envelope which covers the power-law regions of the curves at different times, showing the self-similar features of the system evolution. Consistent with $F_A$, the energy spectra also transition to a $U(A) \sim A^{-2}$ (and correspondingly $U(k) \sim k$) power law as the system evolves (panels b and c). Similarly, this behaviour is independent of the initial slope and the spectra tend to evolve to the $\sim A^{-2}$ ($\sim k$) fixed point. For the collisionless case (not shown), the evolution of $F_A$ and spectra show similar behaviours and evolve to power laws with the same exponents as their resistive-MHD counterparts. This self-similar evolution of $F_A$ and spectra can be demonstrated analytically by applying the hierarchical model of Zhou et al. (Reference Zhou, Bhat, Loureiro and Uzdensky2019) to a power-law distribution of islands under the assumption that coalescence events most often occur between similar size islands – see Appendix C. A derivation of the observed $k^1$ ($A^{-2}$) spectrum can be obtained based on the self-similar properties. However, we are not yet able to prove analytically that the $k^1$ ($A^{-2}$) scaling is an ‘attractor’ of the spectra and $F_A$ that the system is expected to evolve to. One final observation worth making is that, at later times, there is a peak at the centre of the distributions and spectra, which is likely a result of an incomplete merging of islands within that range of $A$, while the rest of the distribution in regions of smaller and larger $A$ still shows a $-2$ power law.

Figure 9. Power-law initial distribution for $R=1,\alpha =1$ and $S_L=10^4$: (a) time evolution of the area distribution function $F(A)$; (b,c) time evolution of spectra $U(A,t)$ and $U(k,t)$, respectively. The dashed lines show reference $A^{-2}$, $k^{1}$ power-laws in the respective plots.