3.1. Wave initialization

We consider the nonlinear interaction between two overlapping, perpendicularly polarized Alfvén waves that counter-propagate in a periodic 3-D domain along a uniform guide field $\boldsymbol {B}_{\boldsymbol {0}} = B_0 \hat {\boldsymbol {z}}$. The waves are described by positive constants representing the components of the wave vector that are perpendicular ($k_{\perp }$) and parallel ($k_{\parallel }$) to the guide field $\boldsymbol {B}_{\boldsymbol {0}}$. We employ a cubic box with $L_{\perp }=L_x=L_y=2{\rm \pi}$, $L_{\parallel } = L_z = 2{\rm \pi}$ and a uniform resolution of $(N_x, N_y, N_z) = (N_{\perp }, N_{\perp }, N_{\parallel })$ cells per initial wavelength. For a scale-free definition of characteristic (wave)lengths, we prescribe $k_{\parallel }=2{\rm \pi} /L_{\parallel }$ and $k_{\perp }=2{\rm \pi} /L_{\perp }$. Then, we examine the waves described by the wave vectors $\boldsymbol {k}^{+}_1=k_{\perp } {\hat {\boldsymbol {x}}} - k_{\parallel } {\hat {\boldsymbol {z}}}$, from now on represented as $(1,0,-1)$ and $\boldsymbol {k}^{-}_2=k_{\perp } \hat {\boldsymbol {y}} + k_{\parallel } \hat {\boldsymbol {z}}$, or $(0,1,1)$. The magnetic field is initialized through a vector potential $\boldsymbol {A} = (-B_0 y, 0, \delta B_{\perp }[\sin ({k_\perp }x+{k_\parallel }z) + \sin ({k_\perp }y-{k_\parallel }z)])$, which represents the initial counter-propagating Alfvén waves with frequency $\omega _0 = k_{\parallel } v_A$ and $\hat {B} = {\delta B_{\perp }}/{B_0}$. The electric field is initialized as $\boldsymbol {E} = (v_A B_y, v_A B_x,0)$ such that the velocity is equal to the drift velocity $\boldsymbol {v} = \boldsymbol {E} \times \boldsymbol {B}/B^2$. We note that the overlapping Alfvén waves, in contrast to a single Alfvén wave, are not an exact force-free MHD equilibrium owing to a small second-order violation $\boldsymbol {E}_{\mp } \boldsymbol {\cdot } \boldsymbol {B}_{\pm } \neq 0$ between the fields of the two waves. These perturbations are instantaneously algebraically cut in ET-FFE and damped on a short timescale in BHAC. Initially, we set a gas to magnetic pressure ratio of $\beta =2p/B_0^2=0.02$ and set the magnetization $\sigma \in [10^{-2};10^{-1};1;10;100]$ (corresponding to $\sigma _\textrm {hot} \in [10^{-2};10^{-1};1;7;20]$), where we vary the density $\rho$ and keep a constant guide field with $B_0=1$. We fix the pressure to $p=0.01$ and employ an adiabatic index $\hat {\gamma }=4/3$ for an ideal relativistic gas. In the force-free MHD case, $\beta \rightarrow 0$, $\sigma \rightarrow \infty$ and $v_A \rightarrow 1$.

Through the three-wave interaction, turbulence can transfer magnetic energy anisotropically from large to small scales. The two initial waves interact to form an inherently nonlinear, purely magnetic mode, physically representing a shear in the magnetic field with wave vector $\boldsymbol {k}^{0}_2 = k_{\perp }(\hat {\boldsymbol {x}}+ \hat {\boldsymbol {y}})$, or $(1,1,0)$, which does not grow secularly in time but oscillates at twice the frequency of the primary waves (see Paper I). The interaction between this secondary mode and the primary Alfvén modes nonlinearly generates two tertiary Alfvén waves whose energy grows secularly in time. The $\boldsymbol {{k_1^{+}}}$ wave transfers energy to an Alfvén mode with $\boldsymbol {k}_3^{+} = 2k_{\perp } \boldsymbol {\hat {x}} + k_{\perp } \boldsymbol {\hat {y}} - k_{\parallel } \boldsymbol {\hat {z}}$, or $(2,1,-1)$, and the $\boldsymbol {{k_1^{-}}}$ wave to $\boldsymbol {k}_3^{-} = k_{\perp } \boldsymbol {\hat {x}} + 2k_{\perp } \boldsymbol {\hat {y}} + k_{\parallel } \boldsymbol {\hat {z}}$, or $(1,2,1)$.

To capture all physical behaviour including nonlinear interactions at play in magnetized turbulence, we solve the full 3-D set of equations. The strength of the nonlinearity is characterized by $\chi = k_{\perp } v_{\perp } / (k_{\parallel } v_A)$, where $k_{\perp } = k_{\parallel }$ in the set-up presented in this section, and $v_{\perp } / v_A \sim \delta B_{\perp } / B_0=10^{-1}$. This gives $\chi = 10^{-1}$, which results in a nonlinear time of $\chi ^{-2} = 100$ wave-crossings $2{\rm \pi} / \omega _0$. The limit of strong turbulence occurs when the nonlinear energy transfer timescale is comparable to the wave period, as to say, when $\chi \gtrsim 1$ (Goldreich & Sridhar Reference Goldreich and Sridhar1995). The nonlinearities governing these dominant three-wave interactions require variations in both directions perpendicular to the equilibrium magnetic field (Howes & Nielson Reference Howes and Nielson2013; Howes Reference Howes2014, Paper I). In the two-dimensional (2-D) limit perpendicular to the uniform field, the dominant nonlinearities governing the three-wave interactions are retained, yet the linear physics of the anisotropic cascade is absent: the linear interactions representing the propagation of the Alfvén waves along the magnetic field are only non-zero when the parallel wavenumber $k_{\parallel }$ is non-zero, which requires a (third) field-parallel dimension.

3.2. Current sheet formation

In this section, we analyse the formation of the sheets, their thickness, and the dynamics until the moment they break-up and magnetic energy is quickly dissipated. As shown in Paper I, the energy transfer to smaller scales is mediated by self-consistently generated $k_\parallel =0$ modes in the $\chi \ll 1$ limit, specifically, at least until the nonlinear time of $\sim \chi ^{-2} = 100$ wave-crossings. If the nonlinearly generated Alfvén waves increase in power, they can generate coherent current sheet-like structures owing to constructive interference with the primary Alfvén modes. Although in the ideal limit, current sheets are Dirac delta-functions, at finite dissipation, they can be represented by many fewer modes such as those in Howes (Reference Howes2016), who showed that the formation of current sheets results from interference between just five complex Fourier modes. To study whether current sheets can form and become thin enough to provide an efficient dissipation channel, we evolve the system for 200 wave-crossings, beyond the nonlinear time of $\sim \chi ^{-2} = 100$ wave-crossings.

Current sheets emerge once the wave interaction reaches a nonlinear state. To track the dynamics of an emerging current sheet, we show the evolution of the normalized in-plane current density $\hat {J}_z = [\boldsymbol {\nabla } \times \boldsymbol {B}]_z / (k_{\parallel }B_0)$Footnote ³ perpendicular to the equilibrium field $\boldsymbol {B} = B_0\hat {\boldsymbol {z}}$ at $z={\rm \pi}$ between 100 and 150 wave-crossings in figure 1, for the ET-FFE run with $512^3$ grid-points. We find that the results from ET-FFE ($512^3$) are comparable to the high-resolution runs conducted with BHAC ($1024^3$). The high-order reconstruction methods employed in ET-FFE capture the smallest scales with very good accuracy (see also § 3.3 and Mahlmann et al. Reference Mahlmann, Aloy, Mewes and Cerdá-Durán2021b). Current sheets are characterized by in-plane magnetic field line reversals (overplotted as black lines) and are clearly distinguishable by a strong current density (red and blue colours). At the interface of interacting eddies, current sheets are compressed into thin layers, as we examine in detail in figure 2. In figure 1, such regions are indicated by an anti-parallel in-plane magnetic field, for example at $x=0.9{\rm \pi}, y=0.9{\rm \pi}$ at $t=125\ [2{\rm \pi} /\omega _0]$. We also detect regions of non-zero current density which are not associated with current sheets owing to the lack of anti-parallel magnetic field, i.e. there is no in-plane magnetic null, as in figure 1 at $x=0.7{\rm \pi}, y=0.7{\rm \pi}$ at $t=110\ [2{\rm \pi} /\omega _0]$.

Figure 1. Normalized out-of-plane current density $\hat {J}_z = [\boldsymbol {\nabla } \times \boldsymbol {B}]_z / (k_{\perp }B_0)$ in the $(x,y)$-plane perpendicular to the guide field showing the formation ($t \lesssim 100\ [2{\rm \pi} /\omega _0]$), thinning ($t \sim 100-120 \ [2{\rm \pi} /\omega _0]$) and break-up of current sheets ($t \sim 120\text {--}130 \ [2{\rm \pi} /\omega _0]$) into developed turbulence ($t \gtrsim 130\text {--}150 \ [2{\rm \pi} /\omega _0]$) with ET-FFE. In-plane magnetic field lines are overplotted. We note that the evolution of the current density is similar for a run at $1024^3$ cells with BHAC. Figure 4 presents a 3-D rendering of the evolving current sheet enclosed by the dashed green rectangles. An animated version of this figure is provided in Media Supplement A (2021).

Figure 2. Current sheets form between elongated eddies. Zoomed-in view of the out-of-plane current density $\hat {J}_z = [\boldsymbol {\nabla } \times \boldsymbol {B}]_z / (k_{\perp }B_0)$ at 125 wave-crossings. Contours of the out-of-plane vector potential $A_z$ (where $\boldsymbol {B} = \boldsymbol {\nabla } \times \boldsymbol {A}$), which represent in-plane magnetic field lines, are overplotted in green, showing that regions of large current density originate from compressed magnetic field (direction indicated by arrows). We select the strongest current sheet in the domain of a BHAC $1024^3$ simulation (a) and an ET-FFE $512^3$ set-up (b). To indicate the respective applied resolutions, we show a fraction of the numerical grid stacked with blocks, where each block consists of $8\times 8$ cells. The current sheets (in bright red) are covered along their thickness $\delta \approx 0.01 \times 2{\rm \pi} \approx 0.06$ by at least 8 cells in both cases. A one-dimensional (1-D) outline of field quantities across the centre of each current sheet (dashed black lines) is presented in figure 3.

The current sheets in this system can be characterized by the aspect ratio in the perpendicular plane between their thickness, $\delta$, and width, $w$. We determine $\delta$ as the full-width at half maximum of the out-of-plane current density $\hat {J}_z$. The sheet becomes thinner between 100 and 127 crossings, until it breaks up. The width of the sheet $w$, namely, its extension along the interface of two eddies, is determined by the perpendicular wavelength $\lambda _{\perp }=2{\rm \pi}$ of the initial waves in the box. It remains approximately constant $w\approx 2{\rm \pi}$ from its formation at 100 crossings until the sheet breaks up at 127 crossings (see also figure 3). The aspect ratio decreases from $\delta /w \gtrsim 0.025$ (ET-FFE) at 100 wave-crossings onward to $\delta /w \approx 0.01$ (ET-FFE/BHAC) at 125 wave-crossings (see also figure 3(e), indicated by the black dashed line). The reconnection rate is expected to scale with the aspect ratio, $v_\textrm {rec}/c \sim \delta /w$, which suggests that the sheet breaks up in the linear regime of the instability (Biskamp Reference Biskamp2000; Ni et al. Reference Ni, Germaschewski, Huang, Sullivan, Yang and Bhattacharjee2010), and in the nonlinear regime, the asymptotic ‘fast reconnection rate’ (Bhattacharjee et al. Reference Bhattacharjee, Huang, Yang and Rogers2009; Uzdensky, Loureiro & Schekochihin Reference Uzdensky, Loureiro and Schekochihin2010; Huang & Bhattacharjee Reference Huang and Bhattacharjee2016) of $0.01$ might be reached. The life-time of the sheets is approximately 25 wave-crossings, and hence it is maintained on timescales (and length scales) of the initial large-scale waves (cf. Howes, McCubbin & Klein Reference Howes, McCubbin and Klein2018). Turbulent fluctuations initiate inside the current sheets after 125 wave-crossings, where energy is quickly dissipated at the grid scale, which is smaller than the current sheet thickness, $\Delta x / \delta \lesssim 0.2$.

Figure 3. Current sheets form similarly across different codes and resolutions. Profiles of normalized magnetic field $\hat {B}_i = \delta B_i/B_0$ (a,b) in minimum variance coordinates (cf. Sonnerup & Cahill Reference Sonnerup and Cahill1967; Howes Reference Howes2016), and out-of-plane current density $\hat {J}_z = [\boldsymbol {\nabla } \times \boldsymbol {B}]_z / (k_{\perp }B_0)$ (c,d) at $t=125\ [2{\rm \pi} /\omega _0]$. The 1-D outlines are extracted across the current sheet (see dashed black lines in figure 2), where the coordinate $s$ is centred at the respective structure. (e) Aspect ratio $\delta /w$ of the current sheet measured as the full-width half maximum of the current density normalized by its approximately constant width $w\approx 2{\rm \pi}$ versus wave-crossings for both BHAC and ET-FFE. The current sheet forms around the nonlinear time $t=100\ [2{\rm \pi} /\omega _0]$ and the thinning process halts at around 125 wave-crossings at $\delta /w\approx 0.01$, which suggests that the sheet breaks up once an asymptotic reconnection rate $v_\textrm {rec}/c \sim \delta /w \approx 0.01$ (indicated by a vertical dashed line) is reached.

The 3-D structure of the current sheet is lucidly illustrated in figure 4. We provide a volume rendering of $\hat {J}_z$ and observe a clear coherent structure that breaks after 127 wave-crossings. Initially, the current structure is, indeed, sheet-like in three dimensions. In the direction parallel to the guide field, the characteristic scale is the current sheet's length, $L$. We measure the length of the sheet along the $z$-axis to be $L\approx 2{\rm \pi}$, which is determined by the parallel wavelength $\lambda _{\parallel }=2{\rm \pi}$ along the equilibrium magnetic field. At the same time, it becomes visually evident that the extension in the plane perpendicular to the guide field is less than the box diagonal, but rather corresponds to $w\approx 2{\rm \pi}$ (see above). At $t=125\ [2{\rm \pi} /\omega _0]$, the current sheet is already very thin (see also figure 1). In the very short time interval of $\Delta T\approx 5\ [2{\rm \pi} /\omega _0]$, the current sheet ripples and breaks up into turbulent structures in the plane perpendicular to the guide field. During this time, we identified that non-ideal electric fields build up in the localized region of the current sheet. Such fields are subsequently damped to maintain a (globally) force-free field configuration.

Figure 4. Volume rendering of the normalized out-of-plane current density $\hat {J}_z = [\boldsymbol {\nabla } \times \boldsymbol {B}]_z / (k_{\perp }B_0)$, which shows the moment of the break-up of the current sheet between ${t=125\ [2{\rm \pi} /\omega _0]}$ and ${t=130\ [2{\rm \pi} /\omega _0]}$, within the reduced volume $x\times y\times z=[2{\rm \pi} /5,6{\rm \pi} /5]\times [2{\rm \pi} /5,6{\rm \pi} /5]\times [0,2{\rm \pi} ]$ at a resolution of $512^3$. An animated version of this figure is provided in Media Supplement B (2021).

3.3. Comparison of current sheet characteristics in ET-FFE and BHAC

The emerging current sheet is analysed in detail in figures 2 and 3 and compared between the employed numerical schemes. We illustrate the current sheet properties and their numerical resolution in figure 2. There, we plot numerical blocks of $8\times 8$ cells in a zoom-in on the strongest currents in the domain, which shows that the current sheet is indeed resolved by more than one block (approximately 11 cells) at the break-up time when the sheet is at its thinnest point, in the case of BHAC with a total of $1024^3$ cells, and by slightly less than one block (approximately 8 cells) in the case of ET-FFE with a total of $512^3$ cells. While such resolution arguments do not imply a convergence of dissipative dynamics in the force-free current sheet, they correspond to the minimum resolution needed to capture their formation in the magnetically dominated regime (as was established by Mahlmann et al. Reference Mahlmann, Aloy, Mewes and Cerdá-Durán2021b, namely, 5 to 10 cells per current sheet width). A notable characteristic of the examined sheet can be drawn from the overplotted contours of the out-of-plane vector potential $A_z$ (green lines): one can see that the regions of high current density originate from magnetic field compression at the edges of the elongated eddies.

The current sheet formation and evolution are very similar for the different codes and resolutions. In figure 3(a,b), we show the normalized magnetic field $\hat {B}_i = \delta B_i/B_0$ and (c,d) show the normalized out-of-plane current density $\hat {J}_z = [\boldsymbol {\nabla } \times \boldsymbol {B}]_z / (k_{\parallel }B_0)$ along a cut perpendicular to the current sheet presented in figure 2 just before it breaks at the time $t = 125\ [2{\rm \pi} /\omega _0]$. In this analysis, we choose minimum variance coordinates (cf. Sonnerup & Cahill Reference Sonnerup and Cahill1967; Howes Reference Howes2016) for comparability. These coordinates combine the guide field direction $\hat {\boldsymbol {m}}=\hat {\boldsymbol {z}}$, as well as the in-plane directions $\hat {\boldsymbol {n}}=(\hat {\boldsymbol {x}}+\hat {\boldsymbol {y}})/\sqrt {2}$ and $\hat {\boldsymbol {l}}=(\hat {\boldsymbol {x}}-\hat {\boldsymbol {y}})/\sqrt {2}$. For both BHAC and ET-FFE, we use the location of the strongest current in the domain for the analysis (indicated by the black dashed line in figure 2). The current sheet shows a typical profile with a peak in the current density and magnetic field reversals of the reconnecting in-plane components $\hat {B}_x$ and $\hat {B}_y$ (i.e. the perpendicular component $\hat {B}_l$). The gradient of the guide field pressure (red line) is non-zero in the current sheet, which compensates for gas pressure that is absent in the force-free limit of MHD, and effectively working to sustain the evolving force-free current sheet by a locally balanced magnetic pressure at the location of the field reversal (similar to the stationary force-free current sheets analysed by Komissarov, Barkov & Lyutikov Reference Komissarov, Barkov and Lyutikov2007; Del Zanna et al. Reference Del Zanna, Papini, Landi, Bugli and Bucciantini2016; Mahlmann et al. Reference Mahlmann, Aloy, Mewes and Cerdá-Durán2021b).

The process of current sheet formation and thinning is independent of the manner in which force-free violations are managed (see also Appendix B) and is likely dominated by the field compression driven by mode interactions. The thickness $\delta = 0.01w$ after 125 crossings is approximately 10 times larger than the smallest resolvable perpendicular scale in the BHAC simulation, $\sim 2{\rm \pi} /1024$ and five times larger than the smallest resolvable scale in the ET-FFE simulation, $\sim 2{\rm \pi} /512$ (figure 3e). Additionally, we measure the thickness of the current sheet in time for progressively increasing resolutions in the right panel, which shows that the thinning rate and aspect ratio in ET-FFE converges for $256^3$ and $512^3$ grid points, while for BHAC, it converges for $512^3$ and $1024^3$ grid points. The thinning rate and aspect ratio also converge between the highest resolution results of the two different algorithms. We demonstrate that the current sheet's thickness is larger than a single numerical cell by measuring the number of cells across the sheet's thickness at 125 wave-crossings just before the sheet breaks up. The observation that the current sheet formation, thinning as well as asymptotic thickness and time of break-up coincide between runs in BHAC and ET-FFE (with a largely different order of convergence and treatment of the force-free violations) is worth mentioning in this context. Small local differences aside, the comparison of results from the two numerical schemes in figures 2 and 3 shows that the field structure is remarkably similar in both cases. The global dynamics of field compressions is captured accurately in the force-free limit of MHD. In Appendix B, we compare the appearance of non-ideal electric fields in both methods, and conclude that despite numerical differences, the results show a remarkable level of similarity both in magnitude and location of the force-free violations.

In summary, the following features let us conclude that the identified structures are 3-D current sheets: (a) strong and localized ordered current structures; (b) field reversals with a compensating magnetic pressure; (c) compression and subsequent thinning of the current structure; and (d) non-ideal electric fields around magnetic null-lines (see Appendix B) in the plane perpendicular to the guide field.

3.4. Weak turbulence development

In this section, we investigate how energy cascades and dissipates at small scales. In force-free MHD, the total electromagnetic energy is conserved in the domain until the turbulent cascade reaches the numerical grid scale. By measuring the total electromagnetic wave-energy $E_\textrm {EM} = (B-B_0)^2+E^2$, we determine the dissipation rate depending on the resolution and the explicit resistivity, $\eta$ (in BHAC runs). For a resolved simulation, i.e. where the current sheet thickness is resolved by multiple cells, the energy $E_\textrm {EM}$ should dissipate as $\textrm {d}E_\textrm {EM}/\textrm {d}t \sim \boldsymbol {E} \boldsymbol {\cdot } \boldsymbol {J} \sim \eta J^2 \sim \eta |\boldsymbol {\nabla } \times \boldsymbol {B}|^2$, where $\boldsymbol {\nabla } \times \boldsymbol {B} \sim B_k/\lambda$ and $\lambda =2{\rm \pi}$ is the wavelength, which implies that the typical dissipation time is $\tau \sim \lambda ^2\eta ^{-1}$, such that we can approximate the total energy as $E_\textrm {EM}(t) \approx E_{\textrm {EM},t=0}(1-\eta ({2{\rm \pi} t}/{\omega _0}))$. Figure 5(a) shows the estimated dissipation rate (dash–dotted lines) for force-free MHD runs with $\eta =10^{-4}$, $5\times 10^{-5}$ and $10^{-6}$. We conclude that for $\eta \leqslant 10^{-6}$ and resolutions of $\gtrsim 512^3$, Ohmic dissipation has a negligible effect on the initial 130 wave-crossings before the current sheets break-up. The energy dissipation rate agrees very well between ET-FFE ($512^3$ cells) and BHAC ($\eta =10^{-6}$) for a resolution of $1024^3$ cells. This result demonstrates that the formation and the thinning process of the current sheets, which occur between 100 and 125 wave-crossings, neither result in a significant energy decay nor in a violation of the ideal force-free constraints that would reduce the electromagnetic energy density (see also Appendix B). After approximately 125 wave-crossings, a steep decay is observed, which corresponds to the break-up of current sheets in figure 1, and energy is quickly dissipated at the grid-scale. We further conclude that the energy decay is a result of the cascade initiated by the wave interactions by showing that a single non-interacting Alfvén wave (green dashed line in figure 5a) shows no significant energy loss over 200 wave-crossings as indicated in figure 5(a). In figure 5(b), we show the MHD runs (BHAC) at $\sigma =1$, $10$ and $100$, which correspond to smaller initial total electromagnetic energy, showing a similar evolution with a larger variability owing to the transfer of electromagnetic energy to fluid (kinetic and thermal) energy.

Figure 5. Evolution of the (normalized) total wave energy $\hat {E}_\textrm {EM} = E_\textrm {EM} / E_{\textrm {EM},t=0}$, with $E_\textrm {EM}=(B-B_0)^2 + E^2$ for the overlapping Alfvén waves in all force-free MHD cases at different resolutions (a). We present ideal MHD cases and a force-free reference case ((b), BHAC) normalized to the initial electromagnetic wave energy of the force-free set-up for 200 wave-crossings, with initial wave-amplitude $\chi =10^{-1}$ and $N_{\perp }=N_{\parallel } = 512$ grid-cells. In (a), we overplot results for BHAC (blue, red and magenta lines for different resistivities; dotted, solid and thick lines for different resolutions) at resolutions $N_{\perp }=N_{\parallel } = 256$, $512$ and $1024$ as well as ET-FFE (black lines) at $256$ and $512$ grid cells to show that the energy evolution and current sheet formation converges. We also show the (conserved) energy density for a single Alfvén wave at 256 cells/wavelength (green dotted line), which does not undergo a cascade (a). The theoretical expectation of the energy dissipation for a fixed resistivity is indicated by grey dash–dotted lines, with the respective value of $\eta$ shown below the respective line.

We analyse the evolution of the spectral magnetic energy $E_B(k)\,\text {d}k=\mathop {\sum }\nolimits _{{\boldsymbol {k}}\in \text {d}k} {\boldsymbol {B}}_{\boldsymbol {k}}\boldsymbol {\cdot }\boldsymbol {B}^{*}_ {\boldsymbol {k}}$ for the resolvedFootnote ⁴ ET-FFE is run at a resolution of $512^3$ cells in figure 6. We observe a magnetic energy spectrum with a $E_{B_\perp } (k_{\perp })\propto k_{\perp }^{-2}$ power law forming after 100 wave-crossings (green line in panel (a,c)) as a result of three-wave interactions. The power spectrum indicates that the turbulence in our simulations remains in the the weak regime $k_\parallel B_0 \gg k_\perp \delta B$ (Ng & Bhattacharjee Reference Ng and Bhattacharjee1997; Bhattacharjee & Ng Reference Bhattacharjee and Ng2001; Kuznetsov Reference Kuznetsov2001; TenBarge et al. Reference TenBarge, Ripperda, Chernoglazov, Bhattacharjee, Mahlmann, Most, Juno, Yuan and Philippov2021). The transition from weak to strong turbulence is expected to occur after at least $\sim \chi ^{-2}=100$ wave-crossings. When the critical balance condition $k_{\parallel } B_0 \sim k_\perp \delta B$ becomes satisfied, such that the turbulence is in the strong regime, a $E_{B_\perp } (k_{\perp })\propto k_{\perp }^{-3/2}$ energy spectrum is expected to form (Perez & Boldyrev Reference Perez and Boldyrev2008; Schekochihin, Nazarenko & Yousef Reference Schekochihin, Nazarenko and Yousef2012; Verdini & Grappin Reference Verdini and Grappin2012). Owing to limited resolution and, hence, a limited inertial range, the wave energy starts to dissipate at the grid scale after ${\sim }120$ wave-crossings in our simulations such that the transition to strong turbulence is not captured here.

Figure 6. Evolution of spectra (a,b) and non-linearity parameter $\chi$ (c,d) during the continuous wave interaction in the force-free limit. (a) $\hat {E}_{B_\perp }$ spectra at every 20 wave-crossings with wave-amplitude $\chi =10^{-1}$ and resolution $N_{\perp }=N_{\parallel } = 512$ grid-cells until 200 wave-crossings in ET-FFE. (b) Spectra at every wave-crossing between 120 and 130 wave-crossings. (c,d) Same analysis in time for the non-linearity parameter $\chi$ (reduced to its maximum value along slices of equal $k_{\parallel }$). There is a clear transition between 120 and 130 wave-crossings, when the current sheets thin and break-up after 127 wave-crossings, and an extended $\hat {E}_{B_\perp } \propto k_{\perp }^{-2}$ spectrum develops. We note that the spectral evolution in the ET-FFE run at $512^3$ resolution is similar to BHAC at $1024^3$ resolution.

When the energy cascades to the highest $k$ (or thinnest structures) after 100 wave-crossings, pronounced thin current sheets of high-amplitude current density form (see figure 1). Following Howes (Reference Howes2016), the emerging current sheets can be described by as few as the five lowest order modes. Only when smaller scale structures form in the thinning and breaking of current sheets, energy is transferred to higher order modes. This effect can be observed in the shift in the spectral evolution between 100 wave-crossings and the break-up point at 127 wave-crossings (see figure 6(b,d) for a zoom into the interval between 120 and 130 wave-crossings). At this point, the inertial range extends for almost two decades of $k_{\perp }$. Note that a current sheet itself results in a spectral index of $-2$ by assuming the Fourier transform of a (near)-step function (Burgers Reference Burgers1948). When the current sheets break-up at 127 wave-crossings, a maximally developed spectrum has formed as any thinner structures cannot exist on the numerical grid. Once the sheets break-up, the smallest scale (highest $k$) is reached by the cascade and the magnetic energy starts to dissipate at the grid-scale (roughly corresponding to $k_\perp >100$). The spectral energy decays from here onward (orange and red lines).

The growth of turbulence at smaller scales and the subsequent decay is illustrated in the analysis of the $\chi$ parameter in figure 6(c,d). The steady-state expectation of $\chi (k_{\perp })$ can be derived by assuming that the parallel cascade is negligible, and the weak turbulence spectrum shows that $\delta B_{\perp }^2 / k_{\perp } \sim k_{\perp }^{-2}$ and hence $\chi = k_{\perp } v_{\perp } / (k_{\parallel } v_A) \sim \delta B_{\perp } k_{\perp } / (k_{\parallel } B_0) \sim k_{\perp }^{1/2} / k_{\parallel } \sim k_{\perp }^{1/2}$, indicated by the dashed line. The nonlinearity parameter clearly indicates that the turbulence remains in the weak regime $\chi < 1$ and develops towards $\chi \sim k_{\perp }^{1/2}$. The increasing nonlinearity parameter with $k_{\perp }$ and the cascade becoming more anisotropic suggests that, for higher resolutions, the transition to strong turbulence would be inevitable (Meyrand, Galtier & Kiyani Reference Meyrand, Galtier and Kiyani2016).

3.5. Local weak turbulence properties

In this subsection, we investigate the local properties of the turbulence to determine the local anisotropy of the spectrum. The Fourier analysis employed for the computation of the spectra cannot take into account the local properties of turbulence. It implicitly assumes the local mean magnetic field remains parallel to the initial field $\boldsymbol {B}_{\boldsymbol {0}}$. However, each small scale eddy responds to the combined, local contribution of larger scale eddies. The weak turbulence limit requires many collisions of counter-propagating waves to substantially deform the initial eddy, and an Alfvénic eddy propagates along the local, mean magnetic field. Accounting for the mismatch of local and global guide field requires one to define locally parallel and perpendicular cascades (Cho & Vishniac Reference Cho and Vishniac2000; Maron & Goldreich Reference Maron and Goldreich2001). To mitigate the effect of mixing parallel and perpendicular directions, we employ second-order structure functions:

(3.1)\begin{equation} \delta B_i^2 ({\boldsymbol{l}}) = \langle \left\lvert B_i({\boldsymbol{r}}+{\boldsymbol{l}})-B_i({\boldsymbol{r}}) \right\lvert^2\rangle, \end{equation}

where $l$ is a separation vector and angular brackets denote averages over all ${\boldsymbol {r}}$. These functions contain the information about the local spectrum $f({\boldsymbol {r}_{\boldsymbol {0}}}, {\boldsymbol {l}})$ at a given point $\boldsymbol {r}_{\boldsymbol {0}}$, where $l\sim k^{-1}$. The local guide field for two points of an eddy of size $l$ can be defined as

(3.2)\begin{equation} {\boldsymbol{B}}_\textrm{guide}({\boldsymbol{r}},{\boldsymbol{l}}) = \frac{{\boldsymbol{B}}({\boldsymbol{r}}+{\boldsymbol{l}})+ {\boldsymbol{B}}({\boldsymbol{r}})}{2}. \end{equation}

The shape of an eddy is given as $l_\parallel =\boldsymbol {l} \boldsymbol {\cdot } \boldsymbol {B}_\textrm {guide}/|{\boldsymbol {B}}_\textrm {guide}|$, $l_{\perp } = (l^2-l^2_\parallel )^{1/2}$, assuming that all eddies are isotropic in the field-perpendicular plane and neglecting effects of dynamical alignment and intermittency (Boldyrev Reference Boldyrev2005,  Reference Boldyrev2006; Chandran, Schekochihin & Mallet Reference Chandran, Schekochihin and Mallet2015; Mallet, Schekochihin & Chandran Reference Mallet, Schekochihin and Chandran2015). The construction of (3.1) requires calculation of a six-dimensional (6-D) integral, and we employ a Monte Carlo method with a total of $N=N_0 n_x n_y n_z$ randomly chosen pairs of vectors to uniformly cover both the space of positions of eddies $\boldsymbol {r}$ and their sizes $\boldsymbol {l}$. Limits for the point-separation vectors $l_i$ are chosen to be $(-n_i/2,n_i/2)$ to include the effects of the periodic boundary conditions in direction $i$. We set $N_0=30$ to decrease the shot-noise and confirm convergence of the results for $N_0=2,5,20$.

The 2-D structure functions for an ET-FFE run at a resolution of $512^3$ cells are shown in figure 7, at time $t=120\ [2{\rm \pi} /\omega _0]$, when thin current sheets have formed and are not yet broken up (see figure 1). All vectors ${\boldsymbol {l}}$ are split into $128$ bins in both parallel and perpendicular directions. The white regions at larger $l_\parallel$ and $l_{\perp }$ do not contain any points as all separations larger than $l_i/2$ are prohibited owing to the periodic grid. The shape of the structure function is defined by the winding of magnetic field lines. A bright over-scaled feature at $l_\parallel \sim 0.5 L_\parallel$ and $l_\perp < 0.5 L_\perp$ is associated with the thinning current sheets: its amplitude grows while current sheets are being formed and disappears when current sheets break-up at $t\approx 128$. This feature appears both in the BHAC and ET-FFE results and is independent of how force-free violations are handled (see Appendix B).

Figure 7. Two-dimensional structure functions for (a) $\delta \hat {B}_\parallel$ and (b) $\delta \hat {B}_{\perp }$ with respect to the local guide field, at time $t=120\ [2{\rm \pi} /\omega _0]$ before the current sheets break-up. Here, $l_\parallel$ and $l_\perp$ are defined with respect to the local guide field. White regions at large $l$ do not have any points owing to periodic boundary conditions. The bright feature at $l_\parallel /L_\parallel \approx 0.5$ corresponds to the current sheets.

The natural ability of this approach to consider local guide field variation allows us to measure a local angular anisotropy of the spectrum. To perform such a calculation, we split all point-separating vectors into 18 angular bins of $5^\circ$ extension, measured with respect to the field-parallel direction ($0^\circ$ is parallel to the guide field and $90^\circ$ is perpendicular to it). In each bin, the intermediate range of $|{\boldsymbol {l}}|$ is approximated by a power-law function $\propto l^g$. The spectral index $\alpha$ is found using the relation $\alpha = -(g+1)$ (Monin & Yaglom Reference Monin and Yaglom1999). The resulting anisotropy is shown in figure 8(a), and we observe that the spectrum is steeper in the parallel direction than in the perpendicular. For $\theta \rightarrow 0$, where $\theta$ is the angle between the local guide field and a chosen direction, the spectral index steepens to $-3$, which is the maximum that can be recovered using second-order structure functions as applied here (Farge & Schneider Reference Farge and Schneider2006). This result implies that there is no parallel cascade. For $\theta \gtrsim 20$, we find that the $E_{B_\parallel }$ spectrum is isotropic, which indicates the expected isotropy of a fast wave $E_{B_\parallel }$ spectrum (Cho & Lazarian Reference Cho and Lazarian2002,  Reference Cho and Lazarian2003; Chandran Reference Chandran2005). In figure 8(b), we present the energy in modes parallel and perpendicular to the local guide field. This allows us to estimate the energy distribution in Alfvén and fast waves. The latter have parallel perturbations of the magnetic field while Alfvén waves have only perpendicular fluctuations $\delta {\boldsymbol {B}} \perp {\boldsymbol {B}}_\textrm {guide}$. Most of the magnetic energy is associated with Alfvén waves, and the fraction of energy in fast waves remains low at 1–5 %. For comparison, we overplot the energy distribution in the parallel and perpendicular components to the global guide field $B_0$: they are shown with dashed lines in figure 8(b). The fact that the energies match closely with respect to local and global fields highlights the nature of weak turbulence: as $\tau _{nl} \gg \tau _A$, eddies predominantly interact with the $global$ guide field. The difference in the energy distribution of the parallel modes can be explained by projection of the energy in perpendicular modes on the parallel direction.

Figure 8. Angular properties for the turbulence at $t=120\ [2{\rm \pi} /\omega _0]$: (a) the spectral index of the power-law spectrum as a function of angle between the local guide field and a chosen direction ($\theta$); (b) the total energy contained in modes perpendicular and parallel to the local (solid lines) and global guide field (dashed lines) modes.

3.6. Weak turbulence and current sheets in relativistic magnetohydrodynamics

In this section, we compare the formation of current sheets in weak MHD turbulence to the force-free results presented in the previous sections. In figure 9, we compare force-free MHD and ideal MHD spectra obtained with BHAC for $\sigma =1$, $10$ and $100$ at $512$ cells per initial wavelength, which shows that in all four cases, a $E_{B_\perp } (k_{\perp })\propto k_{\perp }^{-2}$ magnetic energy spectrum forms and then a transition occurs between 100, shown in panel (a), and 130, shown in panel (c), wave-crossings induced by current sheet formation and break-up. After 120 wave-crossings, shown in panel (b), the total electromagnetic energy starts to decay (see figure 5b). Owing to energy transfer between electromagnetic and hydrodynamic (kinetic and thermal) components, the onset time of the steep decay and therewith the development of the spectrum differs slightly for varying $\sigma$, where $\sigma =100$ corresponds most closely to the force-free MHD result.

Figure 9. Evolution of spectra during the continuous wave interaction for different magnetizations. The $\hat {E}_{B_\perp }$ spectra at $100\pm 2$, $120\pm 2$ and $130\pm 2$ wave-crossings at resolution $N_{\perp }=N_{\parallel } = 512$ grid-cells for runs at $\sigma =100$, $10$ and $1$ and force-free MHD in BHAC. The spectra are averaged over five wave-crossings. For the force-free MHD run, we set $\eta =10^{-6}$. The $\hat {E}_{B_\perp } \propto k_{\perp }^{-2}$ spectra develop at 100 wave-crossings.

We show the structure of the forming current sheet shortly before it breaks-up at $t = 118\ [2{\rm \pi} /\omega _0]$ in the $\sigma =100$ MHD run in BHAC in figure 10(a), similar to the force-free current sheet in figure 1. The sheet breaks up at $t = 123\ [2{\rm \pi} /\omega _0]$. A pressure gradient contributes to the force balance to sustain current sheets in the relativistic MHD simulations. For smaller $\sigma$, the current sheets have a smaller aspect ratio and they have an effectively lower Lundquist number owing to the lower Alfvén speed and length. The thinning rate (shown in panel (c)) of the current sheet is very similar to the force-free result, which converges between $512^3$ and $1024^3$ grid points (see figure 3).

Figure 10. Current sheets form between elongated eddies in the MHD run for $\sigma =100$ with BHAC. (a,b) Out-of-plane current density $\hat {J}_z = [\nabla \times \boldsymbol {B}]_z / (k_{\perp }B_0)$ at 118 and 123 wave-crossings, respectively, before and after the break-up of the current sheet. (c) Current sheet's aspect ratio $\delta /w$ measured as the full-width at half maximum of the current density normalized by its approximately constant width $w\approx 2{\rm \pi}$ versus wave-crossings for both the $\sigma =100$ and force-free runs in BHAC at $512^3$ grid points.

4.1. Packet initialization

To allow for a direct comparison with the results of the previous sections, we only modify the wave initialization described in § 3.1 in two key aspects. First, the cubic box is extended along the direction parallel to the guide field, with $L_\parallel =10{\rm \pi}$. This choice allows us to separate the waves in packets while the choice of $k_\parallel$ re-scales the wave in the elongated domain. Second, the two Alfvén waves are localized along the (parallel) $z$-direction by window functions with a Gaussian profile (cf. Verniero et al. Reference Verniero, Howes and Klein2018):

(4.1)\begin{align} w(z)=\exp\left[-\left(\frac{z-z_0}{\varDelta_z}\right)^p\right]+\exp\left[-\left(\frac{z-z_0+L_\parallel}{ \varDelta_z}\right)^p\right]+\exp\left[-\left(\frac{z-z_0-L_\parallel}{\varDelta_z}\right)^p\right]. \end{align}

We choose a window width of $\varDelta _z = L_\parallel / 4{\rm \pi}$ and the power $p=2$. The two waves are localized at the window centres $z_{0-}=L_\parallel /4$ and $z_{0+}=3\times L_\parallel /4$. Initial phase shifts $\delta _-$ and $\delta _+$ align the waves symmetrically in their respective windows. For the tests presented in this section, we fix the nonlinearity parameter to $\chi = k_{\perp } \delta B_{\perp }/ (k_{\parallel }B_0)=0.25$. Hence, though the wave amplitudes are comparable to the continuous set-up (§ 3), the nonlinear time in the localized set-up is drastically reduced to $\chi ^{-2}=16$ wave interactions. This choice is beneficial for this numerical exploration. In contrast to the overlapping waves in § 3, the elongated domain requires each wave packet to propagate significant distances between the interaction events. Decreasing the nonlinear time thus allows us to reach the relevant stages of turbulence during a reasonable computational time.

We present insight into the collisional dynamics of localized Alfvén waves by combining two specific wave types: ${\mathrm {AW}}_1$, where $k_\parallel = 2{\rm \pi} /L_\parallel$ with an initial wave amplitude of $\delta B_{\perp } / B_0=0.05$ and $k_\perp /k_\parallel =5$; and ${\mathrm {AW}}_3$, where $k_\parallel = 3\times 2{\rm \pi} /L_\parallel$ with an initial wave amplitude of $\delta B_{\perp } / B_0=0.15$ and $k_\perp /k_\parallel =5/3$. Specifically, we examine the interactions $\mathrm {AW}_1+\mathrm {AW}_1$ and $\mathrm {AW}_3+\mathrm {AW}_3$ to understand the evolution of higher order modes (figure 11). We then focus on comparing the interactions $\mathrm {AW}_1+\mathrm {AW}_3$ and $\mathrm {AW}_3+\mathrm {AW}_3$ by closely following the effect of shearing modes and the ensuing dissipation of electromagnetic energy (figure 12). The wave vectors differ to those of the continuous case (§ 3) in their parallel component: $k_\parallel$ is both re-scaled and filtered by the window function, which results in the initial wave packets having a broad spectrum of $k_\parallel$ modes rather than $k_\parallel = \pm 1$. It is, thus, simplest to limit the representation of modes to the components of $k_\perp$, namely, $(k_x L_\perp,k_y L_\perp )=(1,0)$ for $\boldsymbol {k}^{+}_1$ and $(0,1)$ for $\boldsymbol {k}^{-}_1$. The mediator mode $\boldsymbol {k}^{0}_2$ is then $(1,1)$, with tertiary modes $(2,1)$ for $\boldsymbol {k}^{+}_3$ and $(1,2)$ for $\boldsymbol {k}^{-}_3$. The reality condition of the Fourier transform also implies a mirrored energy in the modes $(-1,0)$ for $\boldsymbol {k}^{+}_1$ and $(0,-1)$ for $\boldsymbol {k}^{-}_1$, as well as in their respective superpositions. The interaction of Alfvén wave packets proceeds through the coupling of these modes in (superficial) analogy to the continuous case presented in § 3. The details of this mechanism, however, reveal some notable subtleties.

Figure 11. Evolution of individual modes during the interaction: (a) comparison of the ${\mathrm {AW}}_3+{\mathrm {AW}}_3$ and ${\mathrm {AW}}_1+{\mathrm {AW}}_1$ collisions for the intermediate resolution of $256^2\times 320$; (b) decay of the primary modes $(1,0)$ and $(0,1)$ with sustained growth of the $(1,1)$ mediator mode as well as the secondary modes $(2,1)$ and $(1,2)$. The respective ratio between the total energy stored in the Alfvén and in the fast channel of the $(1,1)$ mode is displayed in panel (a). Exemplary decomposition of the $(1,1)$ mode by its polarization for the ${\mathrm {AW}}_3+{\mathrm {AW}}_3$ collisions: (c) total magnetic field along the centre of the domain; (d) projection into the Alfvén mode; (e) projection into the fast mode. An animation of the $(1,1)$ mode structure and dynamics during the first interaction cycles is provided in Media Supplement C (2021) and Media Supplement H (2021).

Figure 12. Decay of (normalized) total wave energy $\hat {E}_\textrm {EM} = E_\textrm {EM} / E_{\textrm {EM},t=0}$, with $E_\textrm {EM}=(B-B_0)^2 + E^2$ driven by the reconnection zones and current sheets induced by the shearing action of the secondary mode. (a,b,c) Visualization of the out-of-plane current density $\hat {J}_z = [\boldsymbol {\nabla } \times \boldsymbol {B}]_z / (k_{\perp }B_0)$ and field composition of the wave packets at different times in high resolution as indicated in the energy plot. The insets provide an impression of the $z$-averaged non-ideal field marker $\boldsymbol {E}\boldsymbol {\cdot }\boldsymbol {B}/\boldsymbol {B}^2$. Regions of strong (numerical) diffusion of non-ideal electric fields (i.e. violations of the force-free conditions) appear in dark shades. Such locations include Y-points, such as the anchor points of the shearing motion, and current sheets. (d) Evolution of the wave energy (normalized) for interactions of different wave types and resolutions. Animations of the interaction dynamics are provided in Media Supplement D (2021) and Media Supplement E (2021).

4.2. Interaction dynamics

Two localized Alfvén waves create linear fast waves as well as linear Alfvén waves during their interaction, and we exploit the mode structure in Fourier space to conduct the analysis of secondary modes. Such a transformation allows us to select specific coordinates of $k_\perp$ and $k_\parallel$ to identify and distinguish Alfvén and fast modes. As we show in figure 11(b), the interacting primary modes $(1,0)$ and $(0,1)$ transfer their energy to the $(1,1)$ mediator mode and, subsequently, to the tertiary waves $(1,2)$ and $(2,1)$. The energy transfer is strongest during each collision, as one can infer from the step-like growth in figure 11 (b). In analogy to Verniero & Howes (Reference Verniero and Howes2018), we examine two key characteristics of the secondary $(1,1)$ mode: field polarization and wave dispersion. Force-free linear Alfvén waves with dispersion $\omega (\boldsymbol {k})=\pm k_\parallel$ are composed of an electric field along $\boldsymbol {k}_\perp$ and a magnetic field along $\boldsymbol {\hat {z}}\times \boldsymbol {k}_\perp$. For the $(1,1)$ mode, this reduces to the following polarization and group velocity

(4.2a–d)\begin{equation} [{\mathrm{AW}}]\quad B^x ={\pm} E^y \quad B^y ={\mp} E^x\quad \frac{\partial\omega}{\partial k_\parallel}={\pm} 1. \end{equation}

Force-free fast waves with dispersion $\omega (\boldsymbol {k})=\pm |\boldsymbol {k}|$ are composed of an electric field along $\boldsymbol {k}_\perp \times \boldsymbol {\hat {z}}$ and a magnetic field along $\boldsymbol {k}\times (\boldsymbol {k}_\perp \times \boldsymbol {\hat {z}})$. For the $(1,1)$ mode, this reduces to the following polarization and group velocity

(4.3a–d)\begin{equation} {[}FW]\quad B^x ={\mp}\frac{k_\perp}{|\boldsymbol{k}|} E^y \quad B^y ={\pm} \frac{k_\perp}{|\boldsymbol{k}|}E^x\quad \frac{\partial\omega}{\partial k}={\pm} 1. \end{equation}

Both Alfvén waves and fast waves are relevant for the relativistic interaction of localized wave packets, and we present a summary of our detailed investigation in the following. For a fixed $\chi =k_\perp \delta B / (k_\parallel B_0)$, a smaller ratio $k_\perp / k_\parallel$ corresponds to a larger $\delta B / B_0$ such that the system is driven away from the reduced MHD limit $\delta B_\perp / B_0 \sim k_\parallel / k_\perp \ll 1$ (Paper I). Hence, the amount of energy stored in the FW increases for a smaller ratio $k_\perp / k_\parallel$ (as we show in figure 11a). During the initial phase of the interaction, fast waves in the $(1,1)$ mode show a progressive decay in energy after each primary wave collision. In contrast, $(1,1)$ Alfvén waves incrementally increase in energy during each primary interaction event and conserve their energy in between. Fast waves effectively propagate as spherical waves in a periodic domain (see also Media Supplement C 2021; Media Supplement H 2021). They, thus, continuously interact with Alfvén waves and other fast waves; energy is transferred out of the $(1,1)$ mode at a higher rate than for their Alfvénic counterpart. Figure 11(c,d,e) probes the decomposition of the wave fields in the $(1,1)$ mode along the aforementioned characteristic directions. For the Alfvén mode (analysed in figure 11d), it is then straightforward to identify the field polarizations derived in (4.2a–d). However, a quantitative analysis of the group velocity as well as the field polarizations of the fast wave (4.3a–d) is more difficult. This intricacy is rooted in the finite spectrum of $k_\parallel$ owing to the localization of the initial pulses, i.e. $k_\parallel$ is not a constant. A careful decomposition of wave modes for each individual $k_\parallel$ is required to confirm the remaining dispersion properties locally in Fourier space. We conducted such a brute-force validation. It shows that the $(1,1)$ mode is, indeed, the superposition of a linear Alfvén wave and a linear fast wave.

We note that in contrast to the continuous interaction of AWs, the localized $(1,1)$ mediator mode is a superposition of waves with a finite spectrum in $k_\parallel$. More specifically, it is not a nonlinear fluctuation with $k_\parallel =0$. This distinction of spectral components is the key to understanding the localized ${\mathrm {AW}}+{\mathrm {AW}}$ interaction. Counter-propagating Alfvén waves interact resonantly and satisfy the matching conditions

(4.4a,b)\begin{equation} \boldsymbol{k}_1^+{+}\boldsymbol{k}_1^-{=}\boldsymbol{k}_2,\quad\omega_1^+{+}\omega_1^-{=}\omega_2. \end{equation}

For these conditions to result in an outgoing AW, as identified in (4.2a–d), one of the incoming modes is required to have $k_\parallel =0$. Such a combination is possible when one wave packet of a finite $k_\parallel$ spectrum interacts with the $k_\parallel =0$ component of another wave packet (Ng & Bhattacharjee Reference Ng and Bhattacharjee1996). The resulting AW would then have a distribution of $k_\parallel$ that matches the primary mode; we find that this is indeed the case. An interaction between primary modes across the full spectral range of $k_\parallel$ also results in an outgoing FW. Owing to the addition of spectral ranges, the FW should have an accumulation of power at $k_\parallel =0$. We observe this feature consistently in our analysis.

The interaction of continuous counter-propagating AWs with the symmetry of constant and oppositely signed $k_\parallel$ has a unique solution to the resonance condition of (4.4a,b). Namely, it yields the nonlinear fluctuation $(1,1,0)$. The localization of waves in packets extends the spectrum of $k_\parallel$ along a finite range and the symmetry in the resonance is broken. In contrast to the continuous interaction (§ 3), the secondary mode has a linear Alfvén wave character and propagates according to the dispersion relation given above. During each interaction of the wave packets, energy is transferred to the secondary (and higher order) modes. The secondary mode manifests through shearing of the magnetic field. While, initially, wave packets form elongated tubes along their respective $k_\perp$, the persistent shearing stretches these tubes, gradually compressing them (cf. figure 12). The interaction of waves with symmetric $k_\parallel$ induces the $(1,1,0)$ oscillating shear mode. Thus, depending on the combination of initial $k_\parallel$, as to say, of the wave packets ${\mathrm {AW}}_1$ and ${\mathrm {AW}}_3$, the shearing motion can be continuously progressing (for non-symmetric $k_\parallel$) or oscillating (for symmetric $k_\parallel$).

In the case of ${\mathrm {AW}}_1+{\mathrm {AW}}_1$, the oscillation of the shearing process prevents the compression of anti-parallel field lines in between eddies into current sheets. However, elongated eddies, e.g. at $x\approx 3{\rm \pi} /2, y\approx {\rm \pi}/2$ at $t=31\ [2{\rm \pi} / \omega _0]$ in figure 12(a), develop and form pancake-like current sheets, e.g. at $x\approx 4{\rm \pi} /3, y\approx {\rm \pi}/3$ at $t=31\ [2{\rm \pi} / \omega _0]$ in figure 12(b). Prominent non-ideal electric fields form at these locations. The signatures of current sheets at later times (figure 12(b)) can be associated with a mild steepening of the energy dissipation. The shearing and stretching of larger eddies into reconnection layers are a characteristic of both continuous and localized Alfvén wave interactions.

In analogy to the case of continuous interaction (§ 3), elongated eddies can form long current sheets that will eventually break up (cf. § 3.2). As in the previous analysis of continuous interactions, we associate the onset of dissipation of cascaded energy with the emergence and break-up of current sheets. We identify sheets of strong current density, prominent non-ideal electric fields, and field reversals in the ${\mathrm {AW}}_1+{\mathrm {AW}}_3$ case (figure 12(c)) owing to field line compression between two eddies at $x\approx {\rm \pi}, y\in [0,{\rm \pi} ]$ in figure 12(c). They emerge and decay at the same time as we observe a change in slope in the energy evolution.

We note that at the relatively high value of $\chi$ for wave packet simulations, while we have seen a persistence of the $E_{B_\perp } (k_{\perp })\propto k_{\perp }^{-2}$ spectrum despite the presence of intermittency (in the form of current sheets), we begin to see some features of a transition to strong turbulence (Meyrand et al. Reference Meyrand, Galtier and Kiyani2016) which call for further studies.

We conclude this section by studying 3-D profiles of the wave interactions of the ${\mathrm {AW}}_3+{\mathrm {AW}}_3$ as well as the ${\mathrm {AW}}_1+{\mathrm {AW}}_3$ collision in figures 13 and 14. In these figures, we colour a large number of field lines seeded in planes at the centre of the respective wave packets (or their superposition) by the current density $|\boldsymbol {\nabla } \times \boldsymbol {B}|_z/(k_{\parallel }B_0)$ in panel (a,c) and by the local non-ideal electric field marker $\boldsymbol {E}\boldsymbol {\cdot }\boldsymbol {B}/\boldsymbol {B}^2$ in panel (b,d). First, we note that the wave magnetic field is well localized even after several tens of interactions. Second, possible 3-D reconnection sites emerge in the form of bulged eddies (figure 13) or current sheets with varying depth (figure 14). These structures exist individually in each wave packet and interact progressively, enhanced by the shearing action of the secondary mode. The $z$-extensions of non-ideal structures is a fraction of the initial window-size $\varDelta _z$ in the case of the bulged eddies. Current sheet structures, however, have a significant depth of the order of (and even slightly above) the initial window-size $\varDelta _z$. Figure 15 dissects the topology of the current sheet emerging during the ${\mathrm {AW}}_1+{\mathrm {AW}}_3$ collision, which shows some anti-parallel field lines and some reconnecting field lines. We mark the localized layer of magnetic null points with red dots (i.e. points with reversals of all field components, where the initial background field $B_0$ along $\hat {\boldsymbol {z}}$ is removed).

Figure 13. Three-dimensional impressions of the ${\mathrm {AW}}_3+{\mathrm {AW}}_3$ collision during a phase of packet separation ($t=31\ [2{\rm \pi} /\omega _0]$, (a,c)) and the subsequent interaction ($t=31.5 [2{\rm \pi} /\omega _0]$, (b,d)). We colour field lines along $\hat {\boldsymbol {B}}$ by the current density $|\boldsymbol {\nabla } \times \boldsymbol {B}|_z/(k_{\parallel }B_0)$ (a,c) and by the local non-ideal electric field marker $\boldsymbol {E}\boldsymbol {\cdot }\boldsymbol {B}/\boldsymbol {B}^2$ (b,d). The vortex-like anchor points of the shear (cf. figure 12) are extended along the $z$-direction and coincide with regions of non-ideal diffusion. The inset in panel (c) shows a top view of the maximum current density to estimate the depth of the current structures. The different colours in panel (b) distinguish between field lines as seeded in the two central planes of the displayed wave packets. A 3-D animation of the dynamics in the interaction region can be found in Media Supplement F (2021).

Figure 14. Three-dimensional impressions of the ${\mathrm {AW}}_1+{\mathrm {AW}}_3$ collision in analogy to figure 13 ($t=31\ [2{\rm \pi} /\omega _0]$, (a,c)); $t=31.5\ [2{\rm \pi} /\omega _0]$, (b,d)). A current layer separates two equally polarized eddies. It coincides with field reversals, strong currents and significant dissipation by emerging non-ideal electric fields. A detailed outline of the reconnection region (green box) with realistic proportions is presented in figure 15. A 3-D animation of the dynamics in the interaction region can be found in Media Supplement G (2021).

Figure 15. Detailed 3-D outline of the reconnection region during the ${\mathrm {AW}}_1+{\mathrm {AW}}_3$ collision (as presented in figure 14). (a) Magnetic null points (red spheres, where we removed the initial background field $B_0$ along $\hat {\boldsymbol {z}}$) between eddies of the same polarity. Yellow field lines are seeded at the locations of the strongest currents, black field lines have a slight offset to these seeds in the perpendicular plane. (b) Volume rendering of the strongest currents (as we showed for the continuous case in figure 4).

Collisions of localized Alfvén waves are efficient mediators of energy conversion to smaller scales (perpendicular to the guide field) and powerful drivers of nonlinear energy cascades. The mixing of wave packets with different ratios of $k_\perp /k_\parallel$ significantly changes the interaction dynamics by suppressing the oscillatory (1,1,0) shear mode. Mixing effectively triggers episodes of very fast dissipation of electromagnetic energy mediated by the break-up of elongated current sheets. As in the continuous case (cf. figure 4), the current sheet ripples as it is stretched between elongated eddies and eventually breaks up into smaller (turbulent) structures. Mixing waves of different $k_\perp /k_\parallel$ is an appealing case for a future exploration of current sheet formation and the transition into the strong turbulence regime for the interaction of overlapping Alfvén waves in high-resolution runs.