The Dynamic Evolution of Solar Wind Streams Following Interchange Reconnection

Interchange reconnection is thought to play an important role in determining the dynamics and material composition of the slow solar wind that originates from near coronal hole boundaries. To explore the implications of this process we simulate the dynamic evolution of a solar wind stream along a newly-opened magnetic flux tube. The initial condition is composed of a piecewise continuous dynamic equilibrium in which the regions above and below the reconnection site are extracted from steady-state solutions along open and closed field lines. The initial discontinuity at the reconnection site is highly unstable and evolves as a Riemann problem, decomposing into an outward-propagating shock and inward-propagating rarefaction that eventually develop into a classic N-wave configuration. This configuration ultimately propagates into the heliosphere as a coherent structure and the entire system eventually settles to a quasi-steady wind solution. In addition to simulating the fluid evolution we also calculate the time-dependent non-equilibrium ionization of oxygen in real time in order to construct in situ diagnostics of the conditions near the reconnection site. This idealized description of the plasma dynamics along a newly-opened magnetic field line provides a baseline for predicting and interpreting the implications of interchange reconnection for the slow solar wind. Notably, the density and velocity within the expanding N-wave are generally enhanced over the ambient wind, as is the O7+/O6+ ionization ratio, which exhibits a discontinuity across the reconnection site that is transported by the flow and arrives later than the propagating N-wave.


INTRODUCTION
One of the more enduring challenges in solar and heliospheric science is to determine what physical processes are responsible for the disparate properties of various solar wind streams. Compared to the fast solar wind (FSW), the slow solar wind (SSW) exhibits lower typical wind speeds with increased variability (McComas et al. 2000). Various authors have explored the effects of different types of heating profiles and magnetic geometries in field-aligned (1D) models (e.g., Cranmer et al. 2007;Grappin et al. 2011, and others), which can produce a wide range of predicted wind speeds but do not predict rapid fluctuations within a single wind stream. Additionally, the material composition of the SSW shows similarities to the plasma found in active regions, with ionization ratios and elemental abundances that are inconsistent with the coronal-hole regions from which the FSW emanates (Zhao et al. 2017;Laming et al. 2019). This suggests that changes in magnetic topology, which are difficult to capture in field-aligned models, may be critical to the properties of the SSW.
The solar corona is partitioned into various spatial domains that reflect the connectivity of the magnetic field. In regions where the magnetic field maps from the photosphere into the heliosphere and beyond, the associated volume is said to be "open". Conversely, where the field maps between positive and negative polarity domains on the photosphere, the associated volume is said to be "closed". The interfaces between these domains collectively form the open-closed boundary, which is strongly correlated with coronal-hole boundaries in X-ray and EUV images (e.g., Nikolić 2019). According to Alfvén's theorem, an ideal plasma (zero electrical resistivity) is coupled to the magnetic field in such a way that individual fluid elements are always connected to the same field line. Therefore, when the connectivity of the magnetic field is static in time the plasma within the various closed magnetic domains remains confined while the plasma within open domains can stream freely into the heliosphere.
For plasma that originates from within a closed domain to escape across the open-closed boundary into an adjacent coronal hole depends on the process of interchange reconnection (Crooker et al. 2002), whereby field lines from either side of the open-closed boundary undergo a change in connectivity so that a new flux tube is formed that extends from a footpoint within a previously closed region on the solar surface out into the open coronal hole region and into the heliosphere and beyond. This process has been studied extensively in 3D magnetohydrodynamic (MHD) simulations of coronal streamers and pseudostreamers (e.g., Masson et al. 2014;Higginson et al. 2017;Aslanyan et al. 2021;Scott et al. 2021); however, these and similar studies have invoked simplistic fluid approximations for computational efficiency. Therefore, while the magnetic evolution is reasonably well understood, many questions remain concerning the associated plasma dynamics. Additionally, while the ionization ratios of trace elements have been studied in detail for steady-state wind streams (see, e.g., Landi et al. 2012;Gilly & Cranmer 2020, and others) little is known about how these evolve during interchange reconnection. These ratios are often taken as diagnostics of the temperature history from the solar surface up to a given height (the so-called "freeze-in height", Owocki et al. 1983) beyond which the ionization ratios undergo no additional evolution; however, if plasma is released from closed magnetic domains at heights that are comparable to the freeze-in height this could have significant implications for in-situ diagnostics of the solar wind.
The primary difficulty in constructing an accurate model of interchange reconnection in the context of the solar wind stems from the many decades of spatial and temporal scales that must be simultaneously resolved. At the largest scale the global magnetic field and location of the acoustic and Alfvénic critical points require a numerical domain that covers several tens of Mm 2 on the solar surface (i.e., the size of an active region) and extends outward to a height of 20R or more, meaning that the time required for fluid features to transit the numerical domain is typically on the order of tens of hours for the fastest propagating signals (e.g., electron thermal speed and fast magneto-acoustic mode) and significantly longer for slower modes (e.g., slow magnetoacoustic mode and bulk flow). At the opposite extreme, the temperature gradients at the base of the transition region -which are required to accurately capture the subtle balance of coronal heating, thermal conduction, and radiative losses that dictate the mass flux into the corona -require a minimum resolution on the order of a few tens of km, resulting in a numerical timestep that is typically on the order of 10 −2 s. As a result, constructing a 3D numerical model of interchange reconnection in the presence of a self-consistent solar wind has proven to be computationally prohibitive, and none are currently in common use.
Here we present a hybrid model that partially mitigates the computational demands of a full 3D simulation by computing the field-aligned (1D) plasma dynamicsincluding detailed physical processes that structure the transition region -with an empirical mechanism for emulating the effects of magnetic reconnection. Beginning with steady-state solutions for the open-and closed-field regions, we construct an initial condition for the plasma along a newly-opened flux tube such that the fluid properties are discontinuous across the reconnection site, being composed of a transonic wind stream above a hotter and more dense hydrostatic column. This technique builds on previous work by Bradshaw et al. (2011), but has been extended to include critical physics for the appropriate treatment of an out-flowing wind solution, including modification of the distant outer boundary and variability of the flux-tube cross-sectional area. By following the hydrodynamic evolution of this newlyformed, post-reconnection plasma column, we are able explore how magnetic reconnection directly affects the plasma evolution when material from disparate source regions are brought into close proximity as the magnetic connectivity changes.
In the following section we describe our simulation design and numerical model. Then in section 3 we describe the evolution of the post-reconnection plasma and associated time-dependent ionization of oxygen, which we calculate for comparison to in situ studies by Zhao et al. (2017). In section 4 we discuss these findings and their implications for in situ measurements as well as the appropriate interpretation of this model in the context of a volume filling, 3D magnetic field. Finally, we offer concluding remarks in section 5.

Fluid Model
We model the plasma as a fully-ionized twotemperature fluid with mass density ρ, velocity u, and individual ion and electron pressures p i and p e , and we solve the field-aligned hydrodynamic equations in time (t) and distance along the field line (s). The magnetic field geometry is taken to be static in time so that the cross-field dynamics can be ignored. In the limit that the electron mass is negligible compared to the average ion mass (m i m e ), and assuming quasi-neutrality (n e = n i = n), the continuity and momentum equations are and where g and A are the gravitational acceleration and cross-sectional area parameterized along s, P = p i + p e is the total pressure, and σ is the viscous stress. The conservative form of the energy equations for the ions and electrons are likewise given by and where the energies are defined as E i = 1 2 ρu 2 + 1 γ−1 p i for the ions and E e = 1 γ−1 p e for the electrons; T i and T e are the ion and electron temperatures; F i and F e are the associated conductive heat fluxes; and h i,e and R are the imposed coronal heating and empirical radiative losses, respectively.
We model the radiative losses with a piecewise polynomial fit as described in Bradshaw et al. (2011) and we prescribe these to fall off to zero below a minimum temperature of 2×10 4 K, emulating the effect of an optically thick chromosphere. The system is closed through the ideal gas law, p i,e = k B nT i,e , and by the definition of the conductive fluxes and ion viscous stress through the Spitzer-Harm formulation; F i,e = −κ i,e T

5/2
i,e ∂ s T i,e and σ = (4/3)µ i ∂ s u. The dynamic viscosity is given by c and ion-electron collision frequency is defined to be ν ei = 4.82 n ln Λ ie c m e /m i . The Coulomb logarithms for ion-ion and ion-electron collisions (ln Λ ii c and ln Λ ie c ) are defined as in the NRL Formulary (Huba 1998) and discussed by Fitzpatrick (2015).
Specific values for the physical constants are listed in Table 1. These are broadly in line with accepted solar values, with the exception of the the coefficient of dynamic viscosity, which is further modified by a multiplicative factor of 1/(1 + α 2 ), where α = 10 2 × λ mfp i ∂ s ln A. This multiplicative factor guarantees that the dynamic viscosity becomes small as the ion meanfree-path (λ mfp i = 8.4 × 10 3 T 2 i /n) approaches or exceeds the geometric length scale, emulating the transition to a collisionless regime as discussed by Endeve & Leer (2001)  In addition to the above fluid equations we also calculate the time-dependent ionization of oxygen as a trace element subject to the evolving fluid. This evolution is decoupled from the energy balance, and does not affect the radiative loss function, which would require tracking additional elements such as He, Fe, S, C, etc., in order to be calculated self-consistently. The ion charge-state distribution of oxygen, with atomic number Z = 8, evolves according to where D t ≡ (∂ t + u ∂ s ) is the material derivative along the flow and Y j are the individual ionization fractions (i.e., the fraction of a given species with atomic number Z having a specific ionization level j). I j and R j are the temperature-dependent ionization and recombination rates from state j, which are calculated as in Bradshaw & Mason (2003) with updated values from version 8 of the CHIANTI atomic database (Dere et al. 1997;Del Zanna et al. 2015). Prohibiting ionization and recombination beyond the range of allowed states, the sum of Eq. (5) over all ionization levels vanishes by construction and the system is both normalized and regularized so that

Steady-state Conditions
We use the HYDRAD code (Bradshaw & Cargill 2013) to solve equations (1)-(5) on a numerical grid of 512 cells that are exponentially spaced in heliocentric radius, spanning a range of 0 ≤ s ≤ 30R above the solar surface, where R = 700 Mm is the radius of the sun. Beyond this base grid we allow up to 16 levels of adaptive refinement through pairwise splitting and merging of cells. The grid spacing ∆s is checked against the length scale of the primitive variables (λ HD ) after every tenth timestep in order to maintain the resolution requirement 0.05 ≤ ∆s/λ HD ≤ 0.1. During adaptive refinement the mass, momentum, and energy densities are explicitly conserved to numerical accuracy. At the limit of refinement the minimum allowable grid spac-ing is ∆s min ≈ 10 4 cm at the base of the numerical domain; however, this extreme resolution is never required in practice, indicating that the transition region is fully resolved at all times.
In order to determine steady-state conditions for the open-and closed-field regions the calculation is first initialized with temperature and density profiles that are generated from a separate routine assuming zero velocity and uniform heating. The magnetic geometry is spherically symmetric with and where g = −2.74 × 10 4 cm s −2 is the gravitational acceleration at the solar surface. The cross-sectional area is dimensionless and normalized to unity at the solar surface. The density and temperature are set to T 0 = 2 × 10 4 K and n 0 = 10 10 cm −3 at a height of 10 9 cm above the solar surface, which sets the location of the base of the transition region within the domain. Above this height the temperatures of the ions and electrons rise abruptly through the transition region to coronal values, while below this level the temperature is effectively uniform and the density increases exponentially toward the interior of the solar surface. At the beginning of the simulation the initiallyuniform heating is replaced with a superposition of three exponentials of the form h i (s) = h 0 exp(−s/s l ), after which the system is advanced in time and allowed to relax toward a new steady-state solution. Exponential heating functions have been used extensively in solar wind models (see, e.g., Pinto et al. 2009) and are employed here for simplicity; although in the future we intend to incorporate more realistic mechanisms (such as those discussed by Cranmer et al. 2007). Specific values for the three ion heating rates and scale heights are given in Table 2. The electrons are not heated directly (h e = 0), but instead rely on collisional coupling in the lower corona and thermal conduction in the extended corona to maintain their temperature. This particular heating model serves to give realistic values for the electron temperature in the lower-mid corona and for the mass flux (per unit solid angle) into the heliosphere.
The initial relaxation is performed subject to two sets of boundary conditions at the radial outer limit of s max = 30R . In the case of set 1: the mass, momentum, and energy densities are linearly extrapolated across the boundary, and the pressure is systematically reduced in the ghost cells outside of the domain in order to encourage the development of a supersonic outflow. Once this outflow has been achieved the pressure Deposition Region h0[erg cm −3 s −1 ] s l [cm] Transition Region 2.9 × 10 −5 7.0 × 10 9 Lower Corona 8.7 × 10 −7 2.1 × 10 10 Extended Corona 5.0 × 10 −9 7.0 × 10 10 Table 2. Heating parameters for ion energy deposition. These are used for both the hydrostatic and transonic steadystate conditions and the dynamic reconnection runs. Electron heating is set to he = 0 in all cases.
reduction is removed and the temperature gradient is adjusted to prevent the formation of inward heat fluxes across the boundary. Following the transition of the outflow to a supersonic condition, a sonic point forms within the domain and migrates progressively inward until it stabilizes in the lower-middle corona at a height of roughly 7 R . The ultimate location of the sonic point is somewhat higher than might be expected owing to the fact that in our model the wind is accelerated entirely by the thermal pressure gradient, which results in somewhat slower wind speeds than if momentum is injected directly (Leer et al. 1982). Because the typical flow speed is of the same order as the sound speed in this configuration, the system achieves a quasi-steady transonic wind solution on the timescale of the acoustic travel time across the domain, which is of order 10 4 s. This quasi-steady transonic solution is depicted by the teal curves in Figure 1.
In the case of set 2: the mass, momentum, energy, and heat fluxes are all fixed to zero through the outer boundary, forcing the system to be closed. This causes the solution to settle toward a hydrostatic state in asymptotic time; because the slowest mode of the system is governed by the transport of mass, the relaxation time scales as the size of the domain divided by the fluid speed, which diverges as the fluid velocity approaches zero. The quasi-steady hydrostatic solution is depicted by the black curves in Figure 1. While the system never fully achieves a hydrostatic state, we consider it to be sufficiently relaxed when the maximum value of the coronal mass flux has fallen below 1% of the value of the transonic wind solution at the same location. This requirement is significantly more strict than the usual hydrostatic scaling of ρu 2 P (i.e., small Mach number), and is enforced specifically to ensure that the velocities of the hydrostatic and transonic solutions are well separated in the lower corona where even the transonic solution is significantly subsonic.
The equilibration time of the hydrostatic solution is extremely long (approximately 10 8 s) owing to our strict relaxation criteria and the length of the numerical domain, which is larger than any closed-field/hydrostatic configuration that might be found in the solar environ-  ment. The reason for constructing such a solution is to ensure compatibility between the numerical domains of the open-and closed-field configurations, and since it is only the lower portion of the plasma column that is needed from the closed-field initial condition, the unphysical radial extent of this solution is not a significant concern. Were the hydrostatic initial condition to be extracted from a curved field line with its apex in the lower or middle corona, the average volumetric heating rate along the entire field line would be increased, with the energy being distributed within a smaller volume, and the solution would likely be hotter and more dense, thereby exaggerating the differences already present in these two calculations.
For both of the solution profiles depicted in Figure 1 the location of the base of the transition region (where T e and T i begin to rise abruptly from their minimum values of 2 × 10 4 K) has fallen by a few hundred km from its initial location of 10Mm, with the final location being slightly lower for the hydrostatic solution than for the transonic solution. This settling is expected as the structure and location of the transition region are determined by the evolving mass and energy fluxes through the model chromosphere and lower corona, which change in time as the system relaxes toward an equilibrium state. Beyond this, the structure of the transition region and lower corona is very similar in both cases, with temperatures and densities that closely agree up to a height of ∼ 20 Mm above the transition region. Both solutions also exhibit oscillatory behavior in the velocity profiles up to a height of ∼ 100 Mm, above which height the transonic solution then begins to accelerate quickly into the middle-corona.

Dynamic Equilibria of Oxygen Ionization Levels
During the initial relaxation the ionization levels are not computed as these add substantial overhead to the calculation. Following the equilibration of the fluid profiles, we then initialize the ionization levels from a set of previously calculated equilibrium ionization solutions, which are parameterized in electron temperature and interpolated onto the numerical domain using the local plasma temperature. The two simulations are then advanced in time until the time-dependent ionization again converges to a steady state in each case.  We refer to the relaxed ionization profiles as "dynamic equilibria", to distinguish them from the assumed "ther-mal equilibrium" solutions that explicitly ignore advective transport. For the hydrostatic case the two solutions are nearly identical since there is no significant flow so the solution is characterized by an instantaneous balance of all ionization and recombination rates with For the transonic solution, however, the dynamic equilibrium is set by a competition between the instantaneous temperature of the fluid and the advective transport, which informs the temperature history of the fluid and, hence, the rate equations. The dynamic equilibrium ionization is characterized by (10) which exhibits spatially uniform ionization levels above the so-called "freeze-in" height (h f ), where the ionization and recombination timescale (τ j ) becomes long compared to the travel time across the temperature scale height so that In principle this freeze-in effect occurs separately for each ionization level; however, in practice the zero-sum nature of the rate equations and the limited number of ionization states with non-zero populations in the corona causes the effect to occur coherently across these states.
In the following sections we will refer to the ionization fractions Y j by their respective ionization states O j+ , with the understanding that these correspond to the underlying spatial-temporal populations. The final, quasisteady states of both the open-and closed-field solutions are shown in Figure 2, with both the initial (assumed) thermal equilibrium and relaxed dynamic equilibrium ionization levels shown for comparison.

Reconnection Model
While magnetic reconnection in three dimensions (3D) can occur throughout a distributed non-ideal region (Priest et al. 2003), the behavior that we have in mind for this study is the more discrete null-point reconnection (Pontin et al. 2013), in which the magnetic diffusivity is assumed to be negligible everywhere except within a very small volume surrounding an isolated coronal null point. A schematic of this process is shown in Figure  3, which depicts a 2D flux system with four magnetic domains that overlie a region of negative polarity flux within a larger region of positive polarity flux. This geometry is representative of a slice through a coronal pseusostreamer, as might be found above an isolated magnetic bipole within a larger unipolar region on the solar surface.
Two of the magnetic domains (I and II) are open to the heliosphere while the other two (III and IV) are closed, mapping from positive to negative polarity regions on the solar surface. The black field lines that partition these domains represent separatrices (or separatrix surfaces in 3D), which emanate from the magnetic null point (yellow dot). During interchange reconnection field lines from domains I and IV are swept into the null point where they reconnect and subsequently retract into domains II and III. As magnetic flux is processed through the reconnection site the plasma is similarly swept across the domain boundaries so that immediately after reconnection the plasma properties along each segment of the reconnected field lines are inherited from their respective source regions.
We construct the post-reconnection initial condition along a representative field line from segments of steadystate conditions for both open-and closed-field configurations, which we concatenate to form a piecewise steady state with a discontinuity at the reconnection site. This construction assumes that plasma in the inflow domains remains in a steady state, even as field lines are swept toward the reconnection site, so the tacit assumption is that the reconnection rate is sufficiently small for motion of the field lines themselves to have no effect on the field-aligned dynamics. Clearly this assumption ignores certain important aspects of null-point reconnection, including the expansion of the flux-tube cross section in the vicinity of the null-point; however, by suppressing these effects we are able to study the fluid-driven evolution more closely and with fewer confounding influences.
From the post-reconnection initial condition we again use the HYDRAD code to evolve the system in time.
To ensure a monotonic transition between the two fluid states on either side of the reconnection site, the constructed post-reconnection state is initially refined using linear interpolation and suppressing explicit conservation of the mass, momentum, and energy densities. Following this initial refinement, the subsequent integration is performed exactly as described in Section 2.1. We perform six calculations in total, using the same heating profiles as in the steady state calculations, with the reconnection site placed at H r ∈ {R /8, R /4, R /2, R , 2R , 4R }. Time-dependent ionization of oxygen is tracked in addition to the hydrodynamic evolution and state-files are output every 100 s of model time. The total simulation model time is 4 × 10 5 s (about 100 hours) from the instant of reconnection in each case, sufficient for the fluid to settle to III IV I II a quasi-steady state that is indistinguishable from the transonic initial condition.

Initial Riemann Decomposition
Because the numerical domains and magnetic geometries (s, g, and A) are identical between both the openand closed-field steady-states, the post-reconnection initial condition begins in equilibrium everywhere except at the reconnection site. There, the initial discontinuity subsequently evolves as a Riemann problem, decomposing into a shock and a rarefaction wave (as well as ion and electron thermal fronts), as needed for the jump conditions across each feature to collectively describe the total change in each of the fluid variables across the discontinuity. In particular, since the wind solution above the reconnection site exhibits a flow that is directed away from the hydrostatic condition below the reconnection site, the collective response is that of a rarefaction; however, because the wind is of lower density than the hydrostatic solution, any outward propagating feature must be compressive. The Riemann solution is, therefore, composed of a leading shock, which propagates outward at supersonic speed in the rest frame of the expanding wind, combined with a rarefaction wave whose leading edge propagates inward at the sound speed in the rest frame of the hydrostatic column below the reconnection site.
This initial evolution is depicted in Figure 4, for a reconnection event located at H r = 2R . The outward propagating shock and inward rarefaction are clearly visible in the velocity, which increases with height across the rarefaction and then decreases again across the shock, while the linear number density (nA) decreases with height across the rarefaction and then decreases again across the shock. Note that in the rest frame of the shock the fluid appears to be moving inward, so the low-speed/high-density post-shock region lies between the shock and the rarefaction, and this region grows in time as the two features propagate away from each other. As the rarefaction propagates inward it grows in time but continues to connect the quasi-static region below to the post-shock region above, so that while the leading edge propagates downward at the acoustic speed, the trailing edge remains nearly fixed at approximately the location of the initial discontinuity, which behaves as a stationary point for both the velocity and number density during this initial evolution.
The temperatures of the two species behave similarly to the number density during the initial Riemann evolution, but with additional dynamics resulting from viscous heating, thermal conduction, and compressive (adiabatic) heating and cooling. Because the thermal conductivity of the ions is relatively weak, their temperature evolves primarily in response to compressive and advective transport, so it tracks closely with the number density profile over time. However, the highly efficient electron thermal conductivity quickly dissipates the initial structure within the electron temperature profile creating broad heating and cooling fronts that extend far ahead of the shock and rarefaction, cooling the hydrostatic column below the reconnection site and heating the expanding wind. These temperature variations are visible in the adiabatic sound speed c 2 s = γ(p e + p i )/ρ (which reflects the average temperature of  Figure 4. Velocity, linear number density, and sound speed of the evolving Riemann solution for a reconnection event at Hr = 2R . The shock and rarefaction propagate away from the initial discontinuity, and are visible in the velocity and density profiles. The sound speed reflects the mean temperature of the two species, and exhibits dynamics resulting from both the hydrodynamic evolution as well as electron thermal conduction. the two species), as depicted in the lower panel of Figure  4.

Rarefaction Reversal
As the rarefaction propagates inward, the temperature below it is preconditioned by the reduction of the downward heat flux in the hydrostatic region below the reconnection site. This causes a weak downflow to develop in the region between the transition region and the rarefaction, which is quickly subsumed by the accelerating outflow across the leading edge of the rarefaction. This is visible in Figure 5 between t = 30 m and t = 60 m. Meanwhile, an outward mass-flux begins to develop through the transition region as the reduced temperature leads to a reduction in radiative cooling, which therefore requires an outward enthalpy flux (E + P )u in order to balance the local heating rate. Eventually the leading edge of the rarefaction stalls against the up-welling of mass through the transition region as it seeks a new equilibrium that is consistent with the electron temperature profile, which settles rapidly toward a wind-like solution. As the leading edge of the rarefaction stalls, the interior continues to expand over time, and the upper extent of the rarefaction accelerates outward, eventually becoming identifiable as the leading edge of a newly formed (reflected) outward-propagating rarefaction. . Velocity and linear number density during reversal of the rarefaction. Initially, the linear number density decreases with height through the rarefaction, consistent with inward-propagation; however, after reversal of the rarefaction the leading edge propagates outward and overtakes the leading shock, and the linear number density within the rarefaction increases with height.
The reversal of the rarefaction coincides with a reversal in the linear number-density profile, which subsequently increases with height (see discussion in the Appendix), and cannot (in the absence of a second shock) match the density jump required to connect the postshock outflow to the steadily decreasing density at the top of the transition region, which continues to fall as it seeks a wind-like solution. This effect, combined with the steadily increasing outward mass flux through the transition region, leads to the development of another (weak) trailing shock, which develops behind the nowoutward-propagating rarefaction. Meanwhile, the leading edge of the rarefaction propagates outward at the sound speed in the post-shock medium, which is itself subsonic in the rest frame of the leading shock. Eventually the outward-propagating rarefaction overtakes the leading shock, and the three features combine to form a single, coherent structure, known as an N-wave 1 (this being the ordered combination of two shocks connected by an interior rarefaction, as discussed by Friedrichs 1948).
The various stages of this reversal are depicted in Figure 5: the downward propagation of the left side of the rarefaction is clearly visible between t = 15m and t = 60 m, along with the monotonic decrease in the linear number density with height. From t = 120 m to t = 240 m the rarefaction accelerates upward, subsuming the post-shock outflow region below the shock, and the linear number density profile begins to flatten. Eventually, by t = 480 m, the upper extent of the rarefaction has overtaken the shock and the linear number density within the rarefaction now decreases inwardly away from the shock. This reversal in the linear number density is the earliest signature of the newly-formed N-Wave, which is characterized by velocity and density profiles whose gradients are codirectional through the entire structure.

N-Wave Structure and Dynamics
The structure and evolution of the outwardpropagating N-Wave are depicted in Figure 6, which shows the time-evolution of the fluid velocity, linear number density, and fractional ionization populations O 6+ and O 7+ . Note again that while n (not shown) decreases monotonically with height, nA necessarily increases with height in the region between the trailing shock and the leading shock, opposite the behavior of the initially-inward propagating rarefaction shown in Figure 4. The connection of the leading shock, through The velocity and density profiles mimic each other within the N-wave, which is composed of a leading shock, an outward-propagating rarefaction, and a trailing shock. The evolution of the ionization profiles reflects a competition between advective transport and spatially-varying ionization and recombination rates, which depend on the evolving temperature profiles.
the rarefaction, to the trailing shock is also visible in the velocity profiles, which exhibits the eponymous "N" shape, decreasing with height across the trailing shock, then increasing quasi-linearly through the rarefaction to a maximum at the leading shock, before then decreas-ing abruptly across the leading shock onto the pre-shock wind solution.
The ionization fractions, on the other hand, do not participate in the compressive and expansive dynamics of the shocks and rarefaction, despite the fact that the initial discontinuity in the values of Y j is co-spatial with the initial Riemann problem. This follows from the form of the advective component of the rate equation (i.e., the material derivative u ∂ s Y j ) which causes the ionization states to evolve as passive scalars in the absence of thermal effects (when the ionization and recombination rates become small) so that the ionization fractions are simply carried along by the flow but do not increase or decrease as the fluid expands. Accordingly, the initial discontinuity (which is above the freeze-in height in the case of the H r = 2R example) is carried along by the high-speed outflow within the rarefaction, but lags behind the leading shock, which propagates faster than the bulk flow. Thus, while the shock reaches the outermost boundary of 30R in just a bit less than 24 hr of simulation time, the ionization signature has only just passed 10R in that time, and does not arrive at the outer boundary until significantly later.
Below the freeze-in height, the fractional ionization populations of the initially hydrostatic fluid evolve in response to the changing temperature that occurs in both time and space as the fluid begins to expand within the rarefaction and eventually settles onto a new wind solution. For fluid parcels that originate sufficiently low in the corona, the transit time to the freeze-in height is long compared to both the timescales of fluid equilibration and ionization and recombination, so that by the time the fluid reaches the freeze-in height, the ionization ratios are indistinguishable from a steady-state wind solution. Fluid that originates higher up within the hydrostatic column has a thermal history that is more strongly representative of the conditions in the closed corona, having less time to equilibrate before it reaches the freeze-in height. The structure of the ionization profiles following reconnection is therefore given by an abrupt transition from the dynamic equilibrium of the initial wind solution just above the reconnection site to something resembling the dynamic equilibrium of the hydrostatic solution at the same height, followed by a smooth transition back to the quasi-steady wind profile, all of which travels outward at the local fluid speed.

Dependence on Reconnection Height
We explored the effects of reconnection height with six simulations that place the reconnection site in a variety of locations from R /8 to 4R . A comparison of these runs is shown in Figure 7, which depicts the evolution of the fluid velocity and oxygen ionization ratio O 7+ /O 6+ . In each case, the strength and speed of the leading shock and trailing rarefaction depend primarily on the size of the initial discontinuity, which is larger for reconnection sites that are higher in the corona. Therefore, reconnection events that occur near or above a height of R display an obvious velocity (and therefore density) signature, as this is the relevant height-scale for the Mach number of the wind solution to have appreciable size. For reconnection sites that are well below this height (H r ≤ 0.5R ) the initial discontinuity is too weak to create a strong shock and the fluid quickly settles onto a wind solution with little more than a short-lived transient wave disturbance. The initial temperature change across the reconnection site is less strongly affected by changes in H r ; however, for reconnection sites lower in the corona the ion and electron temperatures are more collisionaly coupled, so the dissipation of the initial temperature jump by the electrons similarly smooths the ion temperature profile and the reconnection signature becomes very weak in the far-field.
The ionization signature is similarly height dependent; however, in this case the relevant scale is whether the reconnection occurs above or below the freeze-in height, which occurs at roughly R /3. Reconnection events that occur above this height show a strong closed-field signature in the ionization ratios that is largely unaltered as it propagates into the heliosphere, followed by a slow decay back to the wind-like solution. For reconnection sites that are lower in the corona, the ionization ratios undergo more evolution and the signature of the material from within the initially hydrostatic column becomes progressively weaker, so that for reconnection sites well below the freeze-in height there is little-to-no signature, with the plasma having undergone significant thermal-temporal evolution before reaching the heliosphere. A critical difference, however, is that the ionization signature is not dissipated in the same way that the temperature and even velocity profiles can be, so whatever signature survives up to the freeze-in height is well preserved as it propagates into the heliosphere, as seen in the H r = 0.25 R panel of Figure 7. Note that while the relevant height for both fluid and ionization dynamics is of order R in these simulations, these heights are model dependent, and may change for different energy deposition profiles and ion populations, such as Fe (Z = 26), whose freeze-in height is likely to be higher than that of O (Z = 8). To explore the implications of these calculations for in situ observations we have extracted time series data for the fluid variables and ionization ratio O 7+ /O 6+ at s = 20 R (r = 21 R ) above the solar surface, as shown in Figure 8. This height is roughly coincident with the closest approach of Parker Solar Probe (PSP, Fox et al. 2016) during its sixth perihelion pass, and is at once far enough above the solar surface to be representative of the conditions in the heliosphere and also well within the numerical domain so as to avoid any possible boundary effects. Unfortunately, PSP does not carry instrumentation to detect ionization ratios; however, these are not expected to undergo further evolution beyond 20 R , so the values and timing depicted in Figure 8 should scale readily to heliocentric radii of r > 60 R , at which point they will be detectable by Solar Oribiter (SolO, Müller et al. 2020).

Relevance to In Situ Observations
In the figure, results from the six simulations are indicated by the variously colored curves, which run from dark purple to bright yellow for reconnection sites of in- creasing height. Consistent with the velocity profiles in Figure 7, the signatures in the density, temperature, and ionization ratios show consistently larger amplitudes as the reconnection height is increased. The arrival time of these signatures is also earlier for reconnection events that are higher in the corona, due in part to the reduced travel distance from the reconnection site to 20 R . For the most part, the fluid signatures arrive simultaneously, with the electron heat front being an exception due to the significantly larger electron thermal speed, which makes the transit time from the reconnection site to 20 R difficult to resolve on the scales represented here. For the remaining fluid properties, however, the arrival time is governed by the two remaining character-istic speeds -the sound speed (ion acoustic speed) and the bulk flow speed, both of which are on the order of a few 10 2 km s −1 over the majority of the spatial domain.
In the case of the density, velocity, and ion temperature, the primary signature is the leading shock, whose arrival time reflects the shock propagation speed from the reconnection site, through the lower and middle corona, and into the extended corona and heliosphere. This speed depends on the strength of the shock, being equal to the sound speed for weak shocks and as much as a few times the sound speed for strong shocks. Reconnection events that are closer to the transition region exhibit weaker/slower shocks that must propagate from lower in the corona, where the temperature and sound speed are correspondingly smaller. For that reason the signatures from reconnection events higher in the corona are not only larger in amplitude, but also arrive earlier and last longer by as much as a several hours, with the H r = 4R event spanning from t 12 − 30 hrs while the H r = R /8 event occurs over a much shorter time from t 20 − 22 hrs.
Nonetheless, the qualitative signatures are identical in each case, with an initial enhancement across the leading shock, followed by a slow decay through the rarefaction and finally a small jump across the weak, trailing shock, which arrives later for the higher-amplitude events, owing to the larger extent of the rarefaction. This causes the signatures from reconnection events that are higher in the corona to be not only stronger but also longerlived. Critically, the transient flow speed within the N-wave is systematically higher than the ambient wind speed, so that the net effect of interchange reconnection is to increase the speed of the wind relative to a steady-state solution.
For the ionization ratios we see a similar trend, with reconnection events placed higher in the corona having a stronger and more long-lasting signature, that arrives sooner and takes longer to decay. However, as we have previously shown the ionization ratios do not participate in the shock-compression dynamics, but are instead entrained behind the shock. As a result the ionization signatures arrive systematically later than the shock/rarefaction system, being typically delayed by several hours, enough that in most cases the fluid variables have returned to equilibrium values by the time the ionization signatures arrive. Within the ionization curves we see qualitatively similar behavior in all cases, with an initial enhancement of O 7+ /O 6+ reflecting the population that originated from within the hydrostatic region above the freeze-in height but below the reconnection site, followed by a slow decay reflecting the populations that passed through the freeze-in region before the N-wave had fully formed, and finally a return to the wind conditions from populations that originated near the base of the corona and experienced a thermal history that reflects the fully-relaxed wind solution.
From the study of Zhao et al. (2017), the ratio O 7+ /O 6+ in observed solar wind streams is typically in the range 0.05 to 0.2, being systematically lower for wind streams that are ballistically mapped from coronal hole regions vs active regions. By comparison, the steady-state wind stream that we simulate here exhibits an O 7+ /O 6+ ratio of ∼ 0.11, with enhancements following reconnection of up to 0.5 in the most extreme case. This suggests that reconnection events that occur lower in the corona and below the freeze-in height may be better matched to the observations, as would be the case for the H r = 0.25 R simulation, which exhibits a maximum O 7+ /O 6+ ratio of 0.2.
The relative timing of these signatures similarly reflects the amplitude of the N-wave, with larger shock rarefaction systems having faster flows in the rarefaction region, corresponding to earlier arrivals of the associated ionization profiles. The location of the reconnection sites relative to the freeze-in height is further evidenced in the flatness of the ionization curves between the arrival of the reconnection signature and the slow return to the wind-like signature, as this reflects the extent of the material in the column that originated below H r but above the freeze-in height. By comparison, for H r < R the initially hydrostatic populations have clearly undergone significant ionization and recombination on their way into the heliosphere, with the signatures of their origin in the closed-field domains being progressively weaker for lower reconnection heights.

Interpretation in 2+ Dimensions
To this point we have focused on the time-dependent evolution the plasma along individual field lines; however, if we return to the steady-state picture of interchange reconnection as depicted in Figure 3 we can imagine that all of the field lines in the outflow region(s) are undergoing the same dynamic evolution. By associating the field-perpendicular distance from the reconnection site with increasing time since reconnection we can map the various stages of the field-aligned dynamics along a given field line to different subdomains within the reconnection outflow region. A depiction of this mapping is shown in Figure 9, with the various hydrodynamic regions overlaid onto the magnetic field structure, which has been modified slightly from Figure 3 to reflect the collapse of the null-point into a thin current sheet.
Just as in Figure 3, there are four magnetic domains, two inflow regions (I and IV) and two outflow regions (II and III). The pre-reconnection transonic and hydrostatic fluid domains encompass the two magnetic inflow regions, but also extend into the outflow regions. This follows directly from the details of the field-aligned evolution and also reflects a broader causal connection between magnetic topology and plasma dynamics. Since field lines depend on the global structure of the magnetic field at a single instant in time, changes in magnetic connectivity occur instantaneously; however, plasma dynamics are dictated by local processes and signals that propagate at finite speed, so the effects of reconnection are limited to the causally connected region defined by the characteristics of the system (i.e., the leading shock and inward rarefaction).  Figure 9. Schematic of 2D interchange reconnection with overlaid hydrodynamic sub-domains. Solid orange curves represent shock fronts while dashed orange curves are weak discontinuities at the leading and trailing edges of the rarefaction wave. In the upper-right outflow region, near the reconnection site, the leading shock separates the pre-shock wind from the post-shock outflow above the inward rarefaction, which propagates down into the hydrostatic column. Farther from the reconnection site the rarefaction undergoes reversal, eventually subsuming the post-shock outflow. A second shock separates the now-outward rarefaction from the newly-formed wind below it, and together with the leading shock these three features make up the N-wave.
In the closed-field outflow region (domain III from Figure 3) the reconnection discontinuity connects a segment of the wind solution from the lower corona (domain I) with a hydrostatic plasma column from the closed-field region beneath the null point (domain IV). The resulting dynamics have not yet been simulated within our model; however, because the fluid velocity of the wind solution is directed toward the discontinuity, the resulting evolution will likely exhibit a pair of shocks that propagate away from each other and downward toward the solar surface. These would then rebound from the density gradient at the base of the transition region, driving chromospheric evaporation, siphon flows, and other interesting processes in the lower corona.
In the open-field outflow region (domain II from figure 3) the reconnection event connects an upper segment of the wind solution to a hydrostatic column below it and the dynamics along a given field line proceed as previously described. The various shocks and rarefactions along progressively older field lines form loci that generalize to magnetoacoustic disturbances. The leading shock and inward rarefaction along individual field lines generalize to slow magneto-acoustic fronts that separate the causally-disconnected portions of the outflow region in the far-field from the interior of the structure that develops from the initial reconnection discontinuity. The reversal of the inward-propagating rarefaction, which subsequently subsumes the quasi-steady shock outflow, and the formation of the outward-propagating rarefaction and trailing shock each define subdomains as depicted by the variously colored regions in the figure. The final N-wave structure is defined by the two near-parallel shock fronts and the outward-propagating rarefaction that separates them, with the wind speed being again enhanced relative to the steady state within the interior of the N-wave.
It is noteworthy that in both of the outflow regions (the open-field outflow that we have modeled here and the closed-field outflow that we have not yet addressed) the effect of the fluid discontinuity should be to generate a pair of shock-like fronts that expand as a wedge away from the reconnection site. This closely resembles the structure of a Petschek reconnection outflow, where a pair of slow magneto-acoustic shocks form in order to alter the effective aspect ratio of the current sheet to accommodate an increased reconnection rate over the classic Sweet-Parker scaling (Parker 1957;Sweet 1958;Petschek 1964;Kulsrud 2001). Here the shocks result from an asymmetry in the conditions in the two inflow regions, including a non-zero velocity jump across the reconnection site; however, it seems likely that the two are related and may in fact be limiting cases of the same underlying phenomenon. Other authors have previously observed the connection between one-dimensional Riemann problems and Petscheck reconnection, both symmetric and assymetric (see Lin & Lee 1999, and references therin), but these have generally focused on the cross-field dynamics in configurations where the asymptotic speed is zero in both inflow domains. Clearly, there is more work to be done in unifying these descriptions under more general initial conditions as the resulting behavior will depend on the overall reconnection rate and the extent of the asymmetry.

CONCLUSIONS
We have demonstrated how a post-reconnection plasma discontinuity naturally evolves as a shockrarefaction system in the manner of a Riemann problem, with a leading shock that expands into the heliosphere and an underlying rarefaction that propagates inward toward the solar surface. We have further shown how the initially-downward-propagating rarefaction reverses at the transition region to become an outwardpropagating shock-rarefaction system, which eventually overtakes the leading shock to form a shock-rarefactionshock triplet, or N-wave. This structure propagates coherently into the pre-existing wind in the heliosphere at a few times the local sound speed (or slow magnetoacoustic speed) after which the plasma behind it settles to a new quasi-steady wind solution on a time-scale that is comparable to the acoustic transit time from the solar surface to a given heliocentric radius.
A key finding of our model is that the release of material from the closed corona into the open field necessarily enhances the speed of the reconnected wind stream due to the interaction between the under-dense wind above the reconnection site and the high-pressure column below it. Because the fast solar wind is generally steady while the slow wind exhibits significant variability, our interpretation is that interchange reconnection is more common along field lines that support slow wind streams than those that support fast wind streams. This is consistent with the finding that the fast wind emanates from the interior of coronal holes while the slow wind originates from near coronal hole boundaries, where interchange reconnection is expected to be prevalent. However, if our model is accurate then the specific processes that result in the slower baseline speed of the wind that emanates from coronal hole boundaries are independent of the process of interchange reconnection, whose primary effect is the create intermittent enhancements of the slow wind speed relative to the steady-state.
Additionally, we have shown how the time-dependent ionization of oxygen evolves in the context of the postreconnection fluid, with the ionization ratios being carried along by the flow while also evolving subject to the rate equations in the co-moving fluid frame. We find that the location of the reconnection site relative to the freeze-in height is critical to the predicted temporal evolution of the ionization ratios for in situ measurements, with reconnection above the freeze-in height showing a strong signature of the initially hydrostatic plasma column that originates above the freeze-in height but below the reconnection site. As the reconnection site is moved progressively closer to or below the freeze-in height, the initially hydrostatic fluid undergoes more ionization and recombination in its journey to the heliosphere, and so this signature is weakened considerably.
These observations have strong implications for the strength of the signatures from reconnection events that occur near large magnetic structures, whose vertical extent is comparable to the source surface radius (i.e., H r ∼ 3R ) vs. more compact structures with H r R . Since pseudostreamers are usually on the lower end of this scale, it is likely that the strength of the discontinuity will be similarly on the lower end of the parameter space that we have described; however, because smaller structures are more likely to support the nearhydrostatic configuration on which our model depends, it is possible that the dependence on reconnection height will be weaker than these results suggest.
Despite this reduction in the strength of the ionization signal with decreased reconnection height, the timedelay between the arrival of the hydrodynamic and ionization signatures, which corresponds to the different travel times of the propagating signals (shocks and other waves) and the advected signals (material properties), seems to be robust across all reconnection heights. And while we have not calculated them here, it is likely that disparities in abundances between open-and closed-field plasma will follow a similar trend, being similarly transported by the flow. Confirmation of these features will depend on robust measurements of plasma dynamics, composition, and ionization in the inner heliosphere, and will require joint observation efforts from PSP and SolO, as the former will sample plasma in the appropriate source region, but only the latter possesses the suite of instruments required to make the measurement. However, because the plasma's material composition is unlikely to be altered between 20R and 200R (∼ 1AU) it should be possible to reconstruct these signatures under ideal circumstances.
The applicability of these results depends on the existence of pristine initial conditions both above and below the reconnection site and minimal cross-field dynamics during the subsequent evolution. The latter requirement will depend on the global magnetic evolution, which can be highly variable across disparate regions of the solar corona. However, these conditions are not unreasonable provided that the field-line drift velocity far from the reconnection site is not significantly greater than the acoustic speed and that the Lorentz force is sufficient to resist perpendicular gradients in the energy density of the fluid. That is to say that the coronal geometry should be relatively static and the plasma β and Alfvén Mach numbers should both be small.
The validity of the initial condition depends additionally on the structure and dynamics of the reconnection site, which must be compact and and well defined with adjacent open-and closed-field domains being unambiguously identifiable. In particular, it is important that the thickness of the reconnection layer should not be significantly larger than the ion mean-free-path and that the fluid undergo minimal evolution during the time that it takes for a newly-reconnected field line to emerge from the reconnection site. This will ensure that the fluid conditions on the newly reconnected field line appear as a discontinuity, which can then evolve in the manner that we have described. Accordingly, while the slow solar wind is expected to emanate from both pseudostreamers and helmets streamers, the ambiguous nature of the open-closed boundary along helmet streamers and the continuous opening and closing of loops within the heliospheric current sheet may preclude the application of this model in that context; although a similar construction could be applicable given suitable alterations to the model.
Assuming that the reconnection outflow speed is equal to the Alfvén speed, we can estimate the Alfvén transit time over the length of the reconnection site and compare it to the acoustic transit time over an ion meanfree-path length. Combined with our earlier requirement on the thickness of the reconnection layer we then find that the aspect ratio of the reconnection site should be less than 2/β, where β = 2c 2 s /v 2 a reflects the ratio of the acoustic speed to the Alfvén speed. These assumptions are not universally applicable throughout the solar corona; however, they are sufficiently general as to be applicable in certain cases, and we have described a 2+D geometry that represents the likeliest scenario for observing this phenomenon, that being the case of laminar/coherent reconnection across a null-separatrix system that is representative of a simplified coronal pseu-dostreamer. In this construction the expanding N-wave forms a wedge of high-speed wind within the rarefaction that is bounded by slow magneto-acoustic shocks at its leading and trailing edges. Depending on the structure of the magnetic field ahead of the N-wave this may provide a framework for the formation of magnetic switchbacks vis-à-vis the model of Schwadron & McComas (2021). Moreover, the structure of the reconnection outflow in this model is strongly reminiscent of Petscheck reconnection, with the acoustic shocks and rarefaction described here being naturally generalized to slow magnetosonic shocks and rarefactions.
In reality, the solar corona is always evolving, and that evolution will invariably complicate the behavior that we have demonstrated. The sensitivity of these results to the assumptions within our model remain to be tested; however, the dynamics that we have described here should persist in some form anywhere that the basic magnetic geometry is comparable to our model framework and the reconnection process is not too violent. As such, these dynamics should be viewed as a baseline for comparison so that the disparate effects of coherent reconnection and more complicated dynamic processes can be disentangled when interpreting in situ and remote sensing observations. 6. ACKNOWLEDGEMENTS RBS, SJB and MGL were supported for this project by NASA HSR grant 80HQTR21T0106. RBS and MGL were also supported for this project by the Office of Naval Research and by NASA PSP/WISPR grant NNG11EK11I. We thank the anonymous referee for their helpful comments and careful reading of the manuscript. Thanks go also to Peter Wyper for helpful discussions on the conditions near the reconnection site.