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 (Fuselier et al., 2017; 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 κ ² = 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?

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: B x = B 0 ⁡ tanh ( z / L ) n = n 0 ⁡ cosh − 2 ( z / L ) + n c 1 2 + 1 2 tanh | z | − 2 L 0.5 L .

Densities are normalized to n ₀, times to the inverse ion cyclotron frequency ω c i − 1 , lengths to the ion inertial length d _i = c/ω _pi , velocities to the Alfvén speed v _A , and energies to m i v A 2 . These quantities are based on either or both of n ₀ and B ₀: ω p i = ( n 0 e 2 / m i ϵ 0 ) 1 / 2 , v A = B 0 2 / ( μ 0 m i n 0 ) , and ω _ci = B ₀/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.5 B 0 2 , and the electron plasma-to-cyclotron frequency ratio is ω _pe /ω _ce = 2, giving c / v A = ( m i / m e ) 1 / 2 ω p e / ω c e = 20 . The half width of the Harris current sheet is L = 2d _i . The cold ion density is n _c = 0.2n ₀, 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 out-of-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.

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 (Sitnov et al., 2013). 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-of-plane electric field E _y at the X line (normalized to v _A B ₀) (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. 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 anti-dipolarization 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.5 B 0 2 / 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 0 2 / 3 ≈ 0.16 − 0.33 , where v ₀ ≈ 0.7–1.0 is the speed of the exhaust magnetic field lines that varies slightly throughout the 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.

In this section we investigate the structure of the ion distributions within the exhaust. Figure 3 shows the reduced distributions of initially cold ions f _ic (x, v _x ), f _ic (x, v _y ), f _ic (x, v _z ) at z = [0, 2, 4] between the two dipolarization fronts at time tω _ci = 120. The reduced distributions are integrated over the remaining velocity dimensions, e.g. f(v _x ) = ∫f(v _x , v _y , v _z )dv _y dv _z . 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 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 − y ̂ 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 off-equatorial 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.

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 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 = q A y t = 0 . 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 y t = 0 , Δ A y = A y ( t ) − A y t = 0 , 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 Δ v y = − q m i Δ A y , where Δ v y = v y ( t ) − v y t = 0 . In our simulation Δv _y = v _y (t), since v y t = 0 = 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 y t = 0 , 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 − y ̂ 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 y t = 0 (Drake et al., 2009). We also recall that A _y is defined by both B _x = −∂ _z A _y and B _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.

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 (Norgren et al., 2020). We also note that out-of-plane instabilities active at the separatrices, such as the lower-hybrid-drift 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 ion-acoustic 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 (Liu et al., 2019). 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 z ̂ 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⁻³ (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 temperature seen in Figure 2F, T _i = 0.3 corresponds to about 2000 eV. These energies are comparable to average ion temperatures in the central plasma sheet (e.g., Baumjohann et al., 1989). This suggests that magnetic reconnection can heat cold lobe plasma to large enough energies to replenish the hot plasma sheet.