On the Presence and Thermalization of Cold Ions in the Exhaust of Antiparallel Symmetric Reconnection

Using fully kinetic 2.5 dimensional particle-in-cell simulations of anti-parallel symmetric magnetic reconnection, we investigate how initially cold ions are captured by the reconnection process, and how they evolve and behave in the exhaust. We find that initially cold ions can remain cold deep inside the exhaust. Cold ions that enter the exhaust downstream of active separatrices, closer to the dipolarization front, appear as cold counter-streaming beams behind the front. In the off-equatorial region, these cold ions generate ion-acoustic waves that aid in the thermalization both of the incoming and outgoing populations. Closest to the front, due to the stronger magnetization, the ions can remain relatively cold during the neutral plane crossing. In the intermediate exhaust, the weaker magnetization leads to enhanced pitch angle scattering and reflection. Cold ions that enter the exhaust closer to the X line, at active separatrices, evolve into a thermalized exhaust. Here, the cold populations are heated through a combination of thermalization at the separatrices and pitch angle scattering in the curved magnetic field around the neutral plane. Depending on where the ions enter the exhaust, and how long time they have spent there, they are accelerated to different energies. The superposition of separately thermalized ion populations that have been accelerated to different energies form the hot exhaust population.


INTRODUCTION
Magnetic reconnection is a fundamental plasma process that converts energy stored in the magnetic fields to plasma energy by enabling the reconfiguration of the magnetic field topology. Cold plasma can be abundant in regions of magnetic reconnection, and impact the reconnection process in several ways (for a recent review, see Toledo-Redondo et al., 2021). Examples include, but are not limited to, mass loading Dargent et al., 2020;Tenfjord et al., 2020), introduction of an extra cold ion diffusion region (Toledo-Redondo et al., 2016a;Divin et al., 2016), modifications to the Hall physics (Toledo-Redondo et al., 2015;André et al., 2016;Toledo-Redondo et al., 2018), and plasma heating (Toledo-Redondo et al., 2016b;Toledo-Redondo et al., 2017). In the magnetotail, cold lobe plasma below a few hundreds of eV (e.g., Engwall et al., 2009) can be heated to several keV (e.g., Ergun et al., 2018), and at the dayside cold magnetospheric plasma of a few tens of eV (e.g., Borovsky and Denton, 2008) can be heated to hundreds of eV (e.g., Toledo-Redondo et al., 2016b) as a consequence of the reconnection process. However, where, how, and to what extent cold plasma is heated (or not) by reconnection is still not fully understood.
Cold ions have been observed deep within the exhaust, both at the dayside and nightside. In the nightside magnetotail, in the vicinity of dipolarization fronts, both Nagai et al. (2002) and Xu et al. (2019) found, in addition to a hot component, two cold beams counter-streaming parallel to the magnetic field. They proposed that the cold ions moved along the reconnected field lines directly from the lobes due to the topological change enabled by reconnection. The counter-streaming beams are attributed to Fermi acceleration of ions initially dwelling on flux tubes downstream of the front that are being swept up as the front expands past them. This process is also described by Eastwood et al. (2015), although the authors attributed the source population to pre-existing plasma sheet, not the lobes. Closer to the X line, Alm et al. (2018) also found that cold ions could dominate the density well inside the separatrices, and account for a significant fraction even when the magnetic field amplitude approached zero. Counter-streaming beams observed during nightside reconnection has also been attributed to acceleration by the Hall electric field (Wygant et al., 2005). This is also shown in numerical simulations (e.g., Divin et al., 2016). While some studies suggest that the Fermi and Hall mechanisms may be dominant in the far and near exhaust, respectively, their relation and effectiveness are not fully understood. For example, Drake et al. (2009) noted that the energization by the Hall electric field arose from a potential, and that the oscillatory motion along the normal direction of the current sheet did not lead to a significant net energization. Drake et al. (2009) also presented a model of the ion heating based on a scenario of counter-streaming ion populations, originating from opposite sides of the current sheet and being subject to Fermi acceleration in the exhaust. Haggerty et al. (2015) developed the model further by taking into account a parallel electric field supported by the heated electrons moving downstream, away from the X line. The corresponding electric field potential slows down ions in the frame of the outflow field lines, reducing the acceleration, and thereby the ion heating.
At the dayside magnetopause, Li et al. (2017) observed parallel jets, carried by cold magnetospheric ions, on the magnetosheath boundary of the exhaust, showing that initially cold ions could remain cold while traversing the exhaust. Toledo-Redondo et al. (2016b) showed that cold ions could be found inside the reconnection exhaust far away from the X line. Closer to the X line, they found that the cold ions were heated by waves and large electric field gradients inside the separatrix region. Graham et al. (2017) also found that cold magnetospheric ions could interact with magnetosheath ions that enter the magnetosphere through the finite gyroradius effect. The resulting ion-ion streaming instability lead to the formation of lower-hybrid waves that can heat the cold ions. Schriver and Ashour-Abdalla (1990) also suggested that cold inflow plasma in the plasma sheet boundary layer could be heated through an ion-ion streaming instability due to the cold inbound plasma and the hotter parallel ion beams commonly observed there (e.g., Eastman et al., 1986;Nagai et al., 1998). These findings suggest that the separatrices, and the proximity to the X line, may play a role in how cold ions are heated. Dayside observations have also indicated that cold ions may become more efficiently heated if their relative density is lower (Toledo-Redondo et al., 2016b;Toledo-Redondo et al., 2017). These findings are also consistent with those by Xu et al. (2019), who presented one case where the density of the cold populations were large enough such that the density across the front remained approximately constant, in contrast to the more commonly observed density decrease (e.g., Runov et al., 2011).
The thermal population inside the exhaust of reconnection has been studied by several authors. Using hybrid simulations, Nakamura et al. (1998) showed that ion orbit characteristics vary as a function of the distance from the reconnection line due to the change in magnetic field curvature. Based on the magnetic field curvature parameter κ 2 r B /ρ i (Büchner and Zelenyi, 1989), where r B is the radius of magnetic field curvature, and ρ i is the ion gyroradius, they defined three classes of orbits. Closer to the jet front, in the magnetic pile-up region, the radius of magnetic field curvature typically exceeded the ion gyroradii, and the inflow population could remain magnetized and formed the cold counter-streaming beams mentioned above. In an intermediate region, where the radius of curvature was comparable to the ion gyroradii, the motion became stochastic. Closer to the X line, where the ion gyroradii greatly exceeded the curvature radius, the ions followed typical Speiser orbits (Speiser, 1965), meandering in the field reversal and slowly turning towards the outflow direction over the course of several bounces. The Speiser orbits typically form distributions in the shape of half circles (e.g., Nakamura et al., 1998;Zenitani et al., 2013;Vines et al., 2017) that rotates around the current sheet normal direction with distance from the X line (Arzner and Scholer, 2001;Drake et al., 2015;Nagai et al., 2015). The relative location in the circle is related to the amount of gyroturning around the normal magnetic field, which in turn is related to the number of bounces around the neutral plane the corresponding particles has performed (e.g., Zenitani et al., 2013).
In this work, we take a new integrated look at the acceleration and formation of a thermalized exhaust from initially cold inflow ions. We focus on the distinction between the cold ions being swept up by the exhaust as it expands (e.g., Eastwood et al., 2015) and the cold ions entering the exhaust closer to the X line (e.g., Zenitani et al., 2013). Using a fully kinetic 2.5D particle-in-cell simulation, where a cold inflow is captured by the reconnection process, we address the following questions: 1) When can cold ions potentially appear deep within the exhaust? 2) When do the initially cold ions remain cold, and when are they heated? 3) How are the cold ions heated?

SIMULATION SET-UP
To investigate how the cold ions are accelerated and thermalized, we perform a fully kinetic 2.5D particle-in-cell simulation (Hesse et al., 1999) with symmetric inflow conditions and no guide field.
The simulation is initialized with a Harris sheet equilibrium of hot ions and electrons. In addition, we add cold ions and electrons to the inflow regions: Densities are normalized to n 0 , times to the inverse ion cyclotron frequency ω −1 ci , lengths to the ion inertial length d i c/ω pi , velocities to the Alfvén speed v A , and energies to m i v 2 A . These quantities are based on either or both of n 0 and B 0 : ω pi (n 0 e 2 /m i ϵ 0 ) 1/2 , v A B 2 0 /(μ 0 m i n 0 ), and ω ci B 0 /em i . The ion-to-electron mass ratio is m i /m e 100, the ion-to-electron temperature ratio of the hot species is T i /T e 5, with n 0 (T i + T e ) 0.5B 2 0 , and the electron plasma-to-cyclotron frequency ratio is ω pe /ω ce 2, giving c/v A (m i /m e ) 1/2 ω pe / ω ce 20. The half width of the Harris current sheet is L 2d i . The cold ion density is n c 0.2n 0 , and the temperatures of the cold protons and electrons are initially T ic T ec 0, where the subscript 'c' refers to the cold populations. We choose to set the initial temperature to zero to be better able to follow the evolution of the initially cold populations in phase space. The particles are initialized as three ion and three electron populations: the hot Harris sheet populations, the cold plasma originating from the north (z > 0), and the cold plasma originating from the south (z < 0). As such, we are able to trace the origin of the particles at later times when they are mixed within the exhaust. The simulation box size is 200 × 25d i divided into a grid of 6400 × 1600 cells. The boundary conditions are periodic in x and reflective in z. The outof-plane electric field at the reflective boundary (z ±12.5d i ) is E y 0, which implies the conservation of magnetic flux. To initiate reconnection, a localized perturbation is added to the center of the box.

FORMATION OF DIPOLARIZATION FRONT DUE TO DECREASE IN INFLOW DENSITY
The formation of dipolarization fronts have been attributed to three main mechanisms: jet braking (Birn et al., 2011), transient reconnection (Sitnov et al., 2009;Birn et al., 2011;Fu et al., 2013), and spontaneous formation . In this section, we show that the dipolarization fronts in our simulation form due to transient reconnection associated with the transition from (hot) dense to (cold) tenuous inflow plasma. Figure 1A shows the magnetic field and density profiles along a vertical cut (z) through the X line, at tω ci 5. We note that the deviation of the originally hot plasma n i,hot in Figure 1A from a Harris sheet density profile is due to the initial perturbation. Figure 1B shows the temporal evolution of the density and magnetic field in the inflow, 1 d i above the X line. As the denser Harris sheet plasma is processed by the reconnection process, and moved downstream of the separatrices and X line, the inflow adjacent to the X line is gradually replaced by the more tenuous and cold population, resulting in a density decrease. The magnetic field also decreases, albeit slower than n √ , and the combined effect leads to an increase in the inflow Alfvén speed v A . We note that a significant portion of the decrease of |B| at later times is due to flux exhaustion caused by the finite simulation domain. However, the main sequence of interest is the increase in v A , which is not dominated by this effect. At the same time that v A increases, we observe an increase in the reconnection rate, defined as the out-ofplane electric field E y at the X line (normalized to v A B 0 ) ( Figure 1C), closely related to the changes in v A . We therefore conclude that, consistent with mass-loading scaling, the increase in reconnection rate is due to the transition from a dense to a tenuous inflow plasma.
The increase in reconnection rate leads to the formation of two flux pile-up regions propagating away from the X line. This is seen in Figure 1D where we plot B z (x, t) at z 0. The thicker black line shows the location of the main X line. The thinner black lines are equipotential contours of the magnetic vector potential A y , defined by the in-plane magnetic field (B x , B z ) through B ∇ ×A, and show the motion of field lines in the outflow.

FIGURE 1 | (A)
Magnetic field |B| and densities n i,hot , and n i,cold in the inflow along a vertical cut through the X line at tω ci 5. (B) |B|, n i,hot , n i,cold , n tot , and v A in the inflow 1d i above the X line. (C) Reconnection rate, here defined as the out-of-plane electric field E y at the X line. (D) Magnetic field B z around the equatorial plane as function of time and x. The thicker black line shows the location of the X line, while the thinner black lines show equipotential contours of A y . Newly reconnected field lines appear at the X line and propagate outwards. (E) Magnetic field |B| and densities n i,hot , and n i,cold in the outflow along a horizontal cut through the X line at tω ci 120. Note that the colors of the lines are consistent throughout the figure, and we therefore only define the legend the first time a quantity appear.
Frontiers in Astronomy and Space Sciences | www.frontiersin.org September 2021 | Volume 8 | Article 730061 Newly reconnected field lines appear at the main X line and adjacent to islands, and propagate outwards. The islands are seen as localized regions of B z reversals, numbering five in total, with the first one appearing at tω ci ∼ 85. The slopes of the lines indicate the propagation speeds of the field lines. When the reconnection rate increases, we can see how newly reconnected field lines appear at a greater rate and propagate outwards at a higher speed, such that they stack up and form two regions of stronger magnetic field B z . In magnetotail nomenclature, these flux pile-up regions are commonly referred to as dipolarizing flux bundles (DFB), and the leading edges of the two DFBs as dipolarization fronts (DFs) (e.g., Liu et al., 2013). In the magnetotail, the tailward propagating front is also referred to as an antidipolarization front (e.g., Li et al., 2014). At tω ci 120, these fronts have reached x ≈ 70 and x ≈ 135. In Figure 1E, we can see how the fronts roughly separate the initially hot and cold plasmas. From this, we conclude that the dipolarization fronts can be considered as compressed versions of the initial current sheet edge. The field lines that reconnected faster, and went on to form these flux pile-up regions, were also the ones that separated the hot from the cold plasma. This is why the initially cold plasma end up deep within the exhaust, all the way to the dipolarization fronts. Whether or not this initially cold plasma can remain cold depends on the heating mechanisms.
In the following, we will mainly focus on the time tω ci 120, where the two flux pile-up regions have expanded to x 70 and x 135. In Figure 2, we show the density, pressure, and temperature of the initially hot and cold ions, respectively. The gray contour lines show the in-plane magnetic field. The main X line is located at x 109, with a second and third X line at x 100, and x 91. These three X lines straddle two magnetic islands that are centered at x 95, and x 102. In order to more easily compare the hot and cold ion contributions to density and pressure in the exhaust, we have saturated the color scales. In the equatorial plane, the cold ion density exceeds the hot ion density for x ≳ 71, while the cold ion pressure exceeds the hot ion pressure for x ≳ 72 in the left exhaust. While the hot temperature exceeds the cold temperature in large parts of the exhaust, their effect on the dynamics inside the exhaust is negligible due to their low density. The initially cold ions are heated to temperatures comparable to the initial temperature of the hot Harris sheet ions, where the initial Harris sheet temperature is given by T i 0.5B 2 0 /n 0 (1 + T e /T i ) ≈ 0.42. These values are also comparable to the levels of ion heating predicted by Drake et al. (2009), T i m i v 2 0 /3 ≈ 0.16 − 0.33, where v 0 ≈ 0.7-1.0 is the speed of the exhaust magnetic field lines that varies slightly throughout the FIGURE 2 | Density, pressure, and temperature of initially hot and cold ions at tω ci 120: The initially cold plasma dominate the density and pressure at x ≳ 71 and x ≳ 71, respectively. While the hot temperature exceeds the cold temperature in large parts of the exhaust, their effect on the dynamics inside the exhaust is negligible due to their low density. The initially cold ions are heated to temperatures comparable to the initial temperature of the Harris sheet ions.
Frontiers in Astronomy and Space Sciences | www.frontiersin.org September 2021 | Volume 8 | Article 730061 exhaust (as shown by the slope of the black lines in Figure 1D). At the same time, the hot ions are heated to slightly higher values. Overall, the density, pressure, and temperature of initially cold ions are fairly structured. Concentrations of cold ions are found inside the islands, and adjacent to the front (extending along the magnetic field toward the north and the south, as well as along the equatorial plane towards the X line). The density cavities located at the separatrices are associated with high speed electron flows toward the X line (not shown). To investigate the kinetic structure responsible for the density, pressure, and temperature profiles, we will in the next section study the ion distributions within the exhaust.

KINETIC STRUCTURE OF ION DISTRIBUTIONS WITHIN THE EXHAUST
In this section we investigate the structure of the ion distributions within the exhaust. Figure 3 shows the reduced distributions of 2,4] between the two dipolarization fronts at time tω ci 120. The reduced distributions are integrated over the remaining velocity dimensions, e.g.
At this time, the initially cold ions dominate the density and pressure within the exhaust (see Figure 1E and Figure 2), which is why we do not include the initially hot population in the further analysis. The distributions are sampled over boxes spanning a spatial domain of Δx × Δz 0.5 × 0.5, i.e. x 75 → x 75 ± 0.25, and z 4 → z 4 ± 0.25. The dotted vertical lines show the location of the separatrices, the dashed vertical lines at z 0 show the locations of the DFs, and the black solid lines show the corresponding components of v E×B . The still cold populations are identified by higher phase space densities that occupy a relatively small velocity interval. The initially cold ions that have become heated are seen as (relatively) lower phase space densities that occupy larger velocity intervals.
Heated ions are observed throughout the exhaust. Cold ions are most prominent in the inflow region at z 2 and z 4, but are also found throughout the exhaust, at all z. In an intermediate region, inside the separatrices, the inbound ions (v z < 0) are affected by the convergent electric field associated with the separatrix density cavities. At z 4, this region is clearly seen between x 118 and x 127 for the right exhaust ( Figure 3C). At z 2, this region has extended over a larger range along x (x 113-124), and the inbound ions have become more thermalized ( Figure 3F). At z 0, this region maps to the intermediate exhaust ( Figure 3I). In f(v z , z 0), we see counter-streaming populations throughout the exhaust. The coldest populations are found in the vicinity of the X lines and the DFs (where they are superposed on a hot population). Outbound cold ions (v z > 0) in f(v z , z 2) are only observed in the immediate vicinity of the DFs (70 < x < 75 in left exhaust). This means that only a portion of the cold beams observed at z 0 are actually leaving the equatorial FIGURE 3 | Reduced distributions of initially cold ions at tω ci 120, at (A-C) z 4, (D-F) z 2, and (G-I) z 0. The vertical dashed lines show the location of the separatrix, which at z 0 maps to the X line. The solid line shows v E×B . We do not show these in panels (H) and (I) because they contain B x , which is small there. While the cold ions are most apparent in the inflow, they are also present in the outflow. In an intermediate region, just inside the separatrices, the cold incoming population is affected by the dynamics at the separatrices, and becomes partly thermalized. This is perhaps most prominently seen in f(x, v z ).
Frontiers in Astronomy and Space Sciences | www.frontiersin.org September 2021 | Volume 8 | Article 730061 region with their temperature relatively intact. In the rest of the exhaust (closer to the separatrices and X lines), the outbound plasma at z ≥ 2 is strictly thermal.
In f(v y , z 0), we see a distinct transition in the speed of the cold ion population at around x 82 in the left exhaust (with a corresponding transition in the right exhaust). At x < 82, the cold ions are observed at v y > 0, while at x > 82 the cold ions are observed at v y < 0 (with gradually increasing v y towards the X line). As explained by Divin et al. (2016), the initially v y < 0 is acquired through a process where the inbound ions are first accelerated by the Hall electric field E z leading to an enhanced flow v z towards the equatorial plane. This v z is turned towards −ŷ by the Lorentz force v z B x < 0. They defined the discrete transition in v y as the edge of the ion diffusion region. We note that similar behaviour is also observed inside the separatrices in the offequatorial regions. However, the transition is continuous here as v y < 0 is gradually turned towards v y > 0 by the out-of-plane electric field E y deeper within the exhaust. Closest to the DF, the cold beam seen at v y > 0 approaches v y 0. In the following sections, we will investigate the role of separatrices and magnetic field curvature in shaping the kinetic structure of cold ions inside the exhaust.

The Role of Separatrices in Ion Thermalization
Depending on where the ions cross into the exhaust the inbound ion motion v z , mentioned in the previous section, can either be on average magnetized (v ⊥ ≈ v E×B ) or demagnetized (v ⊥ ≠ v E×B ). This is illustrated in Figure 4, where we plot the reduced distributions f ic (v z ) as a function of z at x 75 and 85. The dotted line (now horizontal) again shows the separatrix location, while the solid lines show v E×B,z . As expected, the ion edge (e.g., Gosling et al., 1990;Lindstedt et al., 2009), i.e. the location at where the first reconnected ions appear (here represented by the thermal population centered on v z > 0) inside of the separatrix, is further inside the separatrix deeper within the exhaust. At x 75 and x 85, B x dominates the magnetic field in the range |z| > 1.5 and |z| > 0.5, respectively. Where B x dominates, v z constitutes a perpendicular component of the ion flow, and a deviation of the cold ions from v E×B,z signifies demagnetization. At x 75, the cold ion population follows v E×B,z until z ≈ 2 (see also Toledo-Redondo et al., 2018), while at x 85, they deviate shortly after crossing the separatrix. At x 85, the Hall electric field E z is quasi-stationary over the time it takes for an ion to cross from the separatrix to the neutral plane. For this location, we have also plotted the speed v ϕ z 2eϕ z /m i corresponding to the Hall electric field potential ϕ z − E z d z . We have grounded the potential to a region just outside the separatrix. The increase in v ϕ z seen outside the separatrices is due to wave effects formed due to the finite box size, and have no effect on the dynamics within the exhaust. Although the inbound beam is gradually thermalized between the separatrix and the neutral plane, it roughly follows v ϕ z all the way to z 0. This indicates that these ions are demagnetized, and accelerated across the magnetic field by E z as explained by Divin et al. (2016). This initial speed, governed by the potential, roughly forms the outer boundary of the thermal ion population within the separatrices at this location. This thermal population consists mainly of ions that have crossed the neutral plane at least once. We note that although the other components can also lead to demagnetization, here we chose to focus on v z to illustrate the clear difference between v ϕ z and v E×B,z at x 85.
To further illustrate, and to better understand the differences between the populations that remain cold deep within the exhaust, and the population inside the separatrices that are slightly more thermalized, we make use of test particles   [90,100,110,120]. We have divided the ions into four groups, based on the ions speed v y , and postition x at t 120; yellow (v y > 0, x ≤ 75), black (v y > 0, 75 < x ≤ 85), white (v y < 0), and green (v y > 0, 75 < x < 90). We note that the out-of-plane displacement of the ions during the integrated time interval is about 4 − 8d i . This corresponds to 0.3-0.5R E (n 0.3 cm −3 , 1d i ≈ 415 km, 1R E ∼ 15d i ), which is less than estimates of the dawn-dusk scale of plasma sheet flows derived from observations (Nakamura et al. (2004) finds a dawndusk extent of about 2 − 3R E ). The 3D structure of the reconnection exhaust in the magnetotail should therefore be large enough to accommodate the level of acceleration seen in our simulations.
The main difference between the yellow/black and the white/ green ions is where they cross the separatrices, or rather, what sort of separatrices they cross. While the separatrices, by definition, span the entire range of x in the simulation, they are highly active closer to the X line, and comparatively inactive, further away.
Here, we define active separatrices as characterized by large amplitude parallel electron flows toward the X line, associated density cavities, and electric fields, electrostatic waves propagating along the magnetic field in the direction of the electron flow, and in general stronger Hall fields. While a Hall electric field can still be observed inside the inactive separatrices, it is typically weaker, and is generally characterized by the absence of deeper density cavities, small to none parallel electron flow, and less wave activity. The white/green ions enter the exhaust at active separatrices, are slightly thermalized at the separatrices, substantially accelerated by E z , and eventually form the hot thermal population. The green ones first enter an island, while the white ones directly enter the open exhaust. The yellow/black populations enter the exhaust at inactive separatrices, or rather, the separatrices expand to encompass them. When they are caught up by the expanding exhaust region, they are picked up by the motional electric field E y , and are accelerated while remaining relatively magnetized.
In addition, the black and white ions also show distinct differences in their z distribution. In Figure 5B, we can see that the white ions overall have larger |v z | at z 0. The white ions, affected by the Hall electric field, gain higher v z , and are able to penetrate deeper into the southern side of the exhaust. This is reflected in the density n top ic . The subdivision into the yellow and black populations is mainly based on their behaviour around, and after having crossed, the neutral plane. While the time period of the test particles is not large enough to properly display their different behaviours, the density structure is. At all times, there is a clear cutoff in n top ic at about z −2 in the intermediate exhaust. Closer to the dipolarization front, the density extends below z −2. At t 125, the black ions are situated at the edge of the density cutoff and are just about to reflect. In contrast, the yellow ions, located approximately at the same z, continue downward. In the next section, we will show that whether ions are reflected or not is largely due to the level of magnetization.

The Role of Magnetic Field Curvature in Ion Thermalization
To illustrate how the ion magnetization changes from directly adjacent the front, to a distance away, we take a closer look at the magnetic curvature k B b · ∇b, whereb B/|B| is the magnetic field unit vector. Büchner and Zelenyi (Büchner and Zelenyi, 1989) defined the curvature parameter κ 2 r B /ρ i , where r B | k B | −1 is the curvature radius and ρ i v i,⊥ /ω ci is the ion gyroradius. For κ ≫ 1, the particle motion is adiabatic and the magnetic moment is conserved. This implies that the change in pitch angle as the ion approaches the neutral plane is reversed when the ion leaves the neutral plane on the other side. Hence, the shape of an ensemble is relatively well preserved after transmission. For κ → 1, the motion becomes stochastic, which enables non-reversible pitchangle scattering during the neutral sheet crossing. Figure 6B shows the 2D map of log 10 r B , while Figure 6C shows the value at z 0. For reference, in Figure 6A we also show the 2D density map of initially cold ions originating from the top. In the equatorial plane, closer to the front, due to the flux-pile up and more dipolarized magnetic FIGURE 5 | Test particle positions at five different times. The initial conditions for each particle are taken at t 120, from the locations in phase space shown in (A, B). (C-G), and (H-L) shows the particle position for five different times, overlaid E z and n top ic , respectively. The subdivision into yellow/black/white/green populations are partly based on the initial conditions, and further explained in the text.
Frontiers in Astronomy and Space Sciences | www.frontiersin.org September 2021 | Volume 8 | Article 730061 field, the magnetic field is stronger, the magnetic curvature is smaller, and consequently the radius of curvature is larger. Since the gyroradii vary significantly over our ensemble of ions due to differences in particle speed, so will κ. In Figure 6C, we have plotted the speed v i,⊥ κ −2 r B ω ci , for κ 1. The speed corresponding to κ 1 (yellow line) peaks at r B ω ci ≈ 2.5 at about x 72. At about x 75, r B ω ci ≈ 1. Figures 6D-H show 2D reduced distributions f (v x , v y ) at z 0, and x [72, 74, 76, 78, 80] (white squares in Figures 6A,B). At z 0, the magnetic field is dominated by B z , and f(v x , v y ) is therefore the distribution in the plane perpendicular to B, such that v i,⊥ v 2 x + v 2 y . The circles overlaid the distributions show v i,⊥ κ −2 r B ω ci , with κ 1, based on the local magnetic field, and centered on v E×B (marked by the dot '·'). The bulk velocity is shown with the cross ('×'). The radii of the circles are given by the local value of v κ 1 i,⊥ (yellow line in Figure 6C). Ions outside (inside) the circle correspond to κ < 1 (κ > 1), or equivalently ρ i > r B (ρ i < r B ). Hence, we expect ions that are well within the circles to follow adiabatic motion and to experience only limited pitch-angle scattering. In contrast, ions that are located close to, and outside, the circles are expected to experience significant pitch angle scattering. Closest to the front (x 72), v κ 1 i,⊥ is relatively large, and even encompass the majority of the thermal population. At x 74, the cold population is still well inside v κ 1 i,⊥ , but a large part of the thermal population is outside. At x 76, 78, 80, the drift of the cold populations, and even their thermal spread becomes comparable to, and eventually exceeds, v κ 1 i,⊥ . Based on this, we expect ions crossing the neutral plane closer to the DF to experience less pitch-angle scattering, and to be able to remain relatively cold while traversing to the other side of the current sheet.
To investigate the evolution of the ion distributions around the neutral plane, we plot the 1D reduced distributions f top ic along two select field lines (marked as white lines in Figures 6A,B). In i,⊥ ≈ 0.5 at z 0. The distributions are plotted as a function of distance along the field line s, where s 0 corresponds to z 0, and, because the direction of the field line runs from north to south, s > 0 (s < 0) corresponds to z < 0 (z > 0). The corresponding components of v E×B are shown as black lines.
We start by looking at the field line located closer to the front ( Figures 6I-K). At about s −10, the cold populations slowly drift with the magnetic field. In f top ic (v x ) a hotter outbound (leaving the neutral plane) beam is also observed at about v x −1. At around s −5, we see clear evidence of ion trapping in f top ic (v x ) and f top ic (v z ). These structures in phase space correspond to the small density enhancements seen at around x 70 and z 5 in Figure 6A. We will discuss the corresponding waves further in Section 4.3. Here we merely note that these structures contribute to the cold ion thermalization already before they enter the central current sheet. The cold inbound ions roughly follow v E×B until s ≈ − 2. Here the rapid rotation of the magnetic field fromx toẑ leads to the deviation of v z from v E×B,z , such that the perpendicular motion becomes parallel as part of the Fermi acceleration (slingshot) process (Northrop, 1963). In f top ic (v y ), the cold population is seen to oscillate a few times close to v y 0. This is due to the gyration around the magnetic field of the relatively well magnetized ions. The level of magnetization is also seen in f top ic (v x ) at z 0, which in this location represents a perpendicular component. Close to z 0, the cold component has a speed v x ≈ v ⊥ ≈ v E×B,x . After the neutral plane crossing (z < 0, s > 0), some of the cold ions are reflected at s ≈ − 2, while some manage to continue outward. While both E z and v x B y play minor roles in the acceleration alongẑ, the main responsible force is − v y B x . The ions that are reflected have a speed v y > 0 and therefore experience an upward force − v y B x > 0, since B x < 0 in the south. The ions that manage to continue outward are the ones that have crossed into v y < 0 during the gyromotion and experience a downward force − v y B x < 0. We also note that we can still identify the cold component in f top ic (v x ) at s 10. In the other two components, the beam is thermalized at about s 7.
We now look at the field line located further away from the front, in the intermediate exhaust ( Figures 6L-N). With exception for the ion trapping seen closer to the front, the behaviour is approximately similar until s −3. One important difference between the two field lines is the behaviour in f top ic (v y ) inside s ≈ − 3. Closer to the front, the cold population was seen to oscillate close to v y as part of their gyromotion. Further from the front, the cold population is instead experiencing a continual increase in v y until s ≈ 2, whereafter it no longer is identifiable in any of the components. The positive v y leads to a reflective (upward) force − v y B x > 0. As a result, the entire cold population is reflected, explaining the density cutoff seen at z ≈ −2 in Figure 6A. The absence of v y -oscillations in this location is related to the level of magnetization. The bounce distance (∼ 4d i ) becomes comparable to the gyroradius such that complete gyromotion is not possible. During the reflection, the beam is also thermalized due to the enhanced pitch-angle scattering.

The Role of Streaming Instabilities Deep Within the Exhaust
Inside the exhaust, we observe two off-equatorial regions (one in the north and one in the south) where electrostatic solitary waves (ESW) form (Figure 7). The ESWs form at around z ±4 and grow both in amplitude and perpendicular (to B) extent as they propagate downstream with the magnetic field at low parallel speeds. To investigate how they form we examine a timestep (tω ci 115) before they have grown to their full strength. Figures  8A-D shows a few quantities plotted as a function of the distance s along a magnetic field line marked in Figure 8E (the same field line is also marked by the thicker black line in Figure 7). s 0 marks the neutral plane z 0, s > 0 is in the region z < 0 (since the magnetic field direction is from north to south), and s < 0 is in the region z > 0. Figure 8A shows the parallel electric field E in which the ESWs are seen as bipolar spikes. Figure 8B shows the reduced distribution function of the initially cold ions as a function of v . The inbound ions move at low parallel speeds further out, and are gradually accelerated as they approach and cross the neutral plane. This acceleration is due to Fermi acceleration, where the initial v y gained by E y is first turned towards z by the magnetic force − v y B x , and thereafter towards − x by v y B z . We note that the final speed gain is consistent with what we expect due to Fermi acceleration: v after −v before + 2v DF , where the speed of the front v DF ≈ 0.5. Oscillations in the cold component indicate that they are affected by the parallel electric field. The results of the wave-ion interaction, for a closely located field line, are also seen in Figures 6I,K where we only plot the ions originating from the north.
The level of plasma mixing, defined by the phase space density ratio of cold ions originating from the top to all cold ions: f mix FIGURE 8 | Thermalization due to ion-acoustic waves within the exhaust. Panels (A-E) show quantities sampled along a magnetic field line marked by the black line in panel (E). They are plotted as a function of the distance along the field line s, where s 0 is the equatorial plane and s > 0 is in the region z < 0 (since the magnetic field direction is from north to south).  Figure 8C. Red (purple) indicates that all ions originate from the top (bottom) while yellow indicates mixed origin. The thin black lines in this panel are added as reference and shows the contour levels log 10 f ic −1. While most of the mixing seem to occur close to the equatorial plane, there is also clear evidence of mixing associated with the ESWs. While the inbound ions can remain cold deep within the exhaust, we see in the electron distribution that the electrons are heated ( Figure 8D). We note that this color scale is saturated to highlight the region of phase space affected by the waves. In the location of the waves, electrons form clear hole-like structures (e.g. immediately to the left of the vertical dashed line). The black line show the electron bulk speed, with a finite antiparallel (parallel) drift to the north (south) side of the equatorial plane.
In Figure 7, we can see that the ESWs convect downstream with the magnetic field. In the frame of the magnetic field, their motion is approximately field-aligned. To investigate what kind of instability may generate the waves that grow into ESWs we therefore solve the one-dimensional electrostatic plasma dispersion equation: where Z is the plasma dispersion function (Fried and Conte, 1961), ω ps again is the plasma frequency, v ds is the bulk drift speed, v ts 2T s /m s √ is the parallel thermal speed, of population s, ω ω r + ic is the complex frequency with real frequency ω, and imaginary frequency c as the growth rate. The input is taken from the location marked by black squares in Figure 7A and Figure 8E, and a dashed vertical line in Figures 8A-D. In the dispersion relation solver, we include six populations seen in Figure 8F Figure 8H. Positive growth rates are found for 0 < kd i ≲ 60, with a peak growth rate c max ≲ 0.5 at k max d i ≈ 16. This wave number corresponds to a wavelength λ ≈ 0.4. At k max , ω ≈ 6, with the corresponding phase speed v ph ∼ 0.4. The wavelength (width of bar) and phase speed (vertical location of bar) at max growth rate are plotted in Figure 8B, and correspond well to the observed wave parameters. The wave growth is attributed to an ion-acoustic instability due to the cold population at v ∼ 0, and the heated drifting electrons. We note that in the absence of the electrons, an ion-ion drift instability would also grow. However, due to the presence of the electrons, it is stabilized. In the next section, we will investigate how the ions that were demagnetized at the separatrices subsequently form the hot thermal exhaust.

FORMATION OF THE THERMALIZED EXHAUST
In this section we investigate the formation of the thermal population within the exhaust. The distribution f ic (x, v y ) in Figure 9 (see also Figure 3) exhibits discrete structures in phase space. One of the most distinct structures is a cold beam ranging between v y ≈ −0.5 at x ≈ 82 and v y ≈ 0.5 at x ≈ 92 (with corresponding structures in the right exhaust). As explained in Section 4, this beam is formed by a combination of forces E y > 0 and v z B x < 0 acting on the ions during their first inbound leg towards the neutral plane from the inflow regions.
Where v y < 0, v z B x is dominating, while where v y > 0, E y has been dominating. Divin et al. (2016) showed that this striation in phase space followed contours of constant canonical momentum p y m i v y + qA y . Echoes of this initial striation are seen in both exhausts. While the first echo is towards higher v y , the striations successively rotate counter-clockwise (clockwise) in FIGURE 9 | Reduced ion distribution f(x, v y ) at z 0, and t 120. The lines show contours of constant canonical momentum p y m i v y + qA y , which well explains the overall phase space structure (x, v y ) of the thermalized exhaust. Since v y (t 0) 0, ions with the same p y originated from the same flux tube, defined by a given A y . Ions with v y > 0 (< 0) are currently downstream (upstream) of their original flux tube.
Frontiers in Astronomy and Space Sciences | www.frontiersin.org September 2021 | Volume 8 | Article 730061 the left (right) exhaust, such that the range of v y eventually covers a large range of positive and negative values. The first echo is comprised of ions that have crossed the neutral plane once, have been turned around by B x , and returned to z 0 again. During this motion, they have been accelerated to slightly higher v y , by the reconnection electric field E y . Once they have acquired a v y > 0, they will start to become turned downstream by v y B z < 0. These two forces will jointly move the initial striation to higher v y and lower x. By plotting contours of constant p y in Figure 9, we show that also the higher order (in terms of bounces) striations approximately follow contours of constant p y . This kind of trajectory was first described by Speiser (1965). That is, the thermal exhaust is formed by a superposition of Speiser orbits with varying initial conditions, related to where they enter the exhaust. As explained by Divin et al. (2016), the change in v y of the zero order striation is due to spatially different acceleration by the Hall electric field E z . Closest to the X line, the integrated electric field is smaller than further downstream, and therefore the initial v z leading to v y through v z B x is larger further away from the X line. We also note that each time the ions cross the neutral sheet, they are scattered in the highly curved magnetic field, smearing out and further blending the superposed structures. By integrating test particles in the dynamically changing fields, we find that p y is constant to a good approximation. Initially at t 0, all of these ions had v y 0, such that p y qA t 0 y . Since a given value of A y defines an in-plane field line (B x , B z ), the value of p y informs us of what flux tube the ion were originally located on.
For the values of p y shown in Figure 9, the corresponding flux tubes were initially located at z ± [4.3, 4.7, 5.1, 5.4, 5.9, 6.5, 7.0, 7.6, 8.2]. Correspondingly, as the simulation progresses, the deviation of A y from its initial value A t 0 y , ΔA y A y (t) − A t 0 y , tells us how much the ion is currently deviating from its original flux tube. We can not measure ΔA y in nature, however, the argument can be extended to v y where Δv y v y (t) − v t 0 y . In our simulation Δv y v y (t), since v t 0 y 0. In nature, comparing, for example, the thermal and bulk speeds of lobe ions to those of exhaust ions, we also typically expect v y (t) ≫ v t 0 y , such that Δv y ≈ v y (t). Therefore, ions with v y > 0 (ΔA y < 0) are currently located downstream of their original flux tube, while ions with v y < 0 (ΔA y > 0) are currently located upstream of their original flux tube. To intuitively understand this, it is perhaps easiest to consider an ion impinging upon the X line from the top or the bottom inflow regions, with small velocity v x . Such an ion would first be accelerated by the Hall electric field E z , across the magnetic field, leading to ΔA y > 0. The finite v z would be turned toward −ŷ by the force v z B x . The increase in ΔA y is offset by the decrease in v y (Δv y < 0), such that the change in canonical momentum Δp y mΔv y + qΔA y ≈ 0. The ion will thereafter dwell for some time in the X line region, performing meandering motion in the B x field reversal, and being accelerated by E y , increasing v y . During this time, when the ion remains close to the X line, the magnetic field will convect past it, such that it is successively moved to more upstream field lines (ΔA y < 0). The decrease in ΔA y is offset by the change in v y (Δv y > 0), again keeping p y constant. The obtained v y will be turned toward ± x by the force v y B z . The larger v y , the larger the force, and resulting v x . An ion that were significantly delayed, instead gained a larger v y , and will be better able to catch up to its original field line A t 0 y (Drake et al., 2009). We also recall that A y is defined by both B x −z z A y and B z z x A y . The information of both magnetic field components in the force terms v z B x and v y B z are thus contained within A y , which also helps us understand the connection between the described forces, and why the structure of the phase space is so well described by p y .
We can also think of ΔA y and Δv y in terms of diffusion. As described above, an ensemble of ions that dwells in the vicinity of the X line, being accelerated by E y , sees the magnetic field lines move past them. Or, in other words, the magnetic field diffuses across the ensemble of ions. Further downstream, when v y is turned to v x , the ions are catching up to their original field lines, and the diffusion is reversed. The primary diffusion is done by the electric field, which performs work on the ions. Since the reversed diffusion is done by magnetic forces, no work is done, and a net energy transfer between magnetic field and plasma is realized.

DISCUSSION AND SUMMARY
In this study, we have investigated the occurrence and heating of initially cold ions during symmetric magnetic reconnection. We showed that ions that enter the exhaust at active separatrices with large gradients and electric fields become heated and thermalized more efficiently than ions that enter the exhaust at inactive separatrices lacking large gradients and electric fields. These findings are consistent with observations at the dayside magnetopause (Toledo-Redondo et al., 2016b). The ions that enter the exhaust at inactive separatrices are typically picked up as the exhaust expands into them, both in the outflow (x) and normal (z) directions. While Eastwood et al. (2015) suggested that ion beams observed close to the dipolarization front originated from the pre-existing plasma sheet, Birn et al. (2017) concluded that they originated from the lobes. We find that the cold beams can originate from an extended range of z's that may include, but also extends well outside, the pre-existing hot current sheet population. The cold ion beams closer to the front originate on lower z's, while those further away originate at higher z's. This is illustrated both in Figure 5 (c.f. original locations of yellow and black ions) and Figure 9 where the cold beam (here seen in f ic (v y )) extend over a range of p y . Recall that p y indicated the original z location of the field line/particle, and that lower p y corresponded to a field line initially located further out in the inflow. How far this range of z's extends likely depend on the stage of reconnection. The further the front has propagated, the larger range of z's can be picked up without being significantly thermalized at the separatrices. In addition, while cold counter-streaming beams could be observed close to the equatorial plane throughout the exhaust, the beams that entered the exhaust at inactive separatrices could remain colder in all three components. Also, closer to the X line, the v z was gained mainly through acceleration by E z , while further downstream, the v z gain was through Fermi acceleration (mediated by E y ).
The separatrix environment, and therefore also the ion thermalization at the separatrices may be sensitive to the reduced mass ratio we employ in this simulation (e.g., Lapenta et al., 2010). However, we note that observations also show divergent electric field at the separatrix electron flow channels. The integrated potential associated with these fields in observations could be comparable or higher than the inflow energy of the cold ions . We also note that out-of-plane instabilities active at the separatrices, such as the lower-hybriddrift instability, or ion-ion cross-field drift instabilities do not develop in our 2D simulation. Based on previous observations at the magnetopause, such instabilities would also act to enhance the ion thermalization at the separatrices (e.g., Toledo-Redondo et al., 2017;Graham et al., 2017).
In our simulation, electrostatic field-aligned waves were also found to impact the ion thermalization in the off-equatorial region deep within the exhaust. The waves were produced by an ionacoustic instability due to the cold inbound ion population and the drifting electrons. Field-aligned electrostatic waves has also been found in observations in the vicinity of dipolarization fronts relatively close to the neutral plane . These waves were likely generated by cold proton, and potentially oxygen, beams. Field-aligned electrostatic waves have also been identified in regions closer to the separatrix in the magnetopause boundary layer at the dayside (Steinvall et al., 2021). Similar to our results, they found the waves to be driven by an ion acoustic instability. Altogether, these results shows that waves may play a role in the thermalization of cold ion populations in many different subregions during reconnection.
We showed that the cold ions that entered the exhaust downstream of active separatrices could remain cold all the way to the neutral plane, and beyond. After having crossed the neutral plane, ions in the intermediate exhaust, reflected back at about a distance of 2d i from the neutral plane after having crossed it. Closer to the dipolarization fronts, due to the higher level of magnetization, the cold ions could also continue outward beyond z ±2d i . At the neutral plane, the incoming cold populations also experienced pitch angle scattering due to magnetic field curvature. As such, at a given location, the inbound beam (approaching the neutral plane) is colder than its outgoing counterpart. This is consistent with observations from the vicinity of dipolarization fronts in the magnetotail where Xu et al. (2019) observed cold beams propagating along the magnetic field. The beam approaching the neutral plane was colder than the beam leaving it. Overall, the region alongẑ in which we expect to observe counter-streaming cold ion beams is more limited in the intermediate exhaust than closer to the front.
Since the majority of the outflow is dominated by the species that were initially cold, the division into a hot and cold ion diffusion region such as explored by Divin et al. (2016), and typically invoked at the magnetopause (Toledo-Redondo et al., 2015), is not applicable here. At the dayside, the multiple ion diffusion regions are typically a result of the multiple inflow populations, e.g. cold magnetospheric plasma and warm magnetosheath plasma (Toledo-Redondo et al., 2015). Although the heated ions seen inside the separatrix in Figure 3 do not follow v E×B , this is not an indication of multiple ion diffusion regions, it is merely a superposition of inbound and outbound populations of the same origin. However, multiple ion diffusion regions could exist in the tail if there are multiple species in the inflow (e.g., oxygen ions (Wygant et al., 2005;Tenfjord et al., 2018), if there is an extended region of overlap of hot and cold plasma in the transition from the plasma sheet to the lobes, if there are overlapping regions with ions of different origins (for example the ionosphere and solar wind), or if there are significant north-south plasma asymmetries, for example due to mantle plasma (e.g., Artemyev et al., 2017) or asymmetric ionospheric outflow (e.g., Glocer et al., 2020).
In this study, we have investigated the acceleration, thermalization, and occurrence of cold ions in the exhaust of symmetric antiparallel magnetic reconnection. The overall temperature increase was comparable to the predictions based on the counter-streaming ion model by Drake et al. (2009). However, we found that several processes conspire to form the eventually hot thermalized exhaust ( Figure 10). The initially cold populations are thermalized due to interactions at the separatrices, around the neutral plane, and due to waves inside the exhaust. At the same time, the individual populations are accelerated (mainly) by the out-of-plane electric field E y . At a given position, we can find populations that have been accelerated to different energies, depending on where they entered the exhaust. Together, the thermalization and superposition of ions that have been accelerated to different energies form the hot exhaust population. The acceleration by E y lead to super-Alfvénic speeds. If we assume an inflow magnetic field of 20 nT, and a plasma sheet density of 0.3 cm −3 (recall that the Alfvén speed in our simulation is based on the inflow magnetic field and central Harris sheet density), the Alfvén speed is v A ≈ 800 km/s. An ion accelerated to 1v A and 2v A thus has an energy of about 3300 eV and 13400 eV, respectively. As for the overall FIGURE 10 | Illustration of thermalization process of cold ions. Individual populations are thermalized by one or several processes, for example separatrix activity, magnetic field curvature scattering, or plasma waves. The hot exhaust is formed by a superposition of these individually thermalized populations that are accelerated to different energies.