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

1 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 (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?

2) When do the initially cold ions remain cold, and when are they heated?

3) How are the cold ions heated?

2 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:

$\begin{array}{ccc}\hfill B_x& =\hfill & B_0\tanh (z/L)\hfill \end{array}$

$\begin{array}{ccc}\hfill n& =\hfill & n_0{\cosh }^{-2}(z/L)+n_c\left[\frac{1}{2}+\frac{1}{2}\tanh \left(\frac{|z|-2L}{0.5L}\right)\right].\hfill \end{array}$

Densities are normalized to n₀, times to the inverse ion cyclotron frequency ${\omega }_{ci}^{-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₀: ${\omega }_{pi}={(n_0{e}^2/m_i{ϵ}_0)}^{1/2}$, $v_A=B_0^2/({\mu }_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.5B_0^2$, and the electron plasma-to-cyclotron frequency ratio is ω_pe/ω_ce = 2, giving $c/v_A={(m_i/m_e)}^{1/2}{\omega }_{pe}/{\omega }_{ce}=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.

3 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 (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 $\sqrt{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_AB₀) (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.

FIGURE 1

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.

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.5B_0^2/n_0(1+T_e/T_i)\approx 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\approx 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.

FIGURE 2

FIGURE 2. Density, pressure, and temperature of initially hot and cold ions at tω_ci = 120: (A)n_ih, (B)n_ic, (C)P_ih, (D)P_ic(E)T_ih, (F)T_ic. 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.

4 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 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_ydv_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.

FIGURE 3

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).

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 $-\widehat{y}$ by the Lorentz force v_zB_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.

4.1 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_{{\varphi }_z}=\sqrt{2e{\varphi }_z/m_i}$ corresponding to the Hall electric field potential ϕ_z = −∫E_zd_z. We have grounded the potential to a region just outside the separatrix. The increase in $v_{{\varphi }_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_{{\varphi }_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_{{\varphi }_z}$ and v_E×B,z at x = 85.

FIGURE 4

FIGURE 4. Reduced ion distributions f(z, v_z) at (A)x = 75 and (B) 85 at tω_ci = 120. The solid, dashed, and dotted lines show v_E×B,z, $v_{{\varphi }_z}=\sqrt{2e{\varphi }_z/m_i}$, where ϕ_z = −∫E_zdz, and the location of the separatrix, respectively. Further away from the X line, at x = 75, the cold inbound ions remain cold and follow v_E×B,z deep within the exhaust. Closer to the X line, at x = 85, the cold inbound ions begin to thermalize and deviate from v_E×B,z shortly after crossing the separatrix. Although the phase space density of the inbound beam decreases away from the separatrix, the remnants of it follow $v_{{\varphi }_z}$. The thermal population consists largely of ions that have already crossed z = 0 one or multiple times.

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 integrated in the dynamically changing electric and magnetic fields. The initial velocities of the test ions are chosen from the cold component of $f_{ic}^{top}(x,z=0,t=120)$, see Figures 5A,B. They are integrated forward in time until t = 125, and backward until t = 90. The location of the ions are shown overlaid E_z (Figures 5C–G) and $n_{ic}^{top}$ (Figures 5H–L) for t = [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⁻³, 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 dawn-dusk 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.

FIGURE 5

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_{ic}^{top}$, respectively. The subdivision into yellow/black/white/green populations are partly based on the initial conditions, and further explained in the text.

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_{ic}^{top}$.

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_{ic}^{top}$ 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.

4.2 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 ${\mathbf{k}}_B=\widehat{b}\cdot \nabla \widehat{b}$, where $\widehat{b}=\mathbf{B}/|\mathbf{B}|$ is the magnetic field unit vector. Büchner and Zelenyi (Büchner and Zelenyi, 1989) defined the curvature parameter κ² = r_B/ρ_i, where r_B = |k_B|⁻¹ 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 pitch-angle scattering during the neutral sheet crossing. Figure 6B shows the 2D map of log ₁₀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 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,⊥ = κ⁻²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.

FIGURE 6

FIGURE 6. (A) Density of initially cold ions originating from the north, $n_{ic}^{top}$. (B) Magnetic field curvature radius   log ₁₀r_B. (C) Magnetic field magnitude |B|, r_B, and $v_{i,\perp }^{\kappa =1}=r_B{\omega }_{ci}$. (D–H) Reduced distributions $f_{ic}^{top}(v_x,v_z)$ at z = 0, and six different x, marked by white squares in (A, B). The cross ’×’ shows the bulk velocity, the dot ’⋅’ shows v_E×B, and the larger circle shows $v_{i,\perp }^{\kappa =1}$, centered on v_E×B. (I–N) Reduced distributions f(s, v) along the two different field lines marked as white lines in (A, B). Panels (I–K) corresponds to the field lines located furthest downstream. The black line shows v_E×B. The s is the length along the field line, where s = 0 corresponds to z = 0, and s > 0 (s < 0) corresponds to z < 0 (z > 0).

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,\perp }=\sqrt{v_x^2+v_y^2}$. The circles overlaid the distributions show v_i,⊥ = κ⁻²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_{i,\perp }^{\kappa =1}$ (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_{i,\perp }^{\kappa =1}$ is relatively large, and even encompass the majority of the thermal population. At x = 74, the cold population is still well inside $v_{i,\perp }^{\kappa =1}$, 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_{i,\perp }^{\kappa =1}$. 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_{ic}^{top}$ along two select field lines (marked as white lines in Figures 6A,B). In Figures 6I–K, we show the highlighted field line closest to the front, in the region where $n_{ic}^{top}$ extends below z = −2 and $v_{i,\perp }^{\kappa =1}\gtrsim 2$ at z = 0. In Figures 6L–N, we show the highlighted field line in the intermediate exhaust, where $n_{ic}^{top}$ has a cut-off at z = −2 and $v_{i,\perp }^{\kappa =1}\approx 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_{ic}^{top}(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_{ic}^{top}(v_x)$ and $f_{ic}^{top}(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 from $\widehat{x}$ to $\widehat{z}$ 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_{ic}^{top}(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_{ic}^{top}(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_xB_y play minor roles in the acceleration along $\widehat{z}$, the main responsible force is − v_yB_x. The ions that are reflected have a speed v_y > 0 and therefore experience an upward force − v_yB_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_yB_x < 0. We also note that we can still identify the cold component in $f_{ic}^{top}(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_{ic}^{top}(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_yB_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.

4.3 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_yB_x, and thereafter towards − x by v_yB_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.

FIGURE 7

FIGURE 7. Electrostatic solitary waves inside the southern exhaust. The arrow recurring in every panel identifies one given structure.

FIGURE 8

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). (A) Parallel electric field E_‖. (B) Reduced ion distribution of initially cold ions, log ₁₀f_ic. (C) Mixing level of ions originating from the north and south, respectively, $f_{mix}=f_{ic}^{top}/f_{ic}$. For f_mix = 1 (red) and 0 (purple), all the ions are originating from the top and bottom, respectively. For f_mix = 0.5 (yellow), half of the ions are originating from the top, and half from bottom. For reference we show log ₁₀f_ic = −1 as black contours. (D) Reduced distribution of originally cold electrons, log ₁₀f_ec. We have saturated the color scale to highlight the features in the region of the waves. (E) Parallel electric field E_‖. (F) Ion, and (G) electron distributions at the location marked by a black square in panel (E), and the vertical dashed lines in (A–D). The dashed green line marks a fit to the data using three Maxwellians. (H) Dispersion relation obtained based on the Maxwellian fits in (F–G). The positive growth rate is due to an ion-acoustic instability. The phase speed v_ph and wavelength λ at max growth γ_max is shown as a black bar in (B).

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}=f_{ic}^{top}/f_{ic}$ is shown in 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 ₁₀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:

$0=\sum _s\frac{{\omega }_{ps}^2}{k^2{v}_{ts}^2}{Z}^{\prime }\left(\frac{\omega -k{v}_{ds}}{k{v}_{ts}}\right)$

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}=\sqrt{2T_s/m_s}$ is the parallel thermal speed, of population s, ω = ω_r + iγ is the complex frequency with real frequency ω, and imaginary frequency γ 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 (ions) and Figure 8G (electrons): the initially hot ions and electrons (Harris sheet populations), the initially cold ions and electrons from the top and bottom, respectively. The parameters are: n_i = [0.05, 0.08, 0.15], n_e = [0.11, 0.13, 0.04], v_di = [0.5, 1.15, 0.05] v_de = [−2.1, 2.7, 0.5], v_ti = [1.5, 0.3, 0.07], v_te = [2.8, 3.0, 4.0]. The initially hot plasma does not change the results, but we include them for completeness. Effectively, there are one electron and two ion distributions: (1) a mix of all electrons combined, (2) the cold inbound ion population at v_‖ ≈ 0 consisting of ions originating from the southern inflow, and (3) the outbound ion beam at v_‖ > 0 consisting of a mixture of ions originating from the top and bottom. The resulting dispersion relation is shown in Figure 8H. Positive growth rates are found for 0 < kd_i ≲ 60, with a peak growth rate γ_max ≲ 0.5 at k_maxd_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.

5 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_zB_x < 0 acting on the ions during their first inbound leg towards the neutral plane from the inflow regions. Where v_y < 0, v_zB_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_iv_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_yB_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_zB_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.

FIGURE 9

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_iv_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.

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}$, $\mathrm{\Delta }{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

$\mathrm{\Delta }{v}_y=-\frac{q}{m_i}\mathrm{\Delta }{A}_y,$

where $\mathrm{\Delta }{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)\gg 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 $-\widehat{y}$ by the force v_zB_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_yB_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 = −∂_zA_y and B_z = ∂_xA_y. The information of both magnetic field components in the force terms v_zB_x and v_yB_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.

6 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 (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 $\widehat{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.

FIGURE 10

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.

Data Availability Statement

The dataset used in this study can be found at UiB Open Research data: https://doi.org/10.18710/0257WR.

Author Contributions

CN planned the study and executed the simulation, performed the data analysis, scientific interpretation and lead the manuscript writing. PT, MH, YX shared in the conception of the study. PT, MH shared in the simulation execution and data analysis. PT, MH, ST, WL, YX, NK, SS, HK shared in the scientific interpretation and discussion of the simulation results. All authors shared in the writing of the manuscript.

Funding

This study was supported by NOTUR/NORSTOR under project NN9496K. CN and PT received support from the Research Council of Norway under contract 300865. MH's work prior to moving to NASA ARC was supported by NASA contract NNG04EB99C at SwRI. STR, WL, and PT acknowledges support from the International Space Science Institute (ISSI) international team ‘Cold plasma of ionospheric origin in the Earth's magnetosphere’. The open access publishing fee is funded by the University of Bergen.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

Alm, L., André, M., Vaivads, A., Khotyaintsev, Y. V., Torbert, R. B., Burch, J. L., et al. (2018). Magnetotail Hall Physics in the Presence of Cold Ions. Geophys. Res. Lett. 45, 10941–10950. doi:10.1029/2018GL079857

CrossRef Full Text | Google Scholar

André, M., Li, W., Toledo-Redondo, S., Khotyaintsev, Y. V., Vaivads, A., Graham, D. B., et al. (2016). Magnetic Reconnection and Modification of the Hall Physics Due to Cold Ions at the Magnetopause. Geophys. Res. Lett. 43, 6705–6712. doi:10.1002/2016GL069665

CrossRef Full Text | Google Scholar

Artemyev, H., Angelopoulos, A. V., and au, V. (2017). Ion Dynamics in Magnetotail Reconnection in the Presence of Density Asymmetry. J. Geophys. Res. Space Phys. 122, 2010–2023. doi:10.1002/2016JA023651

CrossRef Full Text | Google Scholar

Arzner, K., and Scholer, M. (2001). Kinetic Structure of the post Plasmoid Plasma Sheet during Magnetotail Reconnection. J. Geophys. Res. 106, 3827–3844. doi:10.1029/2000JA000179

CrossRef Full Text | Google Scholar

Baumjohann, W., Paschmann, G., and Cattell, C. A. (1989). Average Plasma Properties in the central Plasma Sheet. J. Geophys. Res. 94, 6597–6606. doi:10.1029/ja094ia06p06597

CrossRef Full Text | Google Scholar

Birn, J., Nakamura, R., Panov, E. V., and Hesse, M. (2011). Bursty Bulk Flows and Dipolarization in MHD Simulations of Magnetotail Reconnection. J. Geophys. Res. 116, A01210. doi:10.1029/2010JA016083

CrossRef Full Text | Google Scholar

Birn, J., Runov, A., and Zhou, X. Z. (2017). Ion Velocity Distributions in Dipolarization Events: Distributions in the central Plasma Sheet. J. Geophys. Res. Space Phys. 122, 8014–8025. doi:10.1002/2017JA024230

CrossRef Full Text | Google Scholar

Borovsky, J. E., and Denton, M. H. (2008). A Statistical Look at Plasmaspheric Drainage Plumes. J. Geophys. Res. (Space Physics) 113, A09221. doi:10.1029/2007ja012994

CrossRef Full Text | Google Scholar

Büchner, J., and Zelenyi, L. M. (1989). Regular and Chaotic Charged Particle Motion in Magnetotaillike Field Reversals: 1. Basic Theory of Trapped Motion. J. Geophys. Res. 94, 11821–11842. doi:10.1029/JA094iA09p11821

CrossRef Full Text | Google Scholar

Dargent, J., Aunai, N., Lavraud, B., Toledo‐Redondo, S., and Califano, F. (2020). Simulation of Plasmaspheric Plume Impact on Dayside Magnetic Reconnection. Geophys. Res. Lett. 47, e86546. doi:10.1029/2019GL086546

CrossRef Full Text | Google Scholar

Divin, A., Khotyaintsev, Y. V., Vaivads, A., André, M., Toledo-Redondo, S., Markidis, S., et al. (2016). Three-scale Structure of Diffusion Region in the Presence of Cold Ions. J. Geophys. Res. Space Phys. 121, 001–012. doi:10.1002/2016JA023606

CrossRef Full Text | Google Scholar

Drake, J. F., Swisdak, M., Phan, T. D., Cassak, P. A., Shay, M. A., Lepri, S. T., et al. (2009). Ion Heating Resulting from Pickup in Magnetic Reconnection Exhausts. J. Geophys. Res. (Space Physics) 114, A05111. doi:10.1029/2008ja013701

CrossRef Full Text | Google Scholar

Drake, H., Phan, J. F., Eastwood, T. D., McFadden, J. P., and au, J. P. (2015). Ion Temperature Anisotropy across a Magnetotail Reconnection Jet. Geophys. Res. Lett. 42, 7239–7247. doi:10.1002/2015GL065168

PubMed Abstract | CrossRef Full Text | Google Scholar

Eastman, T. E., DeCoster, R. J., and Frank, L. A. (1986). Velocity Distributions of Ion Beams in the Plasma Sheet Boundary Layer. Wash. DC Am. Geophys. Union Geophys. Monogr. Ser. 38, 117–126. doi:10.1029/GM038p0117

CrossRef Full Text | Google Scholar

Eastwood, J. P., Goldman, M. V., Newman, H., Mistry, D. L., Lapenta, R., and au, G. (2015). Ion Reflection and Acceleration Near Magnetotail Dipolarization Fronts Associated with Magnetic Reconnection. J. Geophys. Res. Space Phys. 120, 511–525. doi:10.1002/2014JA020516

CrossRef Full Text | Google Scholar

Engwall, E., Eriksson, A. I., Cully, C. M., André, M., Torbert, R., and Vaith, H. (2009). Earth's Ionospheric Outflow Dominated by Hidden Cold Plasma. Nat. Geosci 2, 24–27. doi:10.1038/ngeo387

CrossRef Full Text | Google Scholar

Ergun, R. E., Goodrich, K. A., Wilder, F. D., Ahmadi, N., Holmes, J. C., Eriksson, S., et al. (2018). Magnetic Reconnection, Turbulence, and Particle Acceleration: Observations in the Earth's Magnetotail. Geophys. Res. Lett. 45, 3338–3347. doi:10.1002/2018GL076993

CrossRef Full Text | Google Scholar

Fried, B. D., and Conte, S. D. (1961). The Plasma Dispersion Function. New York, NY: Academic Press.

Fu, H. S., Cao, J. B., Khotyaintsev, Y. V., Sitnov, M. I., Runov, A., Fu, S. Y., et al. (2013). Dipolarization Fronts as a Consequence of Transient Reconnection: In Situ Evidence. Geophys. Res. Lett. 40, 6023–6027. doi:10.1002/2013GL058620

CrossRef Full Text | Google Scholar

Fuselier, S. A., Burch, J. L., Mukherjee, J., Genestreti, K. J., Vines, S. K., Gomez, R., et al. (2017). Magnetospheric Ion Influence at the Dayside Magnetopause. J. Geophys. Res. Space Phys. 122, 8617–8631. doi:10.1002/2017ja024515

CrossRef Full Text | Google Scholar

Glocer, A., Welling, D., Chappell, C. R., Toth, G., Fok, M. C., Komar, C., et al. (2020). A Case Study on the Origin of Near‐Earth Plasma. J. Geophys. Res. Space Phys. 125, e28205. doi:10.1029/2020JA028205

CrossRef Full Text | Google Scholar

Gosling, J. T., Thomsen, M. F., Bame, S. J., Onsager, T. G., and Russell, C. T. (1990). The Electron Edge of Low Latitude Boundary Layer during Accelerated Flow Events. Geophys. Res. Lett. 17, 1833–1836. doi:10.1029/GL017i011p01833

CrossRef Full Text | Google Scholar

Graham, D. B., Khotyaintsev, Y. V., Norgren, C., Vaivads, A., André, M., Toledo‐Redondo, S., et al. (2017). Lower Hybrid Waves in the Ion Diffusion and Magnetospheric Inflow Regions. J. Geophys. Res. Space Phys. 122, 517–533. doi:10.1002/2016JA023572

CrossRef Full Text | Google Scholar

Haggerty, C. C., Shay, M. A., Drake, J. F., Phan, T. D., and McHugh, C. T. (2015). The Competition of Electron and Ion Heating during Magnetic Reconnection. Geophys. Res. Lett. 42, 9657–9665. doi:10.1002/2015GL065961

CrossRef Full Text | Google Scholar

Hesse, M., Schindler, K., Birn, J., and Kuznetsova, M. (1999). The Diffusion Region in Collisionless Magnetic Reconnection. Phys. Plasmas 6, 1781–1795. doi:10.1063/1.873436

CrossRef Full Text | Google Scholar

Lapenta, G., Markidis, S., Divin, A., Goldman, M., and Newman, D. (2010). Scales of Guide Field Reconnection at the Hydrogen Mass Ratio. Phys. Plasmas 17, 082106. doi:10.1063/1.3467503

CrossRef Full Text | Google Scholar

Li, S.-S., Liu, J., Angelopoulos, V., Runov, A., Zhou, X.-Z., and Kiehas, S. A. (2014). Antidipolarization Fronts Observed by Artemis. J. Geophys. Res. Space Phys. 119, 7181–7198. doi:10.1002/2014JA020062

CrossRef Full Text | Google Scholar

Li, W. Y., André, M., Khotyaintsev, Y. V., Vaivads, A., Fuselier, S. A., Graham, D. B., et al. (2017). Cold Ionospheric Ions in the Magnetic Reconnection Outflow Region. J. Geophys. Res. Space Phys. 122, 10,194–10,202. doi:10.1002/2017JA024287

CrossRef Full Text | Google Scholar

Lindstedt, T., Khotyaintsev, Y. V., Vaivads, A., André, M., Fear, R. C., Lavraud, B., et al. (2009). Separatrix Regions of Magnetic Reconnection at the Magnetopause. Ann. Geophys. 27, 4039–4056. doi:10.5194/angeo-27-4039-2009

CrossRef Full Text | Google Scholar

Liu, J., Angelopoulos, V., Runov, A., and Zhou, X. Z. (2013). On the Current Sheets Surrounding Dipolarizing Flux Bundles in the Magnetotail: The Case for Wedgelets. J. Geophys. Res. (Space Physics) 118, 2000–2020. doi:10.1002/jgra.50092

CrossRef Full Text | Google Scholar

Liu, C. M., Vaivads, A., Graham, D. B., Khotyaintsev, Y. V., Fu, H. S., Johlander, A., et al. (2019). Ion-Beam-Driven Intense Electrostatic Solitary Waves in Reconnection Jet. Geophys. Res. Lett. 46, 12702–12710. doi:10.1029/2019gl085419

CrossRef Full Text | Google Scholar

Nagai, T., Fujimoto, M., Saito, Y., Machida, S., Terasawa, T., Nakamura, R., et al. (1998). Structure and Dynamics of Magnetic Reconnection for Substorm Onsets with Geotail Observations. J. Geophys. Res. 103, 4419–4440. doi:10.1029/97JA02190

CrossRef Full Text | Google Scholar

Nagai, T., Nakamura, M., Shinohara, I., Fujimoto, M., Saito, Y., and Mukai, T. (2002). Counterstreaming Ions as Evidence of Magnetic Reconnection in the Recovery Phase of Substorms at the Kinetic Level. Phys. Plasmas 9, 3705–3711. doi:10.1063/1.1499117

CrossRef Full Text | Google Scholar

Nagai, T., Shinohara, I., and Zenitani, S. (2015). Ion Acceleration Processes in Magnetic Reconnection: Geotail Observations in the Magnetotail. J. Geophys. Res. Space Phys. 120, 1766–1783. doi:10.1002/2014JA020737

CrossRef Full Text | Google Scholar

Nakamura, M. S., Fujimoto, M., and Maezawa, K. (1998). Ion Dynamics and Resultant Velocity Space Distributions in the Course of Magnetotail Reconnection. J. Geophys. Res. 103, 4531–4546. doi:10.1029/97JA01843

CrossRef Full Text | Google Scholar

Nakamura, R., Baumjohann, W., Mouikis, C., Kistler, L. M., Runov, A., Volwerk, M., et al. (2004). Spatial Scale of High-Speed Flows in the Plasma Sheet Observed by Cluster. Geophys. Res. Lett. 31, L09804. doi:10.1029/2004gl019558

CrossRef Full Text | Google Scholar

Norgren, C., Hesse, M., Graham, D. B., Khotyaintsev, Y. V., Tenfjord, P., Vaivads, A., et al. (2020). Electron Acceleration and Thermalization at Magnetotail Separatrices. J. Geophys. Res. (Space Physics) 125, e27440. doi:10.1029/2019ja027440

CrossRef Full Text | Google Scholar

Northrop, T. G. (1963). Adiabatic Charged-Particle Motion. Rev. Geophys. 1, 283–304. doi:10.1029/RG001i003p00283

CrossRef Full Text | Google Scholar

Runov, A., Angelopoulos, V., Zhou, X.-Z., Zhang, X.-J., Li, S., Plaschke, F., et al. (2011). A THEMIS Multicase Study of Dipolarization Fronts in the Magnetotail Plasma Sheet. J. Geophys. Res. 116, A05216. doi:10.1029/2010JA016316

CrossRef Full Text | Google Scholar

Schriver, D., and Ashour-Abdalla, M. (1990). Cold Plasma Heating in the Plasma Sheet Boundary Layer: Theory and Simulations. J. Geophys. Res. 95, 3987–4005. doi:10.1029/JA095iA04p03987

CrossRef Full Text | Google Scholar

Sitnov, M. I., Swisdak, M., and Divin, A. V. (2009). Dipolarization Fronts as a Signature of Transient Reconnection in the Magnetotail. J. Geophys. Res. 114, A04202. doi:10.1029/2008JA013980

CrossRef Full Text | Google Scholar

Sitnov, M. I., Buzulukova, N., Swisdak, M., Merkin, V. G., and Moore, T. E. (2013). Spontaneous Formation of Dipolarization Fronts and Reconnection Onset in the Magnetotail. Geophys. Res. Lett. 40, 22–27. doi:10.1029/2012GL054701

CrossRef Full Text | Google Scholar

Speiser, T. W. (1965). Particle Trajectories in Model Current Sheets: 1. Analytical Solutions. J. Geophys. Res. 70, 4219–4226. doi:10.1029/jz070i017p04219

CrossRef Full Text | Google Scholar

Steinvall, K., Khotyaintsev, Y. V., Graham, D. B., Vaivads, A., André, M., and Russell, C. T. (2021). Large Amplitude Electrostatic Proton Plasma Frequency Waves in the Magnetospheric Separatrix and Outflow Regions during Magnetic Reconnection. Geophys. Res. Lett. 48, e2020GL090286. doi:10.1029/2020gl090286

CrossRef Full Text | Google Scholar

Tenfjord, P., Hesse, M., and Norgren, C. (2018). The Formation of an Oxygen Wave by Magnetic Reconnection. J. Geophys. Res. Space Phys. 123, 9370–9380. doi:10.1029/2018JA026026

CrossRef Full Text | Google Scholar

Tenfjord, P., Hesse, M., Norgren, C., Spinnangr, S. F., Kolstø, H., and Kwagala, N. (2020). Interaction of Cold Streaming Protons with the Reconnection Process. J. Geophys. Res. Space Phys. 125, e27619. doi:10.1029/2019JA027619

CrossRef Full Text | Google Scholar

Toledo-Redondo, S., Vaivads, A., André, M., and Khotyaintsev, Y. V. (2015). Modification of the Hall Physics in Magnetic Reconnection Due to Cold Ions at the Earth's Magnetopause. Geophys. Res. Lett. 42, 6146–6154. doi:10.1002/2015gl065129

CrossRef Full Text | Google Scholar

Toledo-Redondo, S., André, M., Khotyaintsev, Y. V., Vaivads, A., Walsh, A., Li, W., et al. (2016a). Cold Ion Demagnetization Near the X-Line of Magnetic Reconnection. Geophys. Res. Lett. 43, 6759–6767. doi:10.1002/2016GL069877

CrossRef Full Text | Google Scholar

Toledo-Redondo, S., André, M., Vaivads, A., Khotyaintsev, Y. V., Lavraud, B., Graham, D. B., et al. (2016b). Cold Ion Heating at the Dayside Magnetopause during Magnetic Reconnection. Geophys. Res. Lett. 43, 58–66. doi:10.1002/2016gl069877

CrossRef Full Text | Google Scholar

Toledo-Redondo, S., André, M., Khotyaintsev, Y. V., Lavraud, B., Vaivads, A., Graham, D. B., et al. (2017). Energy Budget and Mechanisms of Cold Ion Heating in Asymmetric Magnetic Reconnection. J. Geophys. Res. Space Phys. 122, 9396–9413. doi:10.1002/2017JA024553

CrossRef Full Text | Google Scholar

Toledo‐Redondo, S., Dargent, J., Aunai, N., Lavraud, B., André, M., Li, W., et al. (2018). Perpendicular Current Reduction Caused by Cold Ions of Ionospheric Origin in Magnetic Reconnection at the Magnetopause: Particle‐in‐Cell Simulations and Spacecraft Observations. Geophys. Res. Lett. 45, 10,033–10,042. doi:10.1029/2018GL079051

CrossRef Full Text | Google Scholar

Toledo‐Redondo, S., André, M., Aunai, N., Chappell, C. R., Dargent, J., Fuselier, S. A., et al. (2021). Impacts of Ionospheric Ions on Magnetic Reconnection and Earth's Magnetosphere Dynamics. Rev. Geophys. 59, e2020RG000707. doi:10.1029/2020RG000707

CrossRef Full Text | Google Scholar

Vines, S. K., Fuselier, S. A., Trattner, K. J., Burch, J. L., Allen, R. C., Petrinec, S. M., et al. (2017). Magnetospheric Ion Evolution across the Low‐Latitude Boundary Layer Separatrix. J. Geophys. Res. Space Phys. 122, 10,247–10,262. doi:10.1002/2017JA024061

CrossRef Full Text | Google Scholar

Wygant, J. R., Cattell, C. A., Lysak, R., Song, Y., Dombeck, J., McFadden, J., et al. (2005). Cluster Observations of an Intense normal Component of the Electric Field at a Thin Reconnecting Current Sheet in the Tail and its Role in the Shock-like Acceleration of the Ion Fluid into the Separatrix Region. J. Geophys. Res. 110, A09206. doi:10.1029/2004JA010708

CrossRef Full Text | Google Scholar

Xu, Y., Fu, H. S., Norgren, C., Toledo‐Redondo, S., Liu, C. M., and Dong, X. C. (2019). Ionospheric Cold Ions Detected by MMS behind Dipolarization Fronts. Geophys. Res. Lett. 46, 7883–7892. doi:10.1029/2019GL083885

CrossRef Full Text | Google Scholar

Zenitani, S., Shinohara, I., Nagai, T., and Wada, T. (2013). Kinetic Aspects of the Ion Current Layer in a Reconnection Outflow Exhaust. Phys. Plasmas 20, 092120. doi:10.1063/1.4821963

CrossRef Full Text | Google Scholar