Generation of Turbulence in Colliding Reconnection Jets

Magnetic reconnection is a fundamental phenomenon in astrophysical plasmas. It consists of the recombination of magnetic field topology due to the violation of the frozen-in law for magnetic field of ideal magnetohydrodinamics (MHD). Due to this recombination, plasma is ejected in the form of jets and magnetic energy is converted into kinetic energy. The presence of magnetic reconnection has been ascertained in the solar corona (Shibata 1998; Sui et al. 2004) and in the Earth's magnetosphere (Retinò et al. 2007; Burch et al. 2016), and it is thought to be responsible for magnetic energy release and plasma heating (Drake et al. 2006). Reconnection has also been identified in the solar wind (Phan et al. 2006; Gosling 2007). A close link exists between the presence of magnetic reconnection and the phenomenology of plasma turbulence (Matthaeus & Velli 2011). In plasmas, similarly to fluid dynamics, magnetic and velocity fluctuation energy cascade from large to small scales, due to nonlinear interactions, following the turbulent phenomenology. This nonlinear cascade produces smaller and smaller scale magnetic shears that can eventually undergo the process of magnetic reconnection. Magnetic reconnection, indeed, can be viewed as an active ingredient of plasma turbulence (Matthaeus & Lamkin 1986; Servidio et al. 2009, 2011; Franci et al. 2017).

On the other hand, magnetic reconnection can act as a trigger of plasma turbulence in the reconnection outflows (Matthaeus & Lamkin 1986; Malara et al. 1991, 1992; Lapenta 2008; Huang & Bhattacharjee 2010; Beresnyak 2016). Numerical simulations have shown that reconnection outflows can comprise different types of plasma instabilities at both fluid (Guo et al. 2014) and kinetic scales (Vapirev et al. 2013; Huang et al. 2015). These instabilities, along with the presence of shears, make reconnection outflows turbulent rather than laminar. The latter is not just numerical evidence, but has been recently confirmed by in situ spacecraft observations (Eastwood et al. 2009; Osman et al. 2015). Contrary to the "classical" view of reconnection energetics, where the conversion of magnetic into kinetic energy happens only in the main reconnection sites (Shay et al. 2007), reconnection outflows are now also believed to be regions where the plasma is heated and particles are accelerated (Daughton et al. 2011; Lapenta et al. 2014, 2015). It has been shown, both in numerical simulations (Leonardis et al. 2013; Pucci et al. 2017) and from the most recent in situ observation (Fu et al. 2017), that, due to turbulence, energy exchange between fields and particles in the outflows is intermittent. This means that the greatest part of energy transfer happens in low volume-filling, very intense current sheets. Recently, numerical simulations (Olshevsky et al. 2016) and observations (Fu et al. 2017) have shown that energy exchange is more efficient in the presence of magnetic null points of spiral topological type (O-point configurations) rather than the radial nulls (X-points). This suggests than energy dissipation is more efficient in a region surrounded by two, or multiple, X-lines. In the classical two-dimensional MHD picture of magnetic reconnection, an O-point is always present between the two X-points (for topological reasons), i.e., in the middle of two counter-propagating reconnection jets, if the X-lines are active.

The physics of the collision of reconnection jets has been studied numerically with both fluid and kinetic models (Karimabadi et al. 1999; Oka et al. 2008; Nakamura et al. 2010; Markidis et al. 2012), sometimes in the frame of the tearing or plasmoid instability (Bhattacharjee et al. 2009; Landi et al. 2015). Recently, the first observations of colliding reconnection jets have been carried out in the presence of both weak (Alexandrova et al. 2016) and strong (Øieroset et al. 2016) guide fields. Observational results show that collisions of jets cause wave activity and are responsible for secondary reconnection events in the outflows, but a clear theoretical picture of the collision of reconnection jets is still lacking.

Even though such observational results have been compared with two-dimensional simulations of magnetic reconnection, it is likely that a three-dimensional description would be closer to the physical reality. In fact, the nature of magnetic reconnection in three dimensions is intrinsically different from its two-dimensional counterpart (Priest et al. 2003). Numerical simulations have shown that 3D instabilities with out-of-plane wave numbers develop during reconnection both in the case of weak (Ricci et al. 2004b; Vapirev et al. 2013) and strong guide fields (Daughton et al. 2011; Fermo et al. 2012). The presence of these instabilities highly enriches the dynamics and the properties of the generated turbulence. For this reason, it is worth also studying the collision of reconnection jets in three dimensions.

In this work, we present a three-dimensional (3D) numerical simulation to study the physics of the collision of reconnection jets, in the presence of a strong guide field, at kinetic scales. We show how the collision of jets gives rise to plasma turbulence in a confined region of the outflows. We provide a numerical method to detect the turbulent region and perform local statistical analysis to characterize the turbulence that is produced in such regions. We believe that these results can be used for the interpretation of present and future observations of collisions of jets at ion and subion scales. Moreover, we provide a theoretical description of turbulence generated by magnetic reconnection in a three-dimensional kinetic model.

We consider a plasma made of ions (protons) and electrons. At the initial time, the plasma is at the Harris equilibrium

$$\begin{eqnarray*}&&{\boldsymbol{B}}={B}_{0x}\tanh (y/\delta ){{\boldsymbol{e}}}_{x}+{B}_{0z}{{\boldsymbol{e}}}_{z},\,n={n}_{b}+\displaystyle \frac{{n}_{0}}{{\cosh }^{2}(y/\delta )},\end{eqnarray*}$$

where B_0z is the guide field, δ is the thickness of the Harris sheet, and n_0b is the background density. The coordinates are chosen such that the equilibrium sheared magnetic field is directed along the x axis and varies along y, causing the initial current sheets to be directed parallel to the z axis, which is also the guide field direction. For both species, the initial distribution function is a Maxwellian with homogeneous temperature. The system obeys the Vlasov–Maxwell equations, which are numerically solved using the semi-implicit particle-in-cell (PIC) code iPIC3D (Brackbill & Forslund 1982; Markidis et al. 2010; Lapenta 2012).

The semi-implicit method gives the advantage of allowing bigger grid size and time step compared with explicit scheme, resulting in a much lower computational cost. The price to be paid for that is that not all relevant spatial and temporal scales within the plasma, such as the electron inertial length and the inverse of the electron plasma frequency, are resolved. However, even though electron scale dynamics is poorly resolved, kinetic effects from both ions and electrons are retained. The semi-implicit method has been compared with explicit methods in previous numerical studies of magnetic reconnection, showing a positive agreement (Ricci et al. 2004a). Moreover, iPic3D has recently been used to study plasma turbulence in the context of magnetic reconnection (Pucci et al. 2017) and of Taylor-Green vortex decay (Olshevsky et al. 2018). These works have shown iPic3D capability to well reproduce the properties of subion scale turbulence as predicted by theoretical models and observed in space plasmas.

The numerical domain is a Cartesian box with size [40, 30, 10] d_i, where d_i is the ion inertial length. We use [512, 384, 128] cells, each one initially populated with 125 particles of each species. A reduced ion/electron mass ratio m_i/m_e = 256 is used, which fixes the spatial resolution to Δx = 1.25d_e, where d_e is the electron inertial length. The initial ion/electron temperature ratio is T_i/T_e = 5, and thermal velocities u_th,i = 0.0063c and u_th,e = 0.045c, where c is the speed of light. The thickness of the initial current sheet is δ = 0.5d_i and the background density is n_b0/n₀ = 1/10. The value of the guide field is set to B_0z/B_0x = 2. The Alfvén to light speed ratio is c_a/c = 0.0217 and the electron plasma to cyclotron frequency ratio is ω_pe/Ω_ce = 2.88. We use an open boundary condition in the y direction and we impose periodicity along x and z. The equilibrium is perturbed by a Gaussian fluctuation of the z component of the vector potential, located at (x = 0, y = 7.5), which initializes the magnetic reconnection process (Lapenta et al. 2010). The system evolution is followed in time up to 29τ_i, where τ_i is the ion gyro-period, using a time step of dt ≈ τ_e/10, where τ_e is the electron gyro-period. Due to open boundary conditions and numerical dissipation associated with the semi-implicit method, the total energy in the system is not conserved. However, these two effects are small and in the final state of the simulation the total energy differs from the initial energy by 3%.

Because of the strong guide field, we expect that current sheets will be preferentially aligned along the z direction. In Figure 1, we represent a time evolution of the reconnecting current layer. The quantity plotted in panels (a–e) is $\langle {J}_{z}{\rangle }_{z}$, the out-of-plane total current density as a function of x and y and averaged along z. The initial Harris sheet is a narrow layer, directed along z, as reported in panel (a), and it has been perturbed at the boundary x = 0. Because of periodicity, two reconnection sites form at opposite sides of the domain in the x direction. Panel (f) shows the evolution of the mean square ion, electron, and total current $\langle {J}^{2}\rangle$. At the beginning, the current is mainly carried by ions, due to the Harris equilibrium, but as the evolution proceeds electron current becomes dominant. The total current reaches a peak value at $84.06{{\rm{\Omega }}}_{{ci}}^{-1}$ and rapidly decreases to a smaller value. Then, the value of $\langle {J}^{2}\rangle$ remains almost constant, slightly decreasing up to the end of the simulation.

Zoom In Zoom Out Reset image size

The X-lines produce two reconnection jets that approach one against the other in the middle of the simulation domain (panel (b)). Notice that the two outflows carry oppositely directed current sheets even before the actual collision. This is due to the pile up of regions with oppositely directed B_y which are generated and squeezed one against the other by the pressure of the two counter-propagating jets. When the two outflows collide, the maximum of $\langle {J}^{2}\rangle$ is reached (panel (c)). At that instant, the two outflows can still be distinguished even though they look like they are caught in each other. The two counter-propagating outflows also carry with them oppositely directed magnetic fields, which makes the system unstable to secondary reconnection. The large-scale current sheets present in panel (c), whose dimensions are of the order of ∼5–10d_i, are fragmented in smaller ones by multiple reconnection events. The effect is visible in panel (d), where a multitude of current sheets at different scales appears in the collision region. This suggests the beginning of a turbulent cascade, that through nonlinear interactions forms smaller current sheets. At later times, the turbulence that develops in the center of the domain continuously generates and disrupts current sheets and the associated magnetic islands. It is worth noting that the area where the current sheets are located decreases in time (panel (e)). This is reasonable considering that the system is not driven externally, but is just relaxing after the collision of jets.

We have just presented a 2D evolution of the system, relying on the fact that the presence of a strong magnetic field favors the formation of perpendicular gradients (Oughton et al. 1994, 2015). To verify this assumption, in Figure 2, we present a comparison of $\langle {J}_{z}{\rangle }_{z}$ with J_z at a fixed z, namely z = 5.0d_i. The plot is taken in the turbulent phase at $t=92.19{{\rm{\Omega }}}_{{ci}}^{-1}$ (panel (d) in Figure 1), right after the collision of jets. $\langle {J}_{z}{\rangle }_{z}$ and $\langle {J}_{z}\rangle (z=5.0)$ in the xy plane are plotted in the first and second columns, respectively. The latter looks like the blurry version of the former. This can be interpreted as follows. The current sheets are elongated in the z direction and, for this reason, they are not ruled out by the average. At the same time, a small level of fluctuations with wavevectors in the z direction are present, which make the picture of the 2D cut blurry. Those fluctuations are ruled out when averaging along z. In the third column, a 1D cut of the two quantities is plotted as a function of x. We can see that several current sheets are encountered moving in the x direction. The strongest ones are found in the central region. There is no big difference between the two variables, confirming that the fluctuations in the z direction are small compared to the average current intensity. In the second and third rows we show two zooms of the two quantities taken in the central part of the simulation. We observe that the two quantities also look similar at smaller scales. Moreover, we can notice that the size of the smallest current sheets are of a fraction (1/5 or less) of the ion inertial length d_i, corresponding to ∼3d_e, where d_e is the electron inertial length.

Zoom In Zoom Out Reset image size

This preliminary analysis has shown that a turbulent behavior is observed after the collision of the two jets. The comparison between $\langle {J}_{z}{\rangle }_{z}$ and ${J}_{z}(z=5.0)$ suggests that strong fluctuations are produced perpendicularly to the guide field, and small amplitude fluctuations are present along the guide field. In order to make the last statement more quantitative, we have analyzed the second-order structure function of the magnetic field. This is defined as ${S}_{B}^{2}({\boldsymbol{r}})=\tfrac{1}{V}\int | {\boldsymbol{B}}({\boldsymbol{x}}+{\boldsymbol{r}})-{\boldsymbol{B}}({\boldsymbol{x}}){| }^{2}{d}^{3}x$, where B is the total magnetic field, r is a displacement in the physical space, and V is the total volume of the simulation box. Similar to the power spectral density, this quantity represents the average of the magnetic fluctuation energy between two points of the domain separated by a lag r. Since we are interested in evaluating the energy of the fluctuations varying in the direction parallel and perpendicular to z, i.e., the guide field direction, we have computed the so-called reduced structure function defined as ${S}_{B}^{2}({{\ell }}_{i})=\tfrac{1}{V}\int | {\boldsymbol{B}}({\boldsymbol{x}}+{{\ell }}_{i}{{\boldsymbol{e}}}_{i})-{\boldsymbol{B}}({\boldsymbol{x}}){| }^{2}{d}^{3}x$, with i = x, y, z. These quantities, plotted in Figure 3, represent the domain average of the magnetic energy fluctuations with wavevector in the i direction at a scale ℓ. For statistical reasons, we consider maximum lags ${{\ell }}_{i}^{\max }={L}_{j}/2$, where L_j is the box size in the j direction. The minimum lag is set to ${{\ell }}_{j}^{\min }={\rm{\Delta }}{x}_{j}=1.25{d}_{e}$. The results confirm what was suggested by the visual analysis of the current. The blue curve, representing ${S}_{B}^{2}({{\ell }}_{z})$, is always below the red and green curves, representing ${S}_{B}^{2}({{\ell }}_{x})$ and ${S}_{B}^{2}({{\ell }}_{y})$, respectively. The magnetic fluctuation energy is larger for wavevectors perpendicular to the guide field and smaller for parallel ones at all scales. Moreover, the fluctuations in the direction of the initial shear, i.e., along y, are slightly stronger than the fluctuations in the x direction. In order to check if this difference is due to the presence of the background Harris field, we subtract its contribution to the total magnetic field, defining the new following variable: ${\boldsymbol{b}}={\boldsymbol{B}}-\langle {{\boldsymbol{B}}}_{x}(y){\rangle }_{x,z}{{\boldsymbol{e}}}_{x}$ $-{B}_{0z}{{\boldsymbol{e}}}_{z}$. The variable b can be seen as the magnetic fluctuation field. ${S}_{b}^{2}({{\ell }}_{y})$ is plotted in Figure 3 as a green dashed line. The fluctuation energy in the y direction is stronger than that in the x direction even when the contribution of the Harris field is not taken into account. Both are at least one order of magnitude stronger than the fluctuations in the z direction. We can also notice that at large scales the structure function tends to saturate only in the z direction. The saturation is reached at the scale of the correlation length. Since large-scale shears are present in both the x and y directions, the corresponding structure function does not saturate for lags as large as half the box size in each direction. Consistently, when the contribution of the Harris field is subtracted from the total field, the structure function in the y direction saturates at around 5d_i. This size is comparable to the size of the turbulent region depicted in Figure 1 (panel (d)) and Figure 2.

Zoom In Zoom Out Reset image size

The analysis of ${S}_{B}^{2}$ contributes to the result that turbulence develops in the box after the collision of jets. This turbulence is mainly two-dimensional, with fluctuation energy concentrated in wavevectors perpendicular to the reconnection guide field. A second smaller anisotropy is found for the in-plane fluctuations, where more energy is present in the direction of the initial magnetic shear. The analyses performed so far are global, in the sense that the full computational domain was used as a single sample. However, the turbulent region seems to be confined in a smaller central region surrounding the location where the outflows have collided. In the next section, we provide a method to isolate that turbulent region and study its properties separately from the rest of the domain.

In this section, we provide a method for isolating the turbulent region and we describe the turbulence analysis performed on the region itself comparing it with the rest of the simulation box.

4.1. The Method

In the previous section, we showed that turbulence develops in the simulated reconnection event, due to the collision of jets, and seems to be localized in the surrounding of the collision site. We have also shown that such turbulence is quasi two-dimensional and develops in the plane perpendicular to the guide field. We build our method bearing on this last feature and considering the system as being two-dimensional, thus neglecting the dependence on the z component. In 2D MHD, the in-plane components of the magnetic field is a function only of the out-of-plane component of the vector potential A. From ${\boldsymbol{B}}={\rm{\nabla }}\times {\boldsymbol{A}}$, one finds that ${B}_{x}=\tfrac{\partial {A}_{z}}{\partial y}$, and ${B}_{y}=-\tfrac{\partial {A}_{z}}{\partial x}$, where xy is the 2D plane and A_z is the z component of A. From the last two relations it is easy to show that $({\boldsymbol{B}}\cdot {\rm{\nabla }}){A}_{z}=0$. The last expression means that A_z is constantly moving along an in-plane field line or, equivalently, that in-plane magnetic field lines are isocontours of A_z. It is also easy to show that

$$\begin{eqnarray}&&-{({\rm{\nabla }}\times {\boldsymbol{B}})}_{z}={{\rm{\nabla }}}^{2}{A}_{z}.\end{eqnarray} \tag{ 1 }$$

This equation will be used later.

Operationally, the two-dimensional ingredient is obtained by averaging the magnetic field B in space, over the z direction. The z-average of the 2D magnetic field is similarly related to the z-average of the full 3D vector potential, that is, the 2D potential ${A}_{z}=\hat{z}\cdot \langle {\boldsymbol{A}}{\rangle }_{z}$. This averaging operation has the effect of eliminating the fluctuations in the z direction, which are, however, much smaller than those in x and y, as shown previously. The map of the 2D potential A_z is plotted in panel (b) of Figure 4, along with its isocontours, which are the in-plane magnetic field lines. In order to compute A_z, we solved Equation (1) using a Fourier method in an expanded (periodic) domain, then plotted the result in the original domain. The magnetic island produced by reconnection in the center of the box is clearly visible. Moreover, we notice that A_z has a minimum in the middle of the island and a saddle point at the location of the main X-line. This is consistent with what is found in incompressible MHD simulations of 2D magnetic reconnection (Matthaeus 1982). It is also consistent with Equation (1), as we explain in the following. The presence of a minimum implies that ${{\rm{\nabla }}}^{2}{A}_{z}\gt 0$, and consequently ${({\rm{\nabla }}\times {\boldsymbol{B}})}_{z}\lt 0$. From Maxwell equations, we know that ${\rm{\nabla }}\times {\boldsymbol{B}}=\tfrac{1}{c}\tfrac{\partial {\boldsymbol{E}}}{\partial t}+\tfrac{4\pi }{c}{\boldsymbol{J}}$. At large scales and slow frequencies, the displacement current term can be neglected and the curl of B becomes proportional to J. Considering that the current in the initial current sheet is negative for construction, the minimum of A_z in the middle of the island is justified. In panel (a) of Figure 4, we plot a 2D map of $\langle {J}_{z}{\rangle }_{z}$ at time $t=92.19{{\rm{\Omega }}}_{{ci}}^{-1}$ and two isocontours of A_z, i.e., A_z = −0.098 and A_z = −0.138. The former represents the value of A_z at the X-point and at the separatrices, the latter represents a manually selected value of A_z associated with a magnetic field line enclosing most of the more turbulent central region and very little of the more quiescent (low current density) region surrounding it. These two isocontours identify three different regions, plotted in panel (c): Region 1 (R1), out of the separatrices; Region 2 (R2), between the separatrices and the turbulent region; Region 3 (R3), the central turbulent region.

Download figure:

Video Standard image High-resolution image

Since we want to reproduce such a distinction for each instant of time we have implemented the following numerical procedure. First, at a given time t, we compute the in-plane magnetic field average along z and the corresponding 2D vector potential. Second, we find ${A}_{z}^{X}$, namely the value of A_z at the X-point, and ${A}_{z}^{\min }$, namely the minimum of A_z located in the center of the domain. We call ${A}_{z}^{* }$ the suitable value of A_z associated with the magnetic field line that encases the turbulent region. Since this value has to lay between ${A}_{z}^{\min }$ and ${A}_{z}^{X}$, we can write it in the following form:

$$\begin{eqnarray}&&{A}_{z}^{* }={A}_{z}^{\min }+a({A}_{z}^{X}-{A}_{z}^{\min }),\ \mathrm{with}\ 0\lt a\lt 1.\end{eqnarray} \tag{ 2 }$$

Third, we find the suitable value of a as that associated to a contour that encases all the strongest current sheets in the turbulent region. We found a = 1/6 to be suitable for time $t=92.19{{\rm{\Omega }}}_{{ci}}^{-1}$. Then, we repeat the procedure for each time step, keeping the selected value of a unchanged. The procedure is applied from the time when the mean square current is maximum, i.e., $t=84.06{{\rm{\Omega }}}_{{ci}}^{-1}$, up to the end of the simulation. Before that time it has no meaning separating the domain into regions, since turbulence has not developed yet. Choosing a fixed value for a does not mean in principle that ${A}_{z}^{* }$ does not change in time, since both ${A}_{z}^{\min }$ and ${A}_{z}^{X}$ can change. It rather means that we are choosing a value of the vector potential whose relative distance from the minimum and the value at the separatrices is constant. Interestingly, this choice gives a correct result for all time steps (see the animation associated with Figure 4). This can be explained considering that, from Maxwell's equations, $\tfrac{\partial {A}_{z}}{\partial t}=-{{cE}}_{z}$, meaning that the value of A_z changes due to the reconnection electric field. After the initial phase when jets are ejected, the value of this field becomes smaller causing the vector potential to be nearly constant in time on the separatrices. This method allows us to identify different regions in the simulation domain in 2D. Since small fluctuations are present in the z direction, we can extend the 2D region along z and perform the statistical analysis on the 3D extension of each region.

4.2. The Analysis

We have identified the turbulent region looking at the out-of-plane component of the total current, thus associating turbulence with the current activity. In Figure 5, we show the mean square current in the three regions R1, R2, and R3, as a function of time (panel (a)). The full box average is plotted as a reference. The mean square current in R3 is much higher than in R1 and R2. Its trend in time looks very similar to the trend of the full box average, shown in panel (f) of Figure 1. Panel (b) shows the filling factor F of each region as a function of time, computed as the ratio between the volume of each region and the volume of the whole box. R3 is the smallest region and occupies around 6% of the total volume. The size of the region does not change significantly in time. This analysis shows that the most intense current sheets are within region R3, which is a small portion of the total domain. However, region R3 is large enough to perform a statistical analysis of fluctuations contained therein.

Zoom In Zoom Out Reset image size

Since the chosen regions are not rectangular, and their boundaries are not periodic, a Fourier-based spectral analysis is not a suitable approach. For this reason, as done for the full box analysis, we compute the second-order structure functions of the magnetic field in the three regions for lags in the x, y, or z directions. In Figure 6, we report the result at time $t=92.19{{\rm{\Omega }}}_{{ci}}^{-1}$. We notice first that the energy for parallel increment lags (panel (c)) is smaller than for perpendicular lags (panels (a) and (b)), in each region (panel (c)). Moreover, different regions have similar parallel fluctuations energy at all scales. Panels (a) and (b) show the second-order structure function of the in-plane magnetic fluctuations. Due to the difference in shape and size of the regions, the curves do not span the same range of lags. The maximum possible lag is in fact defined as half of the sample size. The red curve, associated to region R3, terminates at smaller lags because it is obtained in the smallest region R3. Moreover, the maximum ℓ_x for R3 is bigger than the maximum ℓ_y, since the shape of the region is elongated in the x direction. At large scales, a saturation in the value of the second-order structure function is observed just for R3. The lag at which the saturation is observed gives an estimate of the correlation length of around 5d_i along x and 3d_i along y, which is comparable to the sizes of the region. Such saturation is not observed for R1 and R2, due to the presence of the large-scale magnetic shear whose dimension is comparable to the sizes of those regions. At scales smaller than the correlation length, the red curve (R3) is always higher than green (R2) and blue (R1), ${S}_{B,{\rm{R}}3}^{2}\gt {S}_{B,{\rm{R}}2}^{2}\gt {S}_{B,{\rm{R}}1}^{2}$.

Zoom In Zoom Out Reset image size

The average energy of magnetic fluctuations varying in the plane perpendicular to the guide field is larger in R3, implying that the nonlinear cascade has acted more efficiently in that region. Therefore, this proves that R3 is the actual region where turbulence develops.

In order to confirm that the turbulent cascade proceeds up to the end of the simulation, we have computed ${S}_{B}^{2}$ in each region for different times. At each time we take the value of ${S}_{B}^{2}$ for ${{\ell }}_{j}={d}_{i}$ for j = x, y, z. The result is plotted in Figure 7. In the perpendicular directions the average of the magnetic fluctuation energy at scale d_i is larger in R3, meaning that more energy is cascading to subion scales in R3 compared to the rest of the domain. The energy in the y direction is larger, on average, due to the presence of the Harris field. In the z direction, the values in R3 eventually become larger as the simulation proceeds. At earlier times, the large energy in R2 is due to the formation of an electron instability that develops on the separatrices (Daughton et al. 2011; Fermo et al. 2012). This instability remains confined in R2 and saturates for times larger than $t\sim 100{{\rm{\Omega }}}_{{ci}}^{-1}$.

Zoom In Zoom Out Reset image size

To complete our analysis, we study the Eulerian frequency spectra of the electric and magnetic field in R3. The fields have been sampled during the simulation at a cadence of 3Δt by means of probes ("virtual satellites") placed at predefined locations inside the simulation box. For the analysis we present here, we have used satellites from a region of dimensions $1.25\times 1.25\,\times 10.0\,{d}_{i}^{3}$ located in the middle of R3. The selected region contains 3 × 3 × 24 = 216 equally spaced satellites. We analyze the data in the time frame $92.18{{\rm{\Omega }}}_{{ci}}^{-1}\lt t\lt 187.09{{\rm{\Omega }}}_{{ci}}^{-1}$, extending from the initial development of the turbulence to the end of the simulation. Similar to what happens for in situ spacecraft observations, magnetic field data are not periodic. Moreover, the situation we are considering here is not stationary, since turbulence in R3 is decaying in time. In order to avoid spurious effects coming from signal detrending or windowing, we use the technique proposed by Bieber et al. (1993) to compute the frequency power spectra. This method has been used, e.g., for the analysis of the recent in situ solar wind observations (Bale et al. 2005). The method follows these steps. First, the series of the increments are computed from the actual field. If B(t) is the magnetic field as a function of time, its increments are defined as ${\rm{\Delta }}{\boldsymbol{B}}(t)={\boldsymbol{B}}(t+{\rm{\Delta }}t)-{\boldsymbol{B}}(t)$. Second, the autocorrelation function of the increments is computed as ${R}_{{ii}}^{({\rm{\Delta }})}(\tau )\,=\tfrac{1}{V}\int {\rm{\Delta }}{{\boldsymbol{B}}}_{i}(t){\rm{\Delta }}{{\boldsymbol{B}}}_{i}(t+\tau ){dt}$. Third, from the autocorrelation function, the power spectrum from the increments ${P}_{{ii}}^{({\rm{\Delta }})}(f)$ is computed via Fourier transform. Finally, the spectrum of the actual field P_ii(f) is recovered by filtering the spectrum of the increments:

$$\begin{eqnarray*}&&{P}_{{ii}}(f)=\displaystyle \frac{{P}_{{ii}}^{({\rm{\Delta }})}(f)}{4{\sin }^{2}(\pi /{\rm{\Delta }}t)}.\end{eqnarray*}$$

This procedure suppresses the contribution of the low frequencies unresolved in the chosen sample. Computing power spectral density from the autocorrelation function can yield to nonphysical results such as negative power spectrum values (Blackman & Tukey 1958); therefore, we limit our analysis to lags smaller than τ_max = T/10, where T is the duration of the total sample. With this upper limit, we have enough statistics to estimate the spectrum for all frequencies selected.

Figure 8, reports the power spectra of magnetic (panel (a)) and electric (panel (b)) field as a function of frequency. The spectra are computed for each probe separately and then averaged. In Figure 8, both spectra of the three components and total field are plotted. Characteristic frequencies are marked by vertical lines. Reference lines are shown corresponding to spectral indices of −3 and −1. The computed spectra span frequency ranges that extend from a fraction of the ion cyclotron frequency to beyond the electron cyclotron frequency. The electron plasma frequency is poorly resolved. Note that the range of plotted frequencies does not reach the Nyquist frequency because the fields are saved each 3Δt. Magnetic and electric frequency spectra from fluid to electron scales in plasma turbulence have recently been studied by means of in situ spacecraft observations in the magnetosheath by Matteini et al. (2017). It has been shown that the magnetic and electric spectra have approximately similar behavior in the fluid-MHD like regime, while they show different trends in the kinetic range. In particular, the magnetic field spectrum steepens at ion scales, passing from a spectral index between −5/3 and −3/2 in the MHD fluid range to a roughly −3 spectral index at ion and subion scales. The electric field has an opposite trend, passing from a spectral index −3/2 in the fluid range to roughly −1 in the kinetic range. The familiar theoretical explanation for this phenomenon (e.g., Matthaeus et al. 2008) depends on dominance in the kinetic range of nonideal terms in the generalized Ohm's law. When the turbulence k-vectors are mainly perpendicular to the guide field, this heuristic explanation predicts $\delta {E}^{2}\propto {k}^{2}\delta {B}^{2}$ for an electromagnetic fluctuation at a given wavevector k in the subion range.

Zoom In Zoom Out Reset image size

The above reasoning explains the relationship of observed frequency and wavenumber spectra in cases in which the two are proportional and related by the Taylor hypothesis, appropriate for high speed super-Alfvénic flows. This explains the agreement of frequency spectra coming from single spacecraft observations (Bale et al. 2005; Eastwood et al. 2009; Matteini et al. 2017), with wavenumber spectra obtained from numerical simulations (Franci et al. 2017; Pucci et al. 2017).

The results we report here seem to be, at first sight, in contrast to the above scenario, in that the subion scale magnetic and electric spectra in our analysis do not follow a specific power law and strongly differ from a ∝f⁻³ and a ∝f⁻¹ trend, respectively. The main reason for this apparent discrepancy relies on the fact that Taylor hypothesis is not valid in the system considered here, as we presently demonstrate.

The Taylor hypothesis consists of assuming that time variations at a single observational point are only due to the rapid sweeping of spatial structure past that point. This requires that the large-scale bulk flow of the plasma greatly exceeds all characteristic dynamical speeds that might distort those structures. This includes both wave propagation and turbulent velocity (Perri et al. 2017). When Taylor hypothesis is not verified, there is not a simple connection between the frequency and wavenumber spectra, although there remains the possibility of a statistical similarity between the Eulerian frequency spectra and the wavenumber spectra (Chen & Kraichnan 1989; Matthaeus & Bieber 1999) when the dominant effect is random large-scale sweeping.

The Taylor hypothesis is not valid in the system we are considering here. In order to assess this hypothesis, the bulk velocity of the plasma, computed in the region where the virtual satellites are located, must be compared with characteristic wave speeds and with turbulent velocity fluctuations. In Figure 9, we report the comparison between the bulk velocity V, the Alfvén velocity V_A, and the ion sound speed V_s. The bulk velocity V is defined as the plasma center of mass velocity, the Aflvén velocity as ${{\boldsymbol{V}}}_{{\rm{A}}}={\boldsymbol{B}}/\sqrt{4\pi \rho }$, where ρ is the ion density, and the ion sound speed as ${V}_{S}=\sqrt{\tfrac{\gamma {\sum }_{j}{P}_{{jj}}}{3\rho }}$, where j = (x, y, z), P is the ion pressure tensor, and γ is the adiabatic index, which we have chosen equal to 5/3 for simplicity. These characteristic macroscopic speeds are compared in Figure 9, which shows the x, y, and z components of V and V_A, and the ion sound speed (considered isotropic). The bulk speed is in every case smaller or comparable with the other two characteristic speeds. In particular, the ion sound speed is the larger speed in the directions perpendicular to the guide field, while the Alfvén speed is, as expected, the largest speed in the parallel direction.

Zoom In Zoom Out Reset image size

Next the bulk speed is compared with measures of the turbulent fluctuation velocities in Figure 10. In particular, we compute the longitudinal increments $u({\boldsymbol{r}},{\boldsymbol{\ell }})=\hat{{\boldsymbol{\ell }}}\cdot [{\boldsymbol{V}}({\boldsymbol{r}}+{\boldsymbol{\ell }})-{\boldsymbol{V}}({\boldsymbol{r}})]$, where r is a position in space, and ℓ is the direction along which the velocity increments are computed. The root mean square velocity increments are computed in terms of components along the three Cartesian directions as ${u}_{\mathrm{rms}}({{\ell }}_{i})=\sqrt{\langle {\boldsymbol{u}}({\boldsymbol{r}},{{\boldsymbol{\ell }}}_{i})\cdot {\boldsymbol{u}}({\boldsymbol{r}},{{\boldsymbol{\ell }}}_{i})\rangle }$. Here $\langle \rangle$ indicates the average in the region where the probes are located, and ℓ_i with i = x, y, z indicates the direction along which the velocity increments are computed. We considered a lag size equal to d_i for each direction. In this way, we examine, in particular, whether the Taylor hypothesis is valid for the kinetic features seen in the frequency spectra in Figure 8. Figure 10 shows a comparison of the bulk velocity with these measures of fluctuations. Each panel contains the ith component of the bulk velocity and the corresponding ${u}_{\mathrm{rms}}({{\ell }}_{i})$. First, we notice that ${u}_{\mathrm{rms}}({{\ell }}_{z})$ is smaller than ${u}_{\mathrm{rms}}({{\ell }}_{x})$, and ${u}_{\mathrm{rms}}({{\ell }}_{y})$, in accordance with anisotropic, quasi-2D turbulence developed in the plane perpendicular to the guide field. Second, we notice that V_z the bulk speed along the guide field direction is larger than all u_rms. This large speed could in principle be responsible for a similarity between the k_z spectrum and the frequency spectrum. However, as shown in the previous section, turbulence is mainly 2D due to the presence of the strong guide field, and very little energy is contained in the k_z modes. Moreover, the largest characteristic speed in the z direction is the Alfvén speed (see panel (c) in Figure 9), which allows for the presence of fast parallel propagating waves of Alfvènic nature that can invalidate the Taylor hypothesis.

Zoom In Zoom Out Reset image size

Based on this analysis of characteristic speeds, we can conclude that the Taylor hypothesis is not verified in the location where the Eulerian spectra are computed. Therefore, a similarity between frequency and wavenumber spectra is not expected. The magnetic and electric frequency spectra presented in Figure 8 are therefore to be considered as independent of the spatial (wavenumber) structure. Then the observed features are essentially temporal, and the observed features indicate the presence of wave activity at different temporal scales, which we now briefly describe.

First, we notice a flattening of the magnetic spectrum between the ion and electron cyclotron frequency, around f = 10⁻¹, relative to the ∝f⁻³ trend at lower frequencies. This feature is consistent with the presence of whistler waves and is also found in in situ observations (Matteini et al. 2017). To confirm this hypothesis, we applied a pass band filter at a frequency f_ce/5 to the magnetic field (not shown here) finding a phase shift between B_x and B_y consistent with parallel propagating whistler.

Second, a significant peak at the electron cyclotron frequency is observed in both electric and magnetic field spectra. The peak is found for all three components of the fields probably indicating the presence of both electromagnetic and electrostatic modes. To verify that this peak is not due to a numerical artifact, we ran a new simulation (not shown) with a halved time step and saving the fields data at each computational cycle. No relevant change in the results was observed, either in the peak intensity or in its position in frequency. Therefore, we conclude that the peak is associated with physical electron cyclotron modes.

It is worth noting how the aforementioned presence of wave modes at different timescales is not in contrast with a picture of fully developed turbulence. Previous MHD studies have proven that Eulerian frequency spectra can present intensity peaks at fixed frequencies associated with wave modes even when no signature of such modes is present in the wavenumber spectra (Dmitruk et al. 2004; Dmitruk & Matthaeus 2009). A result similar to the one presented here was found recently in two-dimensional hybrid simulations of kinetic plasma turbulence, where several peaks at frequency larger than the ion cyclotron frequencies have been found in the Eulerian spectra of the magnetic field (Parashar et al. 2010). Here we show such a behavior in a turbulent environment that is self-consistently generated by magnetic reconnection. A detailed description of the modes that contribute to shape the frequency spectra along with the effects due to the lack of validity of Taylor hypothesis goes beyond the scope of this work and will be considered for future studies.

We presented a full kinetic 3D simulation of reconnection in the presence of a strong guide field with reduced mass ratio. The simulation considered periodic boundary conditions in the direction of the outflow allowing counter-propagating reconnection jets to collide. We showed that the collision of reconnection jets coming from two neighboring X-lines drives turbulence in the magnetic island within them.

Turbulence produced is neither homogeneous nor stationary. In the central region, where reconnection jets collide, turbulence is mainly two-dimensional and manifests in the continuous disruption and formation of new current sheets mainly directed along the guide field, whose thicknesses range between the electron and ion inertial lengths. This phenomenology is due to the presence of the strong out-of-plane guide field, given that the contribution of the initial Harris field is small in that region. On the separatrices, an electron instability typically observed in 3D simulations of magnetic reconnection with strong guide field develops with wavenumber directed out of plane (Daughton et al. 2011; Fermo et al. 2012). This instability remains confined on the separatrices and saturates during the evolution before that the turbulent energy in the central region is completely dissipated. This results in having and overall 2D turbulence for later times. How this picture changes with different values of the background field is an interesting issue and it is left for future works. We believe that the effect of the instability would be decreased and that turbulence in the central region would become more isotropic if the value of the guide field is reduced.

For times longer than $\sim 100{{\rm{\Omega }}}_{{ci}}^{-1}$, dynamical activity in the out-of-plane electric current density remains well confined in a small region of the simulation box that occupies a small percentage (6%) of the total volume. Because the system remains fairly close to two-dimensionality we were able to exploit a representation that is strictly applicable in the 2D case and, possibly, will not be applicable in a strongly three-dimensional case. Using an identification method based on the magnetic vector potential, we numerically separated the turbulent region from the rest of the computational domain.

Analysis of spatial second-order structure functions show that the selected central region is where turbulence actually develops allowing an efficient cascade of magnetic energy to subion scales. This cascade remains active for the duration of our simulation. Since there is no forcing except for the initial one due to the collision of jets, the current activity slightly reduces in time, remaining confined in the same region.

Power spectrum analysis of electric and magnetic fluctuations in the turbulent region reveal interesting properties that, at first sight, stand in contrast to a number of previous observations and numerical simulations. However, the assumption of Taylor hypothesis made in previous works do not apply to this case, which explains the apparent discrepancy. Under these conditions, the Eulerian frequency spectra computed in our analysis are quite independent of the spectra computed in the Taylor hypothesis scenario, which are interpreted as wavenumber spectra. There is no simple relationship between spatial structure and temporal structure in the present case. Interestingly, the analysis of Eulerian frequency spectra in the turbulent region reveals the presence of wave activity at different frequencies: parallel propagating whistler between the ion and electron cyclotron frequency, and electromagnetic and electrostatic modes at the electron cyclotron frequencies. The presence of whistlers in the region compressed by counter streaming reconnection jets was recently found in situ observation in the case of the small guide field (Alexandrova et al. 2016). Electrostatic and electromagnetic wave activity at the electron cyclotron frequency has also been detected in previous observations of magnetic reconnection in the vicinity of the reconnection sites (Tang et al. 2013; Viberg et al. 2013; Khotyaintsev et al. 2016). Here, we observed electron cyclotron wave activity in the turbulent region formed by the collision of reconnection jets.

In a recent paper, Øieroset et al. (2016) have shown the first observation of the collision of reconnection jets in the presence of a strong guide field. The dynamical picture presented in that work imagined an event of secondary reconnection activated by the counter-propagating jets both carrying magnetic field lines with opposite polarities. Our simulation shows that the physical situation can be less laminar than predicted and that turbulence develops after the collision of jets producing several current sheets, secondary reconnection events, and wave modes up to the electron cyclotron frequency. Our results stand in analogy to the recent result by Fu et al. (2017) who find turbulent reconnection in a region surrounding a secondary O-point. Like that case, our interpretation of the present result is that energy release can occur near such secondary structure in a turbulent reconnection environment. However, a more detailed analysis of energy exchange proxies, such as ${\boldsymbol{J}}\cdot {\boldsymbol{E}}^{\prime}$, where ${\boldsymbol{E}}^{\prime}$ is the electric field in the electron reference frame, is needed in order to confirm such interpretation. This kind of analysis is left for future work. The present results are likely to be relevant for the interpretation of future observations of collisions of reconnection jets, which may have distinctive signatures in space and astrophysical plasmas.

The present work was supported by the H2020 Project DEEP-ER and DEEP-EST, by the Onderzoekfonds KU Leuven (Research Fund KU Leuven, GOA scheme and Space Weaves RUN project), and by Fonds Wetenschappelijk Onderzoek (FWO) short stay grant for study abroad K218017N. This research used resources of the National Energy Research Scientific Computing Center, which is supported by the Office of Science of the US Department of Energy under Contract No. DE-AC02-05CH11231. Additional computing has been provided by NASA NAS and NCCS High Performance Computing, by the Flemish Supercomputing Center (VSC) and by PRACE Tier-0 allocations.