Turbulence and Particle Acceleration in Collisionless Magnetic Reconnection: Effects of Temperature Inhomogeneity across Pre-reconnection Current Sheet

Magnetic reconnection is an important process in various collisionless plasma environments because it reconfigures the magnetic field and releases magnetic energy to accelerate charged particles. Its dynamics depend critically on the properties of the pre-reconnection current sheet. One property in particular, cross-sheet temperature inhomogeneity, which is ubiquitous throughout the heliosphere, has been shown to increase reconnection outflow speed, energy conversion efficiency, and secondary island formation rate using two-dimensional particle-in-cell simulations. Here we expand upon these findings, considering two cases with a long, thin current sheet, one with homogeneous temperature and one with inhomogeneous temperature across the current sheet. In the inhomogeneous temperature case, numerous secondary islands form continuously, which increases current sheet turbulence (well-developed cascade power spectra) at large wavenumbers. Current density, energy conversion, dissipation, and acceleration of high-energy particles are also enhanced relative to the homogenous temperature case. Our results suggest that inhomogeneous temperature profiles, which are realistic, need to be incorporated into studies of turbulence and particle acceleration in collisionless magnetic reconnection.


Introduction
Magnetic reconnection is a process of magnetic field line rearrangement in plasmas during which magnetic energy is explosively converted to plasma kinetic energy through acceleration and/or heating of charged particles (e.g., Birn & Priest 2007;Yamada et al. 2010;and references therein). It occurs within current sheets, sheet-like plasma regions with discontinuous magnetic fields and consequently strong current density. Such sheets abound in heliophysical and astrophysical systems, such as Earth's magnetosphere (e.g., Petrukovich et al. 2015), the solar corona (e.g., Ko et al. 2003), compact objects (e.g., Melia & Kowalenko 2001), and pulsar magnetospheres (e.g., Arons 2012; Cerutti et al. 2015). Their properties (thickness, plasma beta, magnetic configuration) control the intensity, duration, and extent of magnetic reconnection and the ensuing energy conversion (e.g., Lembège & Pellat 1982;Pritchett 2005).
Using two-dimensional (2D) PIC simulations, Lu et al. (2019) showed that such temperature inhomogeneity increases reconnection outflow speed and expedites secondary plasmoid (island) formation. These effects can change the structure of the reconnection diffusion region, which increases the turbulence of fields and flows there (e.g., Lazarian & Vishniac 1999;Daughton et al. 2011;Karimabadi et al. 2014;Fu et al. 2017). Faster reconnection outflow and the presence of secondary islands can also provide additional particle acceleration (e.g., Drake et al. 2006Drake et al. , 2009Oka et al. 2010aOka et al. , 2010b. Motivated by these earlier findings, in this paper, we demonstrate that turbulence and particle acceleration in reconnection also increase in the presence of temperature inhomogeneity across the pre-reconnection current sheet: in the case with strong temperature inhomogeneity, the small-scale secondary islands form not only expeditiously, but also continuously and numerously. Their evolution and interaction increase the turbulence level and the particle acceleration efficiency of reconnection in the high-k (large wavenumber) and high energy ranges, respectively.

Simulation Model
We employ a 2D PIC model. Electrons and ions are all treated as full particles, and magnetic and electric fields are updated in time by solving Maxwell's equations using an explicit algorithm. The positions and velocities of the electrons and ions are advanced in the fields. The coordinate system is 2D Cartesian in the x-z plane within a rectangular domain, The initial equilibrium consists of a Harris current sheet and background plasma without any initial perturbations. The initial particle number density is given by where n 0 is the Harris current sheet peak density, L is the current sheet half-width, and n b is the background density. The initial magnetic field, which is in the x direction, is described by where B 0 is the asymptotical magnitude of the magnetic field. The initial electron and ion velocity distributions are Maxwellian for both the Harris population and background population. The Harris population has a temperature T 0 , and the background population has a temperature T b , from which the combined temperature profile, The background density is n b =0.2n 0 , and the half-width is L=0.5d i , where d i =c/ω pi is the ion inertial length defined by n 0 . The ion-to-electron mass ratio is m i /m e =100, and the speed of light is c=20V A , where V A is the Alfvén speed evaluated using B 0 and n 0 . The initial ion-to-electron temperature ratio is constant, where Ω i0 =eB 0 /m i is the unit ion gyrofrequency. The size of the simulation domain is L x =409.6d i and L z =25.6d i , with grid numbers N x =8192 and N z =512. About 4.9×10 9 particles are used in each simulation run. Periodic boundary conditions are adopted in the x direction; and perfect conductor boundary conditions are used in the z direction. To examine the effects of temperature inhomogeneity, two cases with different normalized values of background temperature, Case A with T b =T 0 and Case B with T b =0.1T 0 , are considered. According to Equation (3), Case A has a homogeneous temperature profile across the current sheet, and Case B has an inhomogeneous temperature profile for which the temperature decreases significantly from the current sheet's center to its boundaries.

Turbulence Structure
Figures 1, 2, and 3 show the structures of magnetic field, electric field, and current density in Cases A and B at a representative time, t 70 i0 1 = W -. In both cases, magnetic reconnection apparently results in multiple X-lines and large (several tens of ion inertial lengths, d i , wide) magnetic islands (Figures 1 and 2). Upon closer inspection, however, it is evident that Case B also results in numerous smaller (several d i wide) magnetic islands near the X-lines, such as the one at Figure 2. We refer to these larger and smaller islands as primary and secondary islands, respectively. Both primary and secondary islands can be identified using the bipolar signature of magnetic field B z in Figures 1(a), 2(a), 3(a), and (g). The reason why secondary islands are formed numerously in Case B (with inhomogeneous temperature) is the temperature inhomogeneity increases the reconnection outflow speed (Lu et al. 2019), and this fast outflow speed favors the streaming tearing mode instability (Wang et al. 1988;Hoshino & Higashimori 2015), which leads to secondary reconnection that is responsible for the formation of secondary islands in Case B. A quadrupolar structure in the out-of-plane magnetic field, B y , forms at each primary X-line. Inside the primary islands, B y is manifested as wavy structures caused by the Weibel instability Schoeffler et al. 2013). Strong electric fields (even stronger than those at primary islands) form at the secondary islands in Case B. The out-of-plane electric field, E y (normalized to V A B 0 ), indicates the reconnection rate. At the primary X-lines, it has a value of ∼0.1, typical for fast reconnection (e.g., Birn et al. 2001). At both ends of the primary islands, the outflowing reconnected magnetic flux piles up, resulting in a stronger E y than that at the primary X-lines, as described theoretically by Lu et al. (2013a). The secondary islands are embedded in fast outflows (V x ) from the reconnection at the primary X-lines. Therefore, their electric field E y is very strong and bipolar because of their bipolar B z and their fast motion (E y ≈V x B z ). The in-plane electric field components, E x and E z , are formed by charge separation, therefore this electric field is also called the Hall electric field. The Hall electric field is much stronger at the secondary islands than at the primary islands, suggesting that there is a much stronger Hall effect or a much larger charge density at secondary islands.
The region near each X-line, where the electron frozen-in condition is violated, is also called the electron diffusion region. In this region, the out-of-plane current density j y is carried primarily by electrons. The j y at the secondary islands, also carried primarily by electrons, is much larger than that at the X-lines. For example, as shown in Figure 3(j), the secondary island at x≈314d i has a very intense j y ≈15en 0 V A (the value at the initial time is 2en 0 V A ). Note that the electron current density j ey is negative in the primary islands but positive in the secondary islands, as shown in Figures 1(g) and 2(g). This, in addition to the spatial scale, can be used to distinguish secondary islands from primary islands in spacecraft in situ observations. The negative j ey inside the primary islands can be explained by the E z ×B x electron drift (in the fine-structured E z , see Figures 1(e) and 2(e)), which the demagnetized ions cannot follow (analogous to the negative j ey just ahead of dipolarization fronts in Earth's magnetotail; see Lu et al. 2016). The current density j y at both ends of the primary islands, which peaks at about 5en 0 V A (Figure 3(d) and (j)), is self-consistent with the sharp B z spatial gradient along x at these locations (Figures 3(a) and (g)).
To quantify the energy conversion, j·E, and dissipation, j·E′ (where E′=E+V i ×B), by reconnection, their profiles along z=0 are plotted for Cases A (Figures 3(e) and (f)) and B (Figures 3(k) and (l)), respectively. Because of the secondary islands, the energy conversion power and dissipation have larger magnitudes in Case B than in Case A. In Case B the dissipation has a peak value of ∼3 at x≈314d i in the secondary islands (Figure 3(l)), about 10 times higher than the value in the primary islands and several tens of times higher than the value at the X-lines. This peak value in the secondary island corresponds to a relatively large value of ∼0.8 nW m −3 in the context of the Earth's magnetotail. Regarding the total dissipation, the integral across this spike at x≈314d i is ∼7.32, whereas the integral across the highest spike in the primary islands in Case A (Figure 3(f)) is ∼3.24, showing that the secondary islands provide the strongest dissipation in the entire system. These j·E′ spikes are well resolved by at least five data points, so they are not random fluctuations. Figure 4 shows the evolution of the reconnected magnetic field's magnitude B z | | along z=0 (the center of the reconnecting current sheet) for Cases A and B. The early stage of evolution corresponds to the linear growth of the instability of a long, thin current sheet-the collisionless tearing mode (Coppi 1965;Drake & Lee 1977), which saturates at about t 30 i0 1 = W -. Although the linear growth results in small magnitude perturbations of B z (and E y ), these perturbations form numerous current sheet tearing points as potential seeds for subsequent nonlinear growth of fast magnetic reconnection. From t 30 i0 1 » Wto t 50 i0 1 » W -, some seed tearing points become locations of fast reconnection (primary X-lines), where B z grows quickly from B 10 2 0 to ∼B 0 during this time. About eight primary X-lines are formed in both Cases A and B, resulting in about eight primary islands between them. In Case A, after t 50 i0 1 » W -, the electron diffusion region elongates, decreasing the reconnection rate and stopping the fast growth of B z (e.g., Daughton et al. 2006). In Case B, after t 50 i0 1 » W -, numerous secondary islands (filament structures at the primary X-lines in Figure 4(b)) are formed. For example, about eight secondary islands are formed near the primary X-line at x≈200d i . These secondary islands eventually merge into the primary islands through an anti-reconnection process (e.g., Pritchett 2008, also called re-reconnection).
In Cases A and B, turbulence is generated by various mechanisms in different stages of evolution. In the early  -W -, the collisionless tearing instability in the current sheet dominates the spectral structure with k x d i ≈1 ( Figure 5), consistent with the linear theory prediction of k x L≈0.55 (e.g., Pritchett et al. 1991; here L is the current sheet half-width, and L=0.5 for our study). Later on, during the nonlinear stage, t 60 80 i0 1 = -W -, once the X-lines and magnetic islands are fully developed, a power spectrum consistent with the energy cascade in turbulence is seen in both Cases A and B. The B z power peaks at k x d i ≈0.1, corresponding to the eight primary X-lines forming in our simulation domain of length L x =409.6d i . From the early to the late stage, the peak of B z power shifts from kd i ≈1 to kd i ≈0.1, which means that 1 out of 10 tearing seeds grows into fast reconnection X-lines. In both Cases A and B, the B z power shows a power-law spectrum in the low k range, with a spectral index of −1.64 in Case A and a harder spectral index of −1.08 in Case B. Also in both cases, the spectrum has a breakpoint at kd i ≈1-2, above which (in the high-k range) the B z power spectral also follows a power law, and the spectral index becomes ≈−4.5. Figure 6 shows the power spectra of magnetic field, electric field, and current density as functions of k along z=0 at t 60 80 i0 The dominant magnetic field component is B z , so the magnetic field power is mostly contributed by the B z power (see Figure 6(a)), which governs the power of any topological changes in the magnetic field. The electric field in the low k range is mostly contributed by the power of the out-of-plane electric field, which characterizes the changes of the in-plane (reconnecting) magnetic field. In the high-k range, the electric field power is dominated by the electric field E x , which is strongest at the small-scale (high-k) secondary islands (see Figures 2(c) and 3(h)). The current density power is dominated by the j y power in almost all k ranges except the primary X-line k range, around kd i ≈0.1, in which the j x power is also large. This strong j x is produced by the difference between ion and electron reconnection outflow motions (ion outflows are Alfvénic, and electron outflows are super-Alfvénic).
The major difference between Cases A and B is that in Case B, numerous secondary islands with a spatial scale on the order of an ion inertial length d i are formed continuously, which causes spectral differences between the two cases, as seen in Figure 7. The B z power in the high-k range in Case B is higher than in Case A, even though its spectral index is very similar in the two cases. The relative difference in the electric field power is even greater than in the magnetic field power. This is because of the secondary islands and their strong bipolar E x (Figures 2(c) and 3(h)): their electric field power in the high-k range is about five times larger in Case B than A, thus affecting the total power. Because of the secondary islands and their strong j y (Figures 2(f) and 3(j)), the current density power in the high-k range is also higher in Case B than Case A. These higher powers of magnetic field, electric field, and current density show that Case B is more dynamic than Case A, especially in the high-k range, because of the continuous formation, growth, and merging of many secondary islands. For the same reason, a similar behavior is exhibited in energy conversion, j·E, and dissipation, j·E′, as seen in Figure 8  . Therefore, the secondary islands in Case B cause more efficient magnetic energy conversion and dissipation than in Case A. Figure 9 shows the increase in particle kinetic energy that is converted from magnetic energy for Cases A and B. In the nonlinear stage, once fast reconnection begins (at t 40 i0 1 = Wand thereafter), there is an abrupt increase in this kinetic energy for both electrons and ions, and it is about 15% greater in Case B than Case A by t 100 i0 1 = W -. This difference is smaller than that in the power spectral densities of j·E and j·E′ (Figure 8) because the largest spectral differences are observed in the high-k range, which contributes only modestly to the total amount of energy conversion and dissipation. However, as we shall see next, this contribution is significant for energetic particle acceleration, i.e., the energy conversion and dissipation in the high-energy range.

Particle Acceleration
Figures 10 and 11 show plasma information for Cases A and B, respectively. As shown in Figures 10(a) and 11(a), the plasma density is low at the X-lines and high inside the magnetic islands (both primary and secondary). The electron temperature is low at the X-lines and at the very center of the primary islands, but high in the magnetic ribbons that wrap around the primary islands, especially at the islands' both ends. The electron temperature in the secondary islands is high. The ion temperature is higher in the magnetic islands and lower at the X-lines. (The background temperature is much lower in Case B than in Case A, and this is from the initial setup.) Figures 10(d) and 11(d) show the charge density e n n i e r = -( )in the two cases. The charge density is negative at the X-lines, which is consistent with the bipolar Hall electric field E z pointing toward the X-lines. The charge density inside the secondary islands is also negative and even more intense than that at the X-lines. This intense, negative charge density is also evidenced by the strong bipolar E x and E z at the secondary islands (see Figures 2(c), (e), and 3(h)). The reconnection electron outflows are super-Alfvénic (Figures 10(e) and 11(e)), and the reconnection ion outflows are usually Alfvénic (Figures 10(f) and 11(f)). This kind of difference in electron and ion outflow motions forms the current density j x at each primary X-line, which provides a high j x power in the low k range, as shown in Figure 6(c). The electron and ion inflows (see V ez and V iz in the last two panels of Figures 10 and 11) also show differences in magnitude, forming the current j z in the inflow region. Figure 12(a) depicts the electron energy spectrum calculated in the entire simulation domain. By comparing the electron  energy spectra in the two cases, one can see that more highenergy electrons are produced in Case B. Figure 12(b) shows the electron energy spectrum at the X-line for the two Cases, in Regions A1 and B1 (marked in Figures 10(a) and 11(a), respectively). Because of the secondary island in Region B1 at x z d , 314, 0 i » ( ) ( ) , the electron energy spectrum differs significantly from that in Region A1 (which has no secondary island)-many more high-energy electrons are produced in Region B1 than Region A1. The ion energy spectra in the entire simulation domain and in Regions A1 and B1 have a similar feature-many more high-energy ions are produced in Case B (see Figures 12(d)), especially at the secondary island (see Figure 12(e)). The electron and ion energy spectra in primary islands, in Regions A2 and B2 (marked in Figures 10(a) and 11(a), respectively) also differ-more high-energy ions are at the primary islands (Figure 12(f)), and the high-energy electron flux is also slightly higher in Region B2 than in Region B1 (Figure 12(c)).
To further determine where these high-energy electrons are accelerated, Figure 13 plots the contours of high-energy   . In Case A, most relativistic electrons are confined to the magnetic ribbons of the primary islands, suggesting that these islands act as the major venue for high -energy electron acceleration. In Case B, relativistic electron energy flux is also high in the ribbons of the primary islands, higher than that in Case A. More importantly, the relativistic electron energy flux is even higher inside secondary islands, especially the one at x z d , 314, 0 i » ( ) ( ) . This means that the secondary islands provide a more efficient acceleration mechanism for high-energy electrons. Because of the efficient electron acceleration in the secondary islands, in the highest energy range shown in Figure 13, γ e >2, the high-energy electron energy flux integrated over the entire simulation domain is about 3.3 times higher in Case B than A. Figure 14 shows the contours of high-energy ion energy flux, , in different energy ranges. Because even the high-energy ions are not relativistic, the ion energy is normalized to another unit, the ion initial temperature kT i0 . In Cases A and B, the high-energy ions are mostly confined to primary islands. Similar to the electrons, the high-energy ion energy flux at the primary islands is higher in Case B than Case A. The secondary islands in Case B can also accelerate ions. However, both primary and secondary islands cannot confine the ions that are in very high energy ranges, such as kT 40 i i0  > , and in this high-energy range, the high-energy ion energy flux integrated over the entire simulation domain is about 4.2 times higher in Case B than Case A. Figure 15 shows the profiles (along z=0) of the highenergy electron and ion energy fluxes for Cases A and B; the magnetic field component B z is also plotted. Each negativethen-positive B z signature corresponds to an X-line, and each positive-then-negative B z signature corresponds to an island. Two peaks very close to each other (on the order of d i ) indicate a secondary island. Seven secondary islands are identified and denoted by red arrows 1-7. For the high-energy electrons, the most notable feature of the electron energy flux is that it is very high inside each secondary island in Case B, especially the one at x≈314d i , inside which the electron energy flux in the highenergy range, γ e >2 (denoted by red curves in Figure 15(e)), is about 10 −1 . Although about one to two orders of magnitude lower than that inside the secondary islands, the high-energy electron energy flux is still high at the primary islands, especially in their magnetic ribbons in both Cases A and B. In the highest energy range, γ e >2, the electron energy flux is about 10 −3 -10 −2 in the magnetic ribbons of the primary islands, and on average, it is higher in Case B than in Case A (Figures 15(b) and (e)). This may be because the energetic electrons and ions in the secondary islands (numerous secondary islands form continuously in Case B) eventually move into these ribbons after the secondary islands merge/ coalesce into the primary islands. Moreover, this merging/ coalescence process through an anti-reconnection process can further accelerate the electrons and ions (e.g., Pritchett 2008;Oka et al. 2010b;Tanaka et al. 2010Tanaka et al. , 2011Le et al. 2012;Zhou et al. 2014) and thus further increase the energetic electron and ion energy fluxes at the primary islands in Case B.
The ion energy flux is high in both the primary and the secondary islands because ions can be trapped in those islands and accelerated to high energies. Unlike the electron energy flux, the ion energy flux in the secondary islands is lower than that in the primary islands. Also note that in the very highenergy ranges (e.g., kT 40 i i0  > ), the ion energy flux is not well confined to the islands because the magnetic structure of these islands is not strong enough to trap these highly energetic ions, which are much heavier than electrons. For example, in the high-energy range, kT 40 i i0  > , the ion velocity v i is faster than about 5.7V A , corresponding to a gyroradius of at least 6d i (the maximum magnetic field at the secondary islands is about 1.0B 0 , see Figure 3(g)). This scale is comparable or even larger than the radius of secondary islands, so the high-energy ions become demagnetized at the secondary islands and cannot be trapped by them. Like that of the electrons, the high-energy ion energy flux is also higher in Case B than Case A in high-energy channels, such as kT 40 i i0  > . For example, the highest ion energy flux in this energy range is about 10 −1 in Case B and about 3×10 −2 in Case A.
Using the electron and ion velocity distribution functions shown in Figures 16 and 17, we demonstrate the kinetics and acceleration mechanisms of the electrons and ions. Figure 16(a) shows the electron velocity distribution functions at a primary X-line. Because the electrons have a Speiser-type meandering motion (Speiser 1965) near the X-line, their distribution function has two peaks in the v ez direction. During its meandering motion, each time an electron reaches z≈0, where the field is almost zero, it is accelerated by the reconnection electric field E y in the −v ey direction. Electrons with a small v ex can stay near the X-line for a longer time, so they can go through more meandering cycles and be accelerated for more times. Therefore, the electron acceleration in the v ey direction peaks at v ex ≈0 (Pritchett 2006a;Ng et al. 2011;Bessho et al. 2014). Electron velocity distribution in the outer ribbons of a primary island is shown in Figure 16(b), which indicates that the electrons have a higher perpendicular velocity v ex and v ey (the magnetic field here is mainly in the z direction), and a ring-like structure is formed in the v ex −v ey plane. These features suggest that the major acceleration mechanism for electrons in the outer ribbons is betatron acceleration caused by magnetic field enhancement. In the inner ribbons of the primary island, the electrons have a different velocity distribution, which has a high flux in both the parallel (v ez ) and perpendicular (v ex and v ey ) directions, as shown in Figure 16(c). This suggests that although betatron acceleration is also retained in the inner ribbons, Fermi acceleration mechanism plays also an important role in this region through field line shrinkage (e.g., Drake et al. 2006). Note that in the magnetic ribbons of the primary islands (Figures 16(b) and (c)), the electrons are much more energetic than at the X-lines ( Figure 16(a)). The electrons in the secondary island are the most energetic ones, as shown in Figure 16(d). Because the electrons are well confined to the center of the secondary island, they can be efficiently accelerated by the electric field structure there through the island surfing acceleration mechanism proposed by Oka et al. (2010aOka et al. ( , 2010b. The acceleration is biased toward the −v ey direction, and this is because the secondary islands are formed at the primary X-lines, where the electrons are pre-accelerated in the −v ey direction by the reconnection electric field E y . The ion velocity distribution function at the X-line, as shown in Figure 17(a), has a similar feature as the electron's-two peaks in the v iz direction, and a preferential acceleration in the v iy direction at v ix ≈0, which indicates that the ions have a similar acceleration mechanism at the X-line, as described above. The ion acceleration in the magnetic ribbons of the primary island shows a multiply striated structure, suggesting a hybrid acceleration processes therein (Figure 17(b)). Ions in the low-energy range are magnetized and thus accelerated by betatron acceleration. Ions in the high-energy range, which form the multiply striated structure at v iy ∼5V A , have a gyroradius of about 5d i (the magnetic field magnitude is about 1.0B 0 ), comparable to the half-width of the ribbon of the primary island. Therefore, these ions are no longer fully frozenin to the magnetic field lines, and their acceleration depends on how many times they enter this ribbon. Every time they enter this region, they are accelerated by the strong E y (Figures 1(d) and 3(c)) in the v iy direction, and when they re-enter this region, they are accelerated again. This type of ion acceleration during multiple entries forms the multiply striated structure at higher v iy in the ion velocity distribution (Figure 17(b)). Although ions inside secondary islands can also be confined and accelerated through the island surfing mechanism, very high energy ions cannot be confined. Inside the secondary island, ions are preferentially accelerated in the +v iy direction, which is also due to the pre-acceleration in the +v iy direction at the primary X-line before the secondary island is formed. The electron and ion velocity distributions described above demonstrate various acceleration mechanisms that are responsible for production of high-energy electrons and ions in different regions during reconnection.

Discussion
Turbulence in magnetic reconnection has recently been observed in Earth's magnetotail and in the solar wind (e.g., Eastwood et al. 2009;Huang et al. 2012;Vörös et al. 2014). Using observations of the terrestrial magnetotail and the solar flare, Hoshino et al. (1994), Fu et al. (2017), and Cheng et al. (2018 have suggested that such turbulence is related to magnetic islands (or called O-lines, plasmoids, etc.). In addition, Fu et al. (2017) have shown that energy dissipation j·E′ in turbulent magnetic reconnection is particularly strong at O-lines. Our PIC simulations present a similar picture of this turbulence and the dissipation therein. We find that in the realistic, inhomogeneous temperature case, continuous formation of numerous secondary islands results in significant amplification of kinetic-scale turbulence, and energy dissipation is very strong in the secondary magnetic islands. As shown in Figure 3(l), the energy dissipation j·E′ peaks in the secondary island at x≈314d i ; this peak value is about 2.8, corresponding to a large value of about 0.8 nW m −3 in the context of Earth's magnetotail, which is on the same order as Fu et al. (2017) Sitnov et al. (2018);and Liang et al. (2019). The quantification of the kinetic dissipation in turbulent collisionless reconnection by secondary islands needs further investigation using these new parameters.
Secondary islands, which form at reconnection X-lines, substantially modulate reconnection properties (e.g., Daughton et al. 2006Daughton et al. , 2011Drake et al. 2006;Pritchett 2006b;Eastwood et al. 2007;Wang et al. 2010Wang et al. , 2016bHuang et al. 2011;Zhou et al. 2012). In our simulation Case B (with inhomogeneous temperature), numerous secondary islands form continuously. The pronounced high-energy electron and ion energy fluxes in these secondary islands (e.g., in the secondary island 7 at t=70 i0 1 W -, as shown in Figure 15(e)) show that these islands are the most efficient particle accelerators in this system. The particle acceleration mechanism, first proposed by Oka et al. (2010a), operates as follows: the strong in-plane electric field E x and E z (Figures 2(c) and (e)) at the secondary islands balances the Lorentz force acting on the particles (v y B z and v y B x , respectively), and this allows them to "surf" along the out-of-plane direction y. As they do, they are accelerated by the out-of-plane electric field E y . This electron and ion island surfing acceleration is biased toward the −v ey and +v iy directions, respectively, because secondary   islands are formed at the primary X-lines, where electrons and ions are pre-accelerated by the reconnection electric field E y . This type of strong particle acceleration inside secondary islands is consistent with the strong energy dissipation therein (see Figure 3(l)).
As revealed by Figures 16(a) and 17(a), the electron and ion acceleration mechanism at the X-line is E y acceleration during the particles' meandering motion, a finding that is consistent with previous PIC simulations (Fu et al. 2006;Pritchett 2006a;Divin et al. 2010;Ng et al. 2011). More recently, the electron velocity distribution has been found to have a fine structure of multiple striations (Bessho et al. 2014(Bessho et al. , 2018Torbert et al. 2018). (This is not clearly seen in our simulation, most likely due to the limited number of particles used to calculate the velocity distribution functions.) Although an X-line can accelerate particles, its configuration cannot well confine the particles because the particles can easily escape in the x direction. Therefore, the high-energy electron energy flux is low at the X-lines (Figures 13-15), which is also consistent with in situ spacecraft observations in Earth's magnetosphere for most cases (e.g., Wu et al. 2015), except for the unique event found by Øieroset et al. (2001).
Electrons are accelerated preferentially in the magnetic ribbons of primary islands by betatron and Fermi acceleration mechanisms, as revealed by Figures 16(b) and (c). Betatron acceleration and Fermi acceleration are adiabatic: they take place during magnetic field enhancement and field line shrinkage, as the first and second adiabatic invariants, respectively, are conserved (e.g., Drake et al. 2006;Ashour-Abdalla et al. 2011;Fu et al. 2011Fu et al. , 2013Guo et al. 2014;Li et al. 2015;Du et al. 2018;Lu et al. 2018a). High-energy ions are also trapped and accelerated in the magnetic ribbons of the primary islands, but nonadiabatically. Every time the ions enter the ribbons, they gyrate for a half circle during which they are accelerated in the v iy direction by the strong E y (Figures 1(d) and 3(c)), and when they re-enter the ribbons, they are accelerated once again. This type of acceleration by multiple entries forms the multiply striated structure at higher v iy in the ion velocity distribution (Figure 17(b)). Note that the above particle acceleration mechanisms at the X-lines, primary islands, and secondary islands are deduced based on velocity distribution functions (Figures 16 and 17) and previous literature. In future studies, it will be interesting to use particle tracing analysis to better demonstrate these acceleration mechanisms.
The aforementioned particle acceleration mechanisms may vary in the case of guide field reconnection, but this is out of the scope of this study. In guide field reconnection, particles can be magnetized more readily by the guide field, so guidingcenter theory can be used to analyze contributions from different adiabatic acceleration mechanisms, such as Fermi, betatron, and parallel electric field (e.g., Dahlin et al. 2014Dahlin et al. , 2015Wang et al. 2016a;Lu et al. 2018a). For example, in guide field reconnection, electron acceleration in secondary islands is no longer through the nonadiabatic island surfing mechanism (Oka et al. 2010a(Oka et al. , 2010b, but through adiabatic Fermi acceleration and parallel electric field acceleration mechanisms, and when the guide field is very large, the contribution of the Fermi mechanism is negligible compared to the parallel electric field acceleration (Wang et al. 2017). The electron acceleration mechanism at the X-line also changes in the presence of the guide field: the electrons no longer follow the meandering orbits but gyrate in the guide field and are accelerated by the parallel electric field (e.g., Pritchett 2006b).
In this study, we focus on collisionless magnetic reconnection and find that its turbulence and energetic particle acceleration are very much increased by the presence of a temperature inhomogeneity, owning to the continuous formation of many secondary islands. Chains of secondary islands (plasmoids) have also been found formed in collisional magnetic reconnection, which is believed to be crucial to its fast reconnection rate based on magnetohydrodynamic (MHD) theory and simulations (e.g., Loureiro et al. 2007;Bhattacharjee et al. 2009;Samtaney et al. 2009;Pucci & Velli 2014;Tenerani et al. 2015) and later PIC simulations (e.g., Daughton et al. 2009aDaughton et al. , 2009b. These collisional MHD and PIC simulations all used a homogeneous temperature profile as their initial condition. Because the temperature inhomogeneity has significant effects on collisionless reconnection, especially on the formation of the secondary islands, their effects on collisional reconnection are potentially also important and worthy of further investigation.

Conclusions
Using 2D PIC simulations of a long, thin current sheet, we studied collisionless magnetic reconnection and its dependence on temperature inhomogeneity across the pre-reconnection current sheet and found that: 1. Long, thin current sheets, with or without temperature inhomogeneity across them, are unstable to the collisionless tearing mode. Although linear growth of this mode provides small perturbations, it forms numerous current sheet tearing points that act as seeds for subsequent nonlinear growth of fast reconnection. About 1 in 10 of these tearing points grows into fast reconnection X-lines. 2. Once fast reconnection begins, the current sheet system evolves nonlinearly: multiple X-lines and primary islands develop. This is the stage in which the temperature inhomogeneity plays an important role, causing numerous secondary islands to be formed continuously at reconnection X-lines. 3. In the nonlinear stage, the formation, evolution, and interaction of multiple X-lines and primary islands cause the current sheet to become turbulent, with welldeveloped cascade power spectra of magnetic field, electric field, current density, and energy conversion and dissipation. The numerous small-scale secondary islands, produced in the presence of temperature inhomogeneity, and their evolution and interaction make the system more turbulent, which changes the spectral structure and increases the spectral power in the high-k range. 4. Because magnetic energy is efficiently converted/dissipated to particle kinetic energy in the nonlinear stage, electrons and ions are accelerated at X-lines and primary islands through various mechanisms. By having a strong power of energy conversion and dissipation, the secondary islands, which form only in the inhomogeneous temperature case, provide an additional venue to efficiently accelerate particles, especially electrons, so the high-energy particle energy flux is much higher in the inhomogeneous temperature case.