Kinetic Plasma Turbulence Generated in a 3D Current Sheet With Magnetic Islands

In this article we aim to investigate the kinetic turbulence in a reconnecting current sheet (RCS) with X- and O-nullpoints and to explore its link to the features of accelerated particles. We carry out simulations of magnetic reconnection in a thin current sheet with 3D magnetic field topology affected by tearing instability until the formation of two large magnetic islands using particle-in-cell (PIC) approach. The model utilizes a strong guiding field that leads to the separation of the particles of opposite charges, the generation of a strong polarization electric field across the RCS, and suppression of kink instability in the “out-of-plane” direction. The accelerated particles of the same charge entering an RCS from the opposite edges are shown accelerated to different energies forming the “bump-in-tail” velocity distributions that, in turn, can generate plasma turbulence in different locations. The turbulence-generated waves produced by either electron or proton beams can be identified from the energy spectra of electromagnetic field fluctuations in the phase and frequency domains. From the phase space analysis we gather that the kinetic turbulence may be generated by accelerated particle beams, which are later found to evolve into a phase-space hole indicating the beam breakage. This happens at some distance from the particle entrance into an RCS, e.g. about 7d i (ion inertial depth) for the electron beam and 12d i for the proton beam. In a wavenumber space the spectral index of the power spectrum of the turbulent magnetic field near the ion inertial length is found to be −2.7 that is consistent with other estimations. The collective turbulence power spectra are consistent with the high-frequency fluctuations of perpendicular electric field, or upper hybrid waves, to occur in a vicinity of X-nullpoints, where the Langmuir (LW) can be generated by accelerated electrons with high growth rates, while further from X-nullponts or on the edges of magnetic islands, where electrons become ejected and start moving across the magnetic field lines, Bernstein waves can be generated. The frequency spectra of high- and low-frequency waves are explored in the kinetic turbulence in the parallel and perpendicular directions to the local magnetic field, showing noticeable lower hybrid turbulence occurring between the electron’s gyro- and plasma frequencies seen also in the wavelet spectra. Fluctuation of the perpendicular electric field component of turbulence can be consistent with the oblique whistler waves generated on the ambient density fluctuations by intense electron beams. This study brings attention to a key role of particle acceleration in generation kinetic turbulence inside current sheets.

The recent space observations of current sheets in the magnetosphere and heliosphere (Fujimoto & Sydora, 2008;Zhou et al., 2009;Huang et al., 2016;Pucci et al., 2017;Eastwood et al., 2018;Phan et al., 2020) and 2D/3D full kinetic and Hall-MHD simulations (Daughton et al., 2004;Matthaeus & Velli, 2011;Roytershteyn et al., 2012;Boldyrev et al., 2013;Franci et al., 2017;Loureiro & Boldyrev, 2017;Papini et al., 2019;Pezzi et al., 2021) had already pointed to a link between reconnection and turbulence. Current sheets contain a sufficient amount of free energy that is released by instabilities in collisionless plasmas at the smallest, kinetic scales often revealing in both hybrid-kinetic and Hall-MHD turbulence simulations the onset of energy transfer at the smallest scales as soon as reconnection is triggered (see for details Matthaeus and Velli, 2011;Papini et al., 2019;Pezzi et al., 2021, and references therein). The small-scale turbulence in a vicinity of those CSs was usually associated with spectral breaks in the magnetic fluctuation spectra near the ion cyclotron frequency Ω ci . At larger scales (low frequencies), there is the characteristic inertial range of the turbulent cascade, while below ion scales the turbulent spectra show a clear power law with spectral indices close to −2.7-2.8 (Boldyrev et al., 2013;Franci et al., 2017;Loureiro & Boldyrev, 2017;Pucci et al., 2017;Muñoz & Büchner, 2018). Moreover, the power laws and spectral breaks near CSs are very similar to those measured in homogeneous turbulent solar wind plasmas (Chen et al., 2008;Zhou et al., 2009;Huang et al., 2016;Eastwood et al., 2018;Phan et al., 2020).
Also one of the longest-known instabilities connected with reconnection is the lower hybrid drift instability (LHDI) long suspected to play a role in reconnection (process as observed in space (Cattell et al., 2005;Chen et al., 2008;Divin et al., 2015;Artemyev et al., 2016) and in the laboratory (Carter et al., 2001). Such LHDI occurs near the lower hybrid frequency ω lh ωpi 1+ω 2 pe /Ω 2 ce , where ω pe (ω pi ) is the electron (ion) plasma frequency, Ω ce is the electron cyclotron frequency (Muñoz & Büchner, 2018). However, all these observations do not yet have the certain answers regarding which processes of a reconnection contribute to the measured turbulent spectra. The kinetic turbulence in reconnecting current sheets has been extensively investigated (see, for example, Drake et al., 2003;Fujimoto & Machida, 2006;Fujimoto, 2014;Muñoz & Büchner, 2018;Lapenta et al., 2020, and references therein). Cattell et al. (2005) observed the electron holes in the separatrix regions similar to the prediction of 3D PIC simulations (Drake et al., 2003) that are considered to be the nonlinear evolution of the bump-in-tail instability, or Buneman instability (Omura et al., 1996). Lapenta et al. (2020) identified the two regimes of turbulent fluctuations in current sheets: one in the outflow leading to a turbulent regime where the fluctuations involve both fields and particles and the other in the inflow and separatrix region, which involves only the electromagnetic fields, without significantly affecting the particles. The two regimes differ much in practical consequences. The outflow regime is capable of inducing a strong and turbulent energy exchange as well as strong anomalous momentum exchange dominated primarily by the electrostatic term in Ohm's law. The inflow regime, in contrast, does not lead to substantial fluctuations in the field-particle energy exchange nor significant anomalous viscosity or resistivity limiting turbulence to the electromagnetic fields only. However, the authors presented a more intuitive interpretation of the detected turbulence obtained from PIC simulations without linking it to the regimes of particle acceleration during magnetic reconnection in the presence of magnetic islands.
To understand these kinetic instabilities generated in reconnecting current sheets one needs to explore acceleration of particles dragged into the reconnection region and to investigate the turbulence generated by them. For this reason, we need to refresh our views about the properties of accelerated particles gained during their passage through a reconnecting current sheet with a single and multiple X-nullpoints and to explore which of them, if any, can lead to the formation of turbulence and in what locations. Since the plasma turbulence introduced by beam instabilities is, in general, inherently a 3D problem in PIC simulations (Goldreich & Sridhar, 1995;Siversky & Zharkova, 2009;Muñoz & Büchner, 2018), it requires the simulation domain for acceleration of particles in current sheets to be a 3-dimensional one. The theoretical and numerical studies of magnetic reconnection are typically performed using a simplified system of 2D antiparallel reconnecting magnetic fields with an additional out-of-plane guiding magnetic field (B g ) in the third dimension. Such RCSs with a finite B g are observed in Earth magnetopause (Silin & Büchner, 2006) and at the impulsive phases of flares and CME eruptions (Fletcher et al., 2011). Owing to large magnetic field gradients and curvatures surrounding the reconnection sites, combined with strong gradients of the plasma temperature and density, the electromagnetic fields vary dramatically inside reconnecting current sheets (RCSs) (Shay et al., 2016;. Furthermore, thin elongated RCSs formed in the diffusion region between the reversed magnetic field lines are often broken down by tearing instability into multiple islands, or O-type nullpoints separated by X-nullpoints (Furth et al., 1963;Loureiro et al., 2007;Bhattacharjee et al., 2009). The presence of magnetic islands in reconnecting current sheets was demonstrated by magnetohydrodynamic (Biskamp, 1986;Loureiro et al., 2005;Drake et al., 2006;Lapenta, 2008;Bárta et al., 2011) and kinetic simulations (Huang & Bhattacharjee, 2010;Karimabadi et al., 2011;Markidis et al., 2012). Such chains of magnetic islands have been identified in many solar flares Lin et al. (2005); Oka et al. (2010); Bárta et al. (2011);Takasao et al. (2012); Nishizuka et al. (2015) and CMEs (Song et al., 2012), in the in-situ observations in the heliosphere (Zharkova & Khabarova, 2012;Khabarova et al., 2015Khabarova et al., , 2021 and Earth magnetotail (Zong et al., 2004;Chen et al., 2008;Wang et al., 2016).
In the case of full 3D RCSs, the guiding field is accepted varying in time and space. In some configurations of 3D RCSs, the out-of-plane variations of the helical magnetic structures become pretty significant, due to the kink instability, obscuring current sheet structures and making it hard to define clear X-nullpoints Egedal et al., 2012). A strong guiding field B g can suppress the out-of-plane kink instability while leaving the concept of magnetic islands still applicable (Lapenta & Brackbill, 1997;Daughton, 1999;Cerutti et al., 2014;Sironi & Spitkovsky, 2014). Nevertheless, further studies have shown that both the cases do not significantly change the scenarios of energy conversion and particle acceleration in 3D RCSs, because the dominant mechanisms of particle energization remain the same as in the 2.5D scenario (Hesse et al., 2001;Zharkova et al., 2011;Guo et al., 2014;Dahlin et al., 2017).
Depending on magnetic field topologies, the presence of a guiding field in an RCS would cause partial or full charge separation between electrons and ions (Pritchett & Coroniti, 2004;Zharkova & Gordovskyy, 2004) because they gyrate in the opposite directions in a magnetic field. This, in turn, can lead to the preferential ejection of the oppositely charged particles into the opposite semiplanes of CSs, or opposite footpoints of reconnecting loops. It makes the hard X-ray sources spatially separated from the c − ray sources in the opposite footpoints of reconnecting magnetic loops Hurford et al., 2003Hurford et al., , 2006. This charge-separation phenomenon is also confirmed in the laboratory experiments (Zhong et al., 2016).
Furthermore, there is a polarization electric field in RCSs confirmed by 3D PIC simulations (Fujimoto, 2006;Zenitani & Hoshino, 2008;Cerutti et al., 2013;Fujimoto, 2014) but its nature was not clear and sometimes mixed with the parallel electric field of accelerated electrons. Then it was shown that the polarization electric field is induced across the reconnection current sheet midplane by the separation of particles of opposite charges (electrons and protons) during their acceleration in current sheets with a strong out-of-plane guiding field; and its magnitude is much larger (by two orders of magnitude) than a reconnecting electric field itself (Siversky & Zharkova, 2009;Zharkova & Agapitov, 2009). Furthermore, the spatial profiles of a polarization electric field were found dependent on magnetic field topologies because this electric field is induced by the separated electrons and protons (Siversky & Zharkova, 2009;Zharkova & Agapitov, 2009;Zharkova & Khabarova, 2012). The presence of polarization electric field is shown to explain the in-situ observations of ion velocity profiles during spacecraft crossings of the heliospheric current sheet, which are found to follow closely the profiles of polarization electric field (Zharkova & Khabarova, 2012. Therefore, the ambient plasma feedback to a presence of accelerated particles during their passage through reconnecting current sheets is very important for the particles of opposite charges. However, the particles of the same charge entering the 3D RCS from the opposite edges would also lead to different energy gains by the particles with the same charge (Siversky & Zharkova, 2009;Zharkova & Khabarova, 2012;Khabarova et al., 2020). The particles that enter the RCS from the side opposite to that, to which they are to be ejected, are classified as "transit" particles, while the particles entering the RCS from the same side where they are to be ejected to, are classified as "bounced" particles. The transit particles gain significantly more energy because they become accelerated on their way to the midplane where the main acceleration occurs, while bounced particles lose their energy while they approach the midplane, thus, gaining much less energy in the current sheet (Zharkova & Gordovskyy, 2005;Siversky & Zharkova, 2009;Zharkova & Agapitov, 2009;Zharkova & Khabarova, 2012).
The energy difference between the transit and bounced particles creates the particle beams with "bump-in-tail" velocity (energy) distributions, which could trigger different two beam instabilities (Buneman, 1958) and naturally generate plasma turbulence. Although, strong turbulence very often appears in the off-plane guiding field direction at the very early stages of 3D PIC simulations of magnetic reconnection Egedal et al., 2012) that obscures any other types of turbulence present in the simulations at later times. And, of course, the kinetic turbulence generated in current sheets can also contribute to particle acceleration by modifying the parameters of accelerated particles (Zharkova & Agapitov, 2009;Drake et al., 2010;Matthaeus & Velli, 2011;Fujimoto, 2014;Muñoz & Büchner, 2016;Huang et al., 2017;Trotta et al., 2020).
The goal of the current research is to explore kinetic turbulence generated by accelerated particles in reconnecting current sheets with multiple X-an O-nullpoints based on the specifics of particle acceleration on 3D magnetic field topologies. As one can note, the accelerated particles definitely gain non-Maxwellian (power-law) distributions during their acceleration in current sheets. Hence, we will attempt to explore the conditions in the phase and frequency domains for energetic particle beams to maintain the pressure anisotropy (Le et al., 2013) and their effects on instabilities generated due to asymmetric acceleration by a reconnection electric field. In addition, we wish to explore anisotropy of the electric and magnetic field fluctuations in turbulence along and perpendicular to the local mean magnetic field B m0 (Howes et al., 2008;Boldyrev et al., 2013) for different locations inside a reconnection region.
The simulation model and magnetic field topology are described in Section 2, the results of simulations of energetic particles and generate turbulence for a current sheet with single and multiple X-nullpoints are presented in section 3 and the general discussion and conclusions are drawn in section 4.

Magnetic Field Topology
In the current article, unlike our previous simulation (Siversky & Zharkova, 2009;, we do not separate the original and induced electromagnetic fields, and adopt the self-consistent 3D PIC simulation to investigate particle acceleration in magnetic islands generated by a magnetic reconnection. However, we will use the previous results  about particle acceleration in the similar reconnection scenarios to evaluate possible mechanisms of the recorded kinetic turbulence. We extend the 3D simulation region to a larger domain compared to the previous 2.5D studies (Siversky & Zharkova, 2009;Muñoz & Büchner, 2016).
The simulations start with a Harris-type current sheet in the x − z plane: where d cs is the half thickness of RCS. The B 0 is the initial guiding field, which is perpendicular to the reconnection plane. In the presented simulation b g B 0y /B 0z 1.0. The initial density variation across the CS is: where n 0 is the ambient density in a current sheet, n b is the density of an accelerated particle beam, and d cs is a current sheet thickness.

Particle Motion Equations
The motion of a charged particle in an electromagnetic field E and B is computed by the relativistic Lorentz equations: where V( p/mc) and r are the particle velocity and position vectors, q and m are the charge and the rest mass of the particle. p is the momentum vector and c is the corresponding Lorentz factor defined as c 1/ 1 − V 2 /c 2 √ . E and B are calculated from the initial electro-magnetic fields and the ones induced by accelerated particles as described in section below.

The Plasma Feedback
Similarly to the early article , in the initial PIC approach we split the electromagnetic field E and B into two components, the background E static and B static , and the local self-consistentẼ andB induced by the particle motions (Eq. (4)): B B static +B, and E E static +Ẽ. Then the fluctuation fields are calculated by the Maxwell solver: where j e and j p are the current densities of electrons and protons updated by the particle solver. The Maxwell's equations are solved by standard finite-difference time-domain method (FDTD) numerically. This approach can help us to identify the effect of the ambient particles that are drifted into a current sheet and accelerated. Then we rerun the 3D PIC simulations by relaxing all the electromagnetic fields and following the reconnection process until the certain time when maximal turbulence is formed.

Numeric Method
After clarifying the accelerated particle dynamics by splitting the electro-magnetic fields as above, we rerun the PIC simulations with VPIC code by relaxing electromagnetic fields of particles and allowing them to interact together with the initial electromagnetic field to reflect a reconnection process initiated by some perturbation. PIC simulations were carried out using the fully relativistic 3D VPIC code (Bowers et al., 2008). Our setup is somehow similar to the one employed in Muñoz & Büchner (2018) with some essential differences. The RCS thickness was d cs 0.5d i (versus 0.25 d i by Muñoz & Büchner, 2018), where d i is the ion inertial length. We chose a mass ratio m i /m e 100, a temperature ratio T i /T e 5, a background plasma density n b /n 0 0.2 versus n b / n 0 1.0 accepted by Muñoz & Büchner (2018), and a frequency ratio ω pe /Ω ce 1.5, where ω pe is the electron plasma frequency and Ω ce is the electron gyro-frequency. Plasma beta is estimated as β e β i 2μ 0 n 0 k B T i /B 2 0 ≈ 0.012 versus 0.016 in Muñoz & Büchner (2018).
Following the approach discussed by Siversky & Zharkova (2009), for the current sheet thickness equal to the ion inertial length, d i , we select the number of cells across the current sheet in a PIC simulation to be di λD cm i /(kT), which is 3 10 3 for the solar corona temperature or 3 · 10 4 for the magnetosphere. To reduce this number, Drake et al. (2006) used a reduced magnitude for the speed of light c 20V A 6·10 6 ms −1 , where V A is the Alfven velocity. Another way to reduce the number of cells was used in the PIC simulation carried out by Karlický (2008), who considered the high-temperature electron-positron plasma, for which the ratio d i /λ D was as low as 10.
The simulation box size is L x × L y × L z 12.8d i × 1.6d i × 51.2d i with grid number 512 × 64 × 2048 using 100 particles per cell. To avoid the problem with the small Debye length λ D , only a small fraction of the plasma particles (with a density of 10 12 m −3 10 6 cm −3 ) is included in the current PIC simulation. This makes the ratio λD di in the current simulations the order of 0.0192, e.g. the mesh step ratio d/λ D 1.3 that is close to that of 1.4 used by Daughton et al. (2011) for the same VPIC code. Hence, this mesh is safe and does not require any corrections on possible numerical stabilities of the explicit PIC code using the linear shape function (Birdsall & Langdon, 1991).
Along the direction x, the conducting boundary condition for the electromagnetic field and open boundary condition for particles are used. The periodic boundary conditions are applied along the z-and y-directions (in the current sheet midplane X 0) to the electromagnetic field and particles. We use a real speed of light without scaling it to Alfven speed, while using a reduced mass ratio between protons and electrons, like Siversky & Zharkova (2009) did. This approach is valid for the coronal magnetic fields only while the density would need to be modified if applied to current sheets in the magnetosphere or heliosphere as the applied setting can lead to larger than real Alfven velocities in the Earth magnetosphere.
To trigger a magnetic reconnection in the plane with magnetic islands, we introduce a small perturbation at the beginning of the simulation, which is written in terms of (δB 0 . . .) in Eq. 1, where δB 0 0.03B 0z . It comes from an out-of-plane vector potential, δB 0 ∇ × δA y , where δA y ∝ cos(2π z−0.5L z L z )cos(π x L x ) satisfying ∇·A 0. This spatial distribution helps us to set the fast reconnection to occur near the center of the simulation box in Figures 1A-D, similar to that reported earlier .
We will gather the kinetic turbulence in the whole simulation region at the particular moment when turbulence is stabilized (experiment 1). Also we will collect the kinetic turbulence data by a hypothetical spacecraft sampling the simulation domains at a few particular points with respect to the local mean magnetic field B m0 (experiment 2). Because the streaming instabilities are often observed in the separatrices (current sheet midplanes) and at the exhaust regions (Cattell et al., 2005;Lapenta et al., 2011;Markidis et al., 2012;Zhang et al., 2019;Lapenta et al., 2020), the positions of the virtual spacecraft are to be simultaneously located in the three points close to the separatrices at different distances away from the X-nullpoints inside the current sheet structure that forms a magnetic island.
Given the relativistic velocities of accelerated particles, which generate the turbulence within a very short timescale after the acceleration start, we can safely assume that any Doppler shifts in the frequencies of turbulence induced by accelerated particles caused by the motion of the ambient plasma particles inside a current sheet are negligible, because the motion of charged particles in an RCS strictly follows rigidly the magnetic field topology completely forgetting its initial velocity at the entry (Zharkova & Gordovskyy, 2004, 2005Dalla & Browning, 2005;Wood & Neukirch, 2005;Siversky & Zharkova, 2009;Xia & Zharkova, 2018;.

Single X-Nullpoints
To understand the physical nature of the turbulence generated inside RCSs with magnetic islands, let us use the models described in our previous papers (Xia & Zharkova, 2018;, which compared particle acceleration in a single X-nullpoint and in coalescent and squashed magnetic islands. The current sheet with a single X-nullpoint was described by the set of equations with the following magnetic field components: ; B y −B 0 ξ y , and a reconnection electric field E y 250 V/m with the current sheet plane to be x − z plane, where d is a current sheet thickness and a is its length (Xia & Zharkova, 2018).
In the PIC approach, there is also a feedback of the ambient plasma considered to the presence of accelerated particles by calculating the electric and magnetic fields induced by accelerated particles as described by Eqs. 5, 6 in Section 2.4. Similarly to Siversky & Zharkova (2009), in the PIC code the authors  introduced the initial (static) background electric and magnetic fields (Verboncoeur et al., 1995;Bowers et al., 2008) and then followed particle acceleration as well as their induced electric and magnetic fields in the current sheets with the single or multiple X-nullpoints (with magnetic islands).
This approach can help us to separate the original magnetic field configuration of the magnetic reconnection from that induced by the plasma feedback due to the presence of accelerated particles. This separation helps to discover potential triggers of plasma turbulence inside these complex magnetic configurations.

Polarization Electric Field
The trajectories of electron and protons calculated in the RCS near a single X-nullpoint for a strong guiding field B y reveal a significant difference between the acceleration paths of the particles with opposite charges. The particles with different charges are shown separated into the opposite sides from the RCS midplane and then ejected to the opposite semi-planes (Siversky & Zharkova, 2009;Xia & Zharkova, 2018;. For a given magnetic topology, energetic electrons can primarily be ejected to the x > 0 semi-plane, while protons to the x < 0 semiplane. One important outcome of this separation is the polarization electric fields induced by the separated particles with opposite charges across the current sheets. This polarization electric field δE x shown in Figure 2 is perpendicular to the RCS midplane, and it is much larger than the reconnecting electric field E y0 induced by the magnetic reconnection process. A polarization electric field was first reported in the 2D PIC simulations by Arzner & Scholer (2001); Fujimoto (2006) and was assigned to particle's inertia motion. However, the particles passing through 2D current sheets do not gain much energy (Litvinenko & Somov, 1993;Litvinenko, 1996) and, as a result, the polarization electric field induced by these accelerated particles owing to separation by inertia would have low magnitudes, in comparison with the reconnection electric field magnitude accelerating particles. Only later by considering acceleration of particles in 3D current sheets with a strong guiding field (Pritchett & Coroniti, 2004;Zharkova & Gordovskyy, 2004;Pritchett, 2005;Zharkova & Gordovskyy, 2005), this polarization electric field was shown to be enforced by significant energy gains by all particles and the separation of electrons from protons/ions across the current sheet midplane. This separation of very energetic electrons and protons generates a significant polarization electric field exceeding by up to two orders the reconnection electric field magnitude (Siversky & Zharkova, 2009;Zharkova & Agapitov, 2009;Zharkova & Khabarova, 2012).
In our further simulations, the plasma density is accepted to vary as 10 8 m −3 and 10 12 m −3 . The polarization electric field distributions are found sensitive to the ambient plasma density as shown in Figure 3B If the density is low, the particle separation is more distinguishable in the phase space as shown in Figure 3A. However, the polarization electric field induced in the more rarified ambient plasma is lower than in the dense plasma. This happens, we believe, because the gradient of magnetic field (the first term in Eq. (5)) remains the same while being much smaller than the currents of accelerated electrons and protons, which are increased for more dense plasma, thus making higher the resulting electric field E x induced by these accelerated particles in denser plasma.
Besides, there is a bump-in-tail at high-energy electrons in the spectrum of Figure 3D which is clearly seen for lower density plasma. When the polarization electric field, E x , becomes larger with a larger density (the charged particle density should also increase) as shown in Figure 3B, the preferential ejection becomes less clear, and the bump-in tail in the particle energy spectrum is smoothed out. Although, this does not change the maximum energy gains by particles as shown by the spectra in Figure 3D, which still remain of the same order of magnitude for all the simulations with different plasma densities.

Plasma Turbulence Generated by Two Beams
Because of the bump-on tail distribution of the energy spectra of accelerated particles shown of Figures 3A,C, there is turbulence formed by Buneman instability (Buneman, 1958), or the electron two stream instability, which, in addition to the background electro-magnetic fields, leads to fluctuations of electric δE x , δE y , δE z and magnetic field vectors |δB x /B x,0 |, |δB y /B y,0 |, |δB z /B z,0 | < 1.0 × 10 -4 in the diffusion region.
The fluctuations of magnetic field are rather small as shown in Figure 4 (right column), while the electric field shows very strong fluctuations (see the left column in Figure 4). Moreover, the fluctuations of δE x are found to be larger than δE y , δE z by an order of the magnitude. The small fluctuations of magnetic field can be understood in terms of the gradient of E x to occur along the x-axis, which shows from Faraday's law, zB/zt −∇×E that E x would not change the magnetic field, as demonstrated by δB pictures in Figure 4. Figure 4, the electric field fluctuations propagate along the z-and y-directions rather than along the x-direction following the trajectories of accelerated particles. The E z component represents Langmuir waves oscillating at ω −1 ≈ 1.3 × 10 -7 s, which is close to the electron plasma frequency ω pe for the plasma density of 10 12 m 3 accepted in this simulation (Siversky & Zharkova, 2009).

Reconnection With Multiple Magnetic Islands
As a result of the simulation setting described in section 2, we present simulations for four different times up to 32Ω − ci 1 when the reconnection reaches the maximum rate similar to Muñoz & Büchner (2018) and the turbulence is stabilized, as shown in Figure 1, that achieved later in time because our current sheet is twice thicker (d s 0.5d i ). There are multiple small magnetic islands formed at the start, which are later merged into the large island in the left across the periodic boundary and two smaller islands on the right-hand side as shown in the density and energy Owing to the periodic boundary conditions at both the ends of the z-axis, the simulation domain represents the RCSs with a chain of magnetic islands, rather than a single X-nullpoint geometry with open exhausts. The energy distributions of electrons at t 24, 32Ω −1 ci (Ω ci is the ion gyrofrequency) show clear asymmetry of particle distributions with respect to the midplane, due to the presence of a strong guiding field. The accelerated particle beams of the same charge gain the two-peak energy distributions that naturally trigger two-stream instabilities leading to the formation of either Langmuir or Bernstein waves depending on the locations where these kinetic instabilities are generated (Siversky & Zharkova, 2009;Muñoz & Büchner, 2016. It has to be noted that our model thickness of the 3D current sheet is twice the thickness used by Muñoz & Büchner (2018), but it has a much smaller beam density n b , or the plasma β, inside the diffusion region. This explains the occurrence of kinetic turbulence in our simulations while it does not appear for the current sheet with the thickness used for the simulations by Muñoz & Büchner (2018).

Suppression of Kink Instability
The reconnection process is shown to be weakly affected by the kink instability at a later time, as evidenced in the isosurface of the electron energy distribution in Figure 5. The distributions are similar in the different x − z planes along the y − direction.
If the guiding field is weak and polarization electric field is weak as well, the reconnecting magnetic fields would be strongly perturbed by turbulence as reported previously Egedal et al., 2012). For example, in the B g 0 case, we observed a twist of the magnetic flux ropes in the simulation box caused by kink instability after the same running time shown in Figure 5B. However, with the increase of the guiding field and the polarization electric field induced by separated electrons and ions, the twists are suppressed shown in Figure 5A.
Thus, the locations and the sizes of magnetic islands in different x − z planes would change, which makes it hard to make statistical analysis depending on the distance from the X-nullpoint on different x − z planes along the y − direction. Therefore, to concentrate on the turbulence other than kink instability, we should stick to the cases with a strong guiding field (b g 1), to avoid this complication.

Evaluation of Generated Turbulence
In our simulation, the ion-scale magnetic islands were formed during magnetic reconnection events as shown in Figures 4A-H. The size of the largest magnetic island reached ∼ 36d i after t 32Ω −1 ci in Figures 4G, when the reconnection reaches the maximum rate and the turbulence is stabilized. Thus, it allows us to study the plasma turbulence developed in the downstream >15d i from the X-nullpoint. As described in section 3.2.2, a strong guiding field (b g 1) is implemented to suppress the out-of-plane kink instability and to keep only the turbulence induced by accelerated particles in the geometry quasi-similar on each x − z plane.
It allows us to get statistical results of turbulence power spectrum collected in the full 3D simulation box including 64 grid points along the y-direction. The isotropized 1D power spectra, similar to the one proposed by Franci et al. (2017), are calculated in the 2D Fourier x − z-plane and averaged/ summed over the y-direction. The power spectra of electric (magnetic) fields of the whole box are measured at t 32Ω −1 ci as |E| 2 (k) (|B| 2 (k)) in the Fourier space from the whole 3D simulation region and presented in Figure 6, where k stands for the wavenumber in the reconnection plane.
In this model, the wave-number spectrum of the magnetic field formed a quasi-stable range from kd i 1 down to above kd e 1. A least-square fitting of |B| 2 (k) ∝ k α over this range indicates the slope α ≈ − 2.7 suggesting that at this moment there is quasi-stable turbulence built up. Hence, in this large 3D simulation box, the turbulent magnetic field power spectrum in the RCS formed a steady spectral slope ∝ k −2.7 near the ion inertial length, and a steeper cascade at electron scales at t 36Ω −1 ci . This is The power spectrum of the electric field drops significantly at the spatial scale close to the electron inertial scale (the solid line, k d e (n 0 ), and dashed line, k d e (n b ), on the right side of the spectra are calculated from the RCS density and background density). This suggests that during the selected time the large-scale turbulent structures are quasi-stable. It looks like the dominant fluctuations in the whole region have rather long periods (or low-frequencies, ≪Ω ce ), which are produced by ion beams, while the spectra show that the electromagnetic energy is strongly damped at the electron characteristic spatial scale (see Figure 6).
Also, in the simulations obtained by Muñoz & Büchner (2018) the 1D turbulence about the X-nullpoint obtained along the z-direction has spectral indices varying in time, which can be explained by stochastic acceleration of particles near X-nullpoint (Zharkova & Gordovskyy, 2004;Dalla & Browning, 2005;Wood & Neukirch, 2005). We understand this shifting index can be caused by the fact that the "bump-it-tail" positions in the velocity spectra of accelerated transit particles near X-nullpoint are constantly changing  and so does the turbulence, which this beam produces. In contrast at the time of maximum reconnection rate in Muñoz & Büchner (2018) the accelerated particles of the same charge (transit and bounced) gain the maximal energy close to the critical one that causes quasi-stable turbulence with noticeable power-law distribution in the wavenumber domain.

Phase Space Distributions
Now let us consider the final reconnection configuration with the two large magnetic islands separated by the X-nullpoint and explore with instant virtual spacecrafts the turbulence generated in the three locations A, B, C within the magnetic island (A), close to its edge (B), and close to X-nullpoint (C) in the current sheet x − z plane shown in the upper plot of Figure 7.
To establish a link between the turbulence and accelerated particles in the locations of these points, let us examine the changes of accelerated particle characteristics in the associated plane x − y perpendicular to the current sheet plane shown by the vertical lines in the upper plot of Figure 7 in the locations of points A and B. This gives a complete 3D presentation of the current sheet, and shows that the accelerated particles have very specific trajectories in the magnetic topology of a current sheet. In the bottom row of Figure 7 we present the particle velocity distributions in the FIGURE 5 | Upper plot: Isosurface of the electron energy distribution (the 35% contour of the max energy) for a strong guiding field (b g 1) in the simulation box of Figure 1 at t 28Ω −1 ci . Bottom plot: Isosurface of the electron energy distribution after the same running time from a similar simulation using no guiding field, e.g. b g 0.
Frontiers in Astronomy and Space Sciences | www.frontiersin.org November 2021 | Volume 8 | Article 665998 x − y plane, e.g. the x − v y phase space for both ions (bottom left) and electrons (bottom right) along the direction perpendicular to the reconnection midplane at points A and B far away from the X-nullpoint.
From the phase space analysis we can speculate that the kinetic turbulence is mainly generated by accelerated particle beams, which are later found to evolve into a phase-space hole indicating their breakage: this happens at the distance from the particle entrance in an RCS of about 7d i for electron beams and at the distances about 12d i for proton beams, where d i is an ion inertial depth. This was consistent with the previous numerical findings for simulations in different reconnecting regimes (Drake et al., 2003;Muñoz & Büchner, 2016) and the observations in the Earth's magnetotail (Khotyaintsev et al., 2010).
The particle distributions demonstrate clear non-Maxwellian features in electron beam distribution shown in location B in Figure 7C: at z 15d i (or Δz ∼ 7d i away from the main X-nullpoint). There are clearly seen two beams at the distance x ≈ 3.5d i : one with lower velocities and another one moving with much higher velocities while revealing a clear fragmented structure. In addition, there are electron holes formed in the phase space between x − 1.5d i to 1.5d i , which can be triggered by the beam-driven lower hybrid instability discussed in section 3.3.
As the inspecting spacecraft moves deeper into the magnetic island to location A, there is also the perturbation in the ion phase space found at z 10d i (or Δz ∼ 12d i away from the X-nullpoint) in Figure 7B, with the three quasi-parallel arcs located in the region between x 0 to 2d i and a very bright blob of very energetic protons located at x 0 representing the different groups of the ion beams formed during acceleration. At this instance there were no electron beams at location A, because the electron beams dissipated at the distance 7d i closer to location B ( Figure 7C), so there should be only the proton ones present  at the point A and any turbulence generated in this location has to be produced by proton beams and their interaction with the ambient plasma (Kucharek et al., 2000;Gomberoff et al., 2002). There are no any clear ion holes in the phase space, but these few arcs are found to quickly disappear further in the downstream of the beam that suggests the ion beams become scattered by the plasma turbulence generated by them that is discussed in section 3.3. Therefore, the particle velocity distributions suggest that accelerated electron or ion beams move away from the X-nullpoint until gaining the critical energy to break from this current sheet. The accelerated ions and electrons form different types of two-beam velocity distributions at different regions of the current sheet, thus producing different types of instabilities (Buneman, 1958;Kucharek et al., 2000;Gomberoff et al., 2002;Siversky & Zharkova, 2009;Muñoz & Büchner, 2018).

Frequency Analysis
Now let us study the plasma turbulence introduced by the beam instabilities using electric and magnetic fluctuations in the frequency domain.

Wavelet Analysis
After we identified the instability signals in the particle phase space, let us utilize the discrete wavelet transform, which is a powerful tool to analyze time-series data collected by a pinpoint in the domain (Farge, 1992). The signals at different grids along the y-direction were transformed to the wavelet power spectra using Morlet wavelet for the simulation domain and time up to 80 Ω ci . The turbulent fields were approximated by a short-time Fourier transform using a sliding Tukey window with an appropriate overlap. Then the results were averaged along the direction of the out-of-plane y-axis and presented at the instances in the positions of virtual spacecrafts located on the grid points along the y-direction at some given (x, z) coordinates (measured in the units of a proton inertial length d i ).
Then we record the fluctuations of electric and magnetic fields in the hypothetical locations of probes A, B, and C during the acceleration of particles in the RCS. The signals from different probes are separately transformed to the wavelet power spectra using Morlet wavelet. Then the results are averaged over all the probes with the same (x, z) coordinates.
The wavelet power spectra of both the electric and magnetic field components shared the similar features at the electron plasma frequency as expected from the results presented in section 3.1.2 and Figure 4. For example, Figure 8 shows the results using the data of the B x component recorded at point B (z 15d i , x 0.25d i ), where the electron holes were observed in the phase space in Figure 7C for a period of 5Ω −1 ci . Comparing the wavenumber spectra of electromagnetic fields from the whole region (section 3.2), the wavelet analysis confirmed that the dominant fluctuations have long periods (or low-frequency, ≪Ω ce ) (strips 1 and 2), which can be produced either by fast electron or by ion beams. This point we discuss further in section 3.3.2.
Furthermore, the wavelet transform revealed wide purple features in the high-frequency region. Figure 8 depicts several high-frequency signals represented by a wide purple strip 3 below and wide purple strip four above the electron plasma frequency ω ce ). Thus, the electromagnetic fields spectra, presented via the wavenumbers and via the wavelet transform, both indicate the important role of electrons in plasma turbulence developed in the given location B of the current sheet between its X and O-nullpoints.

Frequency Spectra of Electromagnetic Fields
We assume that the virtual spacecraft was placed simultaneously at the three different locations: A, B, and C in Figure 7 with the selected points C → A being further away from the X-nullpoint. The selected turbulent magnetic fields are collected in the surveyed boxes of the size of ΔL x ( 0.2d i ) × L y × ΔL z ( 0.2d i ) surrounding the selected points in Figure 7. The values of turbulent fields were averaged in space and time over 5Ω −1 ci using the Fourier transform. Now let us explore the resulting turbulent components of electro-magnetic fields, B and E, in every grid point of the selected locations (A, B, C) by projecting them onto the background field B m0 . This will allow us to get the parallel and perpendicular components of the turbulent field and to evaluate more accurately the turbulence nature in these locations. Note that the distributions presented in Figure 7 are taken from the lefthand side of the X-nullpoint. They are the same as those found at the similar distances on the right-hand side because the model is symmetric with respect to the X-nullpoint. The results are presented in Figure 9 for the parallel (left column) and  Figure 7) at z 15, x 0.25 (in the units of the ion inertial depth d i ) of the time series of B x components, using Morlet wavelet. Note that X and Z as in other plots are measured. The solid dark curve encloses the regions of >95% confidence. By using as a base the X-axis of the frequency spectra shown in Figure 9, the lower-hybrid frequency ω lh can be roughly drawn just above the period of 2 3 (in the units of ω −1 pi ) where the two strongest lower-hybrid frequency strips 1 and 2 (marked by yellow and red colors) are occurred at the initial times ≤36ω pi . There are also high-frequency strips detected between the electron gyrofrequency Ω ce (between 2 0 ω −1 pi and 2 1 ω −1 pi ) and the electron plasma frequency ω pe (near 2 −2 ω −1 pi ): the wide purple strip 3 of the highfrequency turbulence is located below the period of 2 −1 ω −1 pi , while the another purple strip four of this high-frequency turbulence is detected above the electron plasma frequency, just below the period mark of 2 −3 ω −1 pi .
Frontiers in Astronomy and Space Sciences | www.frontiersin.org November 2021 | Volume 8 | Article 665998 perpendicular (right column) components of the turbulent electric and magnetic fields.
In the sub-high-frequency region, Ω ce < ω < ω pe , we found several distinct spikes in all the turbulent fields at three locations marked by blue, purple, and yellow curves. In fact, there are the two small peaks occurring at higher amplitude turbulence at the frequencies below the lower-hybrid frequency, which are specifically well seen in E ⊥ , and the another two stronger peaks appearing at lower amplitude turbulence in parallel and perpendicular electric and magnetic fields at the frequencies between Ω ce and ω pe . Considering that the periodic boundary condition along the z-axis stands for simulating a chain of magnetic islands, it suggests that the magnetic island pool is fulfilled with these electromagnetic fluctuations above Ω ce . Furthermore, both high-frequency fluctuations of δE and δB are mainly perpendicular to B m0 . In the very high-frequency part (≥ω pe ), we first noticed that the perpendicular electric field E ⊥ at f > ω pe is damped significantly as it moves away from the X-nullpoint. In the other words, these high-frequency waves represented by E ⊥ are only observable near X-nullpoints (points B and C), which are also clearly seen in the wide purple patterns (strips 3 and 4) shown in the wavelet plot at these frequencies (see Figure 8).
This high-frequency turbulence is likely to be generated by two-beam instability of electron beams with "bump-in-tail" distributions in the vicinity of X-nullpoint producing Langmuir waves with the wavelength of 2 m (or 2d i in the current setting) and a speed of propagation of (1.7-2.0)· 107 m/s (or about 0.07c) with the period of 1.5·10 -7 s (close to ω −1 pe ) as reported for current sheet parameters in the solar corona in section 4.5 of Siversky & Zharkova (2009). Although, as one can observe from Figure 7C, in some locations electron beams start moving across the magnetic field lines producing, thus, Bernstein waves that are well reflected in the peaks of the perpendicular components of the turbulent fields. Both types of these plasma waves (Langmuir and Bernstein) contribute to the significant peak of high-frequency turbulence seen as in parallel so in perpendicular components. We believe that significant contribution to the broadband kinetic turbulence can appear from the electron shear flow instability suggested by Muñoz & Büchner (2018) which contributes to the perpendicular components of the turbulent electro-magnetic fields.
The most puzzling features in the current evaluation are in the low-frequency part: right below Ω ce , we found a large enhancement in the amplitude of B ⊥ (and a spike in E ) in the point A. Further down in the lower frequency region, the FIGURE 9 | The spectra of different E and B components at selected points (marked in the corresponding colors in Figure 7) as functions of the frequency (normalized to ω pi ): B , E , B ⊥ , E ⊥ with respect to the local mean magnetic field in 3D. The characteristic lower-hybrid frequency ω lh , electron gyro frequency Ω ce , and electron plasma frequency ω pe are labeled as vertical dotted lines.
Frontiers in Astronomy and Space Sciences | www.frontiersin.org November 2021 | Volume 8 | Article 665998 amplitudes of B , B ⊥ , and E ⊥ are much larger over a wider range. The small bump near ω lh (especially in the parallel electric fields) near point A at z 10, x 1, measured in the units of d i , where ω lh is the lower hybrid frequency represents the lower hybrid waves. Since in this location A we recorded only a very intense proton beam shown in Figure 7B while electron beams in the vicinity of point B were broken and formed an electron hole as shown in Figure 7C; it is safe to assume that at this instance the turbulence in point A is generated by ion/proton beams (Kucharek et al., 2000;Gomberoff et al., 2002). There is a noticeable increase of the turbulence close to the lower-hybrid frequency in the parallel components at point A and in the perpendicular components in points A and B (see Figure 9) in the parallel B in the left top plot and perpendicular electric field E in the right bottom plot). The lower-frequency turbulence is also seen in point B shown in the wavelet plot in Figure 8 as very bright strips 1 and 2 that could be driven the field-aligned drifts of highly accelerated electrons (Drake et al., 2003). This turbulence is seen in locations B and C revealing initially a growth of parallel turbulence and strong levels of oblique lower hybrid (LH) waves at later times (for which we recorded the turbulence) coinciding with a substantial parallel electron acceleration. In low-β plasmas with intense parallel currents and both with or without parallel E fields, LH waves are shown to grow even for electron distributions stable to the parallel Buneman instability, or the electron two-stream instability, and to accelerate electrons parallel to B very rapidly (McMillan & Cairns, 2006;Fujimoto & Sydora, 2008). This instability may be seen as the oblique limit of the ion acoustic and Buneman two-stream instabilities at the location where electrons beam eventually fully dissipates ( Figure  7A, point B).
Moreover, Fujimoto (2014) has shown that the intense electron beams can trigger the electron two-stream instability (ETSI) and the beam-driven whistler instability (WI). The ETSI generates the Langmuir waves, while the WI gives lower hybrid waves. This is, we believe, what is observed in the perpendicular components of turbulence in locations A and B as shown in Figure 9, right column, where strong intense accelerated beams propagate (see Figures 7B,C).
As shown in the b and c plots of Figure 7, the particle densities in these points A and B have well recorded inhomogeneities of particle densities clearly seen in Figures 7B,C, which could attribute to the generation of whistler waves in the region near these points as suggested by Zudin et al. (2019). This suggestion is also confirmed by studies of McMillan & Cairns (2007) showing that in plasmas with low beta (as we use in our model) the most unstable mode is not occurring at parallel propagation, but may be at intermediate and very oblique angles. The simulations (McMillan & Cairns, 2007) demonstrate that the very oblique lower hybrid (LH) waves can also arise. The oblique whistler waves are sometimes observed at the lower hybrid frequency in thin current sheets in the heliosphere (Zhou et al., 2009;Artemyev et al., 2016).
Also for point A one can also add generation of the rightpolarized resonant instability by very intense proton beams (Kucharek et al., 2000;Gomberoff et al., 2002). In addition, a kinetic branch of Kelvin-Helmholtz instability can be also enhancing the plasma turbulence near the lower-hybrid frequency since we clearly detected in locations B and A shown in Figure 7 the flows of protons traveling from the X-nullpoint to the O-nullpoint.
These turbulent electro-magnetic field enhancements near lower-hybrid frequency f ≈ ω lh , f < Ω ce , and at higher frequencies Ω ce < f < ω pe are also consistent with the dark horizontal stripes in the wavelet power spectrum shown in section 3.3.1. Evidently, by splitting the electromagnetic fluctuations into the parallel and perpendicular directions, we managed to identify the differences between these striped signals in the frequency analysis, which also appeared in the wavelet analysis reported in section 3.3.1. This allows us to assume that the detected turbulence signals could be the real features.

DISCUSSION AND CONCLUSION
In this article we investigate kinetic turbulence generated by accelerated particles in a reconnecting current sheet (RCS) with X-and O-nullpoints and explore the kinetic turbulence spectra in the wavenumber and frequency domains. We consider reconnection in a thin current sheet with 3D magnetic field topology using 3D particle-in-cell (PIC) approach and carry out the simulations or magnetic reconnection affected by tearing instability. In this simulation we set a larger 3D simulation domain, in which the magnetic reconnection generates two large magnetic islands each ∼ 32d i long. A strong guiding field B g is implemented to suppress the out-ofplane kink instability and to keep the geometry quasi-similar on each x − z plane. It allows us to get statistical results by averaging the data collected from the 64 grid points along the y-direction.
We reiterated our previous findings (Siversky & Zharkova, 2009; that during a magnetic reconnection in the presence of a guiding magnetic field, the particles of the same charge drifting into the RCSs from the opposite boundaries would gain different energies, higher for the transit particles and lower for the bounced particles. As a result, the high-energy accelerated particles of the same charge form non-Maxwellian distributions with the "bump-in-tail," which leads to Buneman instability (Buneman, 1958) or the electron two-stream instability, and generates the observed turbulence (Jaroschek et al., 2004;Siversky & Zharkova, 2009;Drake et al., 2010;Muñoz & Büchner, 2016. The turbulent magnetic and electric fields generated in the RCS gathered in the large 3D simulation box at the time of t 36Ω −1 ci reveal the turbulent power spectra in the wavenumber space to have a steady spectral slope ∝ k −2.7 near the ion inertial length, and a steeper cascade at electron scales, which is consistent with the other 3D PIC simulations of kinetic turbulence (Muñoz & Büchner, 2018;Li et al., 2019) and analytical estimations (Boldyrev et al., 2013;Loureiro & Boldyrev, 2017).
The characteristic waves produced by either electron or proton beams can be identified from the energy spectra of electromagnetic field fluctuations in the phase and frequency domains and compared with the particle energy gains. We selected the specific point inside the simulated 3D current sheet close to X and O-nullpoints to explore the frequencies of generated turbulence in these particular locations. We inspected the phase space of accelerated particles at this selected time, and identified the two regions with clear non-Maxwellian distributions: close to the X-nullpoints related to drift instabilities produced by accelerated electrons and away from X-nullpoints related either to drift instabilities produced by ions.
From the phase space analysis we gather the kinetic turbulence and speculate that it can be generated by accelerated particle beams seen in these locations. These beams are later found to evolve into the phase-space hole indicating their breakage: this happens at the distance of about 7d i from the particle entrance in an RCS for electron beams and at the distances of about 12d i for proton beams, where d i is the ion inertial depth. This demonstrated that in some locations of current sheet the turbulence can be generated by accelerated electron beams, while in others by proton beams. In addition, there is electron-ion hybrid instability, the kinetic branch of Kelvin-Helmholtz instability, which can also enhance the plasma turbulence near the lower-hybrid frequency since there are clearly detected flows of proton/ions traveling from the X-nullpoint to the O-nullpoint, This was consistent with the previous numerical findings for simulations in different reconnecting regimes (Drake et al., 2003;Muñoz & Büchner, 2016) and the observations in the Earth's magnetotail (Khotyaintsev et al., 2010).
To explore the kinetic turbulence in more detail, we distinguish the parallel and perpendicular components of the electric and magnetic turbulent fields (Boldyrev et al., 2013;Loureiro & Boldyrev, 2017) that reveals different levels of turbulence in the presence of a strong magnetic field. By analyzing the changes in the electric and magnetic fields in the frequency domain at different locations, we can connect non-Maxwellian features in the particle phase space with distinct fluctuations of turbulence. This frequency analysis of the generated turbulence was carried out inside the simulated current sheets: close to X-nullpoint (point C), far away from X-nullpoint (point B) and inside O-nullpoint (point A). The frequency analysis was also supported by Morlet wavelet analysis carried out in point B over the timescale of 80 Ω −1 ci . The particle distributions in points A-C clearly demonstrate non-Maxwellian features in particle distributions, e.g. the electron beam distribution in location B in Figure 7C at z 15d i (or Δz ∼ 7d i away from the main X-nullpoint) reveals two beams at the distance x 1.5: one with lower velocities and another one moving with much higher velocities while revealing a clear fragmented structure. In addition, there are electron holes formed in the phase space between x −1.5d i to 1.5d i , which can be triggered by the beam-driven lower hybrid instability. Also we show that in point A inside the magnetic island there are a few proton beams observed with arc-type structure and a break in the flow that can also produce a welldefined turbulence.
The electron beams introduced high-frequency electromagnetic fluctuations above Ω ce , which were observed in Figure 9 in the frequency spectra of the turbulence generated by beams in the surveyed points (B-C) shown in Figure 7 and also confirmed by the two wide purple strips below and above the electron plasma frequency seen clearly in the wavelet spectra in Figure 8 calculated in point B.
These rapid signals appear as distinct spikes near the highfrequency tail of the power spectra of electric and magnetic fields in Figure 9. These fluctuations are spread from the electron gyro frequency to the electron plasma frequency. This highfrequency turbulence is likely to be generated by two-beam (Buneman), or two-beam instability, of electron beams with "bump-in-tail" distributions in the vicinity of X-nullpoint as indicated by some other simulations (Siversky & Zharkova, 2009;Muñoz & Büchner, 2018) producing Langmuir waves. Although, as one can observe from Figure 7C, in some locations electron beams start moving across the magnetic field lines producing the enhanced ultra-high frequency fluctuations in the E ⊥ component, or Bernstein waves (Bernstein, 1958;Gusakov & Surkov, 2007). The similar signals were found in the inflow region close to the X-nullpoint by Lapenta et al. (2020).
Such high-frequency harmonics above Ω ce have been recently discovered by the MMS satellites near the electron diffusion region in the magnetopause (Dokgo et al., 2019). On the other hand, Li et al. (2020) reported the signals in E ⊥ and B ⊥ power spectra peaks at the harmonics of nΩ ce , where n 1, 2, 3, . . . near an electron diffusion region in the magnetotail and they were attributed to the electron Bernstein waves. One difference in the observation is that ω pe /Ω ce ≈ 27 in the magnetosphere, which keeps those two signals well separated. But this ratio is much lower in most PIC simulations including ours (ω pe /Ω ce 15), so we could not distinguish them clearly.
While in location A deeper into the magnetic island there is seen perturbation in the ion phase space at z 10d i (or Δz ∼ 12d i away from the X-nullpoint) in Figure 7B, with the three quasi-parallel arcs located in the region between x 0 to 2d i and a very bright blob of very energetic protons located at x 0 representing the different groups of the ion beams formed during acceleration. These few arcs are found to disappear quickly further in the downstream of the beam that suggest the ion beams become scattered by the plasma turbulence. Thus, the ion beams would also be quickly suppressed by two-stream instabilities. The difference between the electron and ion phase space suggests that to understand the full picture of plasma turbulence due to magnetic reconnection, it requires the simulation size to be much bigger than the diffusion region (Eastwood et al., 2018;Zhang et al., 2019).
Although, there is a noticeable increase of the turbulence close to the lower-hybrid frequency in the parallel components at point A and in the perpendicular components in points A and B (see Figure 9, parallel B in the left top plot and the perpendicular electric field E in the right bottom plot). As shown in the lower plots of Figure 7, the particle densities in these points A and B have well-recorded inhomogeneities of particle densities clearly seen in Figures 7B,C. The lower-hybrid waves can be generated by two-stream instabilities as shown in the energy distribution of Figure 7B (Papadopoulos & Palmadesso, 1976;Fujimoto & Sydora, 2008;Zhou et al., 2014;, or due to the strong density gradient near the separatrices and in the outflow (Drake et al., 2003;Scholer et al., 2003;Divin et al., 2015;Zudin et al., 2019).
In the current simulation the lower-hybrid waves are clearly seen in both the frequency and wavelet analysis applied to the gathered kinetic turbulence. The wavelet power spectrum showed that the low-frequency fluctuations at the lower-hybrid frequency have largest amplitudes and, thus, dominate in the region. These turbulent electro-magnetic field enhancements near f ≈ ω lh , f < Ω ce , and Ω ce < f < ω pe are well consistent with the bright yellow and red stripes in the wavelet power spectrum shown in section 3.3.1.
The field-aligned drifts often drive instabilities (Drake et al., 2003) revealing a growth of parallel propagating turbulence initially and strong levels of oblique lower hybrid waves at later times coinciding with substantial parallel electron acceleration (Fujimoto & Sydora, 2008). In low-β plasmas with intense parallel currents and both with or without parallel E fields, LH waves are shown to grow even for electron distributions stable to the parallel Buneman instability and to accelerate electrons parallel to B very rapidly (McMillan & Cairns, 2006). This instability may be seen as the oblique limit of the ion acoustic and Buneman instabilities (McMillan & Cairns, 2007). The low-frequency waves in the current model dominate the turbulence in the regions located further away from the X-nullpont (points A and B) since accelerated particle beams become more intense (Fujimoto & Sydora, 2008) and amplitudes of the fluctuations are increased near the lower-hybrid frequency (Rogers et al., 2000).
This suggestion is also consistent with the other study (McMillan & Cairns, 2007) showing that in plasmas with low beta the most unstable mode is not occurring at parallel propagation, but may be at intermediate and very oblique angles that are observed in the perpendicular components of turbulence in locations A and B shown in Figure 9. Evidently, by splitting the electromagnetic fluctuations into the parallel and perpendicular directions, we managed to identify the differences between these striped signals, confirming them to be the real features since the oblique whistler waves are sometimes observed in thin current sheets (Zhou et al., 2009;Artemyev et al., 2016).
Also, further investigation is required of the kinetic turbulence generated in reconnecting current sheets with different magnetic field topologies and scenarios of reconnections and their links to the specific acceleration paths of the ambient particles dragged into a current sheet with a given magnetic field topology. This dual approach to investigation of kinetic turbulence combining investigation of accelerated particle paths and distributions with the turbulence they can generate can help to uncover more accurately the mechanisms for the generation of kinetic turbulence during magnetic reconnections and its effect on accelerated particles and the whole reconnection process.
In summary, we have identified the plasma turbulence in the RCS with magnetic islands and linked the characteristic fluctuations to the non-Maxwellian distributions of particles in the phase and frequency spaces. The observed waves are found to vary as a function of the distance away from the X-nullpoint. The high-frequency perpendicular fluctuations damp quickly out of the electron diffusion region, while the lower-frequency lowerhybrid (possibly whistler) waves are developing because of the streaming instabilities generated by two electron or two proton beams.
Identifying these characteristic signals in the observation could indicate the existing scenarios of local particle acceleration during their passage through magnetic reconnection regions in the solar wind. These results can be potentially beneficial for the in-situ observations of RCSs near the Sun obtained with the Parker Solar Probe, which has already detected some reconnection sites during its first encounter (Phan et al., 2020).

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.