Heliospheric Compression due to Recent Nearby Supernova Explosions

The widespread detection of 60Fe in geological and lunar archives provides compelling evidence for recent nearby supernova explosions within $\sim 100$ pc around 3 Myr and 7 Myr ago. The blasts from these explosions had a profound effect on the heliosphere. We perform new calculations to study the compression of the heliosphere due to a supernova blast. Assuming a steady but non-isotropic solar wind, we explore a range of properties appropriate for supernova distances inspired by recent 60Fe data, and for a 20 pc supernova proposed to account for mass extinctions at the end-Devonian period. We examine the locations of the termination shock decelerating the solar wind and the heliopause that marks the boundary between the solar wind and supernova material. Pressure balance scaling holds, consistent with studies of other astrospheres. Solar wind anisotropy does not have an appreciable effect on shock geometry. We find that supernova explosions at 50 pc (95 pc) lead to heliopause locations at 16 au (23 au) when the forward shock arrives. Thus, the outer solar system was directly exposed to the blast, but the inner planets -- including the Earth -- were not. This finding reaffirms that the delivery of supernova material to the Earth is not from the blast plasma itself, but likely is from supernova dust grains. After the arrival of the forward shock, the weakening supernova blast will lead to a gradual rebound of the heliosphere, taking $\sim100$s of kyr to expand beyond 100 au. Prospects for future work are discussed.


INTRODUCTION
The local neighborhood of the Sun is an ever-changing environment, as a result of our residence in a starforming galaxy with a dynamic interstellar medium. The heliosphere must therefore evolve with time in response (Müller et al. 2006(Müller et al. , 2009Frisch & Slavin 2006;Frisch et al. 2011). Indeed, events on Galactic scales may have an impact on Earth. Terrestrial ice ages have been linked to our passage through the Galaxy's spiral arms (Gies & Helsel 2005), and our vertical motion through the Galactic disk may have affected the cosmic ray flux on Earth and have corresponding biological signatures (Medvedev & Melott 2007). The Sun's passage through a dense cloud could compress the heliosphere to 1 au or smaller (Yeghikyan & Fahr 2003; it has even been suggested that a recent such event could have occurred, and the corresponding compression could have ultimately had an effect on human evolution (Opher & Loeb 2022). In this work, we explore the hydrodynamic effects of near-Earth supernovae, in which the blast wave moves at a velocity much greater than the Sun's typical speed through the local interstellar medium (ISM).
There is now abundant evidence of multiple recent near-Earth supernovae within 100 pc. The most compelling is the discovery of live (undecayed) samples of 60 Fe (t 1/2 = 2.6 Myr). This radioactive isotope is found in geological records such as ocean sediments and crusts (Knie et al. 1999(Knie et al. , 2004Fitoussi et al. 2008;Wallner et al. 2016;Ludwig et al. 2016;Wallner et al. 2021), Antarctic snow (Koll et al. 2019), and even lunar samples returned by Apollo astronauts (Fimiani et al. 2016). 60 Fe has also been found in cosmic rays (Binns et al. 2016), which show a low-energy excess Fe flux that could be evidence for a recent nearby source (Boschini et al. 2021). Wallner et al. (2021) also detected 244 Pu coincident with the 60 Fe signals. In addition, Korschinek et al. (2020) reports a detection of 53 Mn, and 26 Al has been studied but not yet separated from the overwhelming terrigenic component (Feige et al. 2018).
Further evidence of a nearby supernova comes from the proton, antiproton, and positron cosmic-ray spectra and anisotropy (Kachelrieß et al. 2018;Savchenko et al. 2015), the existence of the Local Bubble (Smith & Cox 2001;Frisch & Dwarkadas 2017), and the observed distribution of Galactic 26 Al (Fujimoto et al. 2020). By tracing back the motion of runaway stars and neutron stars, Neuhäuser et al. (2020) suggested stars in a binary progenitor system for these supernovae; similarly, Tetzlaff et al. (2013) also suggested a progenitor for the nearby Antlia supernova remnant (McCullough et al. 2002). The short 60 Fe lifetime demands that it was produced recently and thus nearby. These data are consistent with a nearby supernova ∼3 Myr ago, and the 60 Fe abundance implies a supernova distance of 60 − 130 pc (Fry et al. 2015).
The observed 60 Fe came to us in the form of dust; 60 Fe ions would have been a component of the plasma that gets deflected by the heliosphere. The dynamics of dust grains in the outer heliosphere have received considerable theoretical study (e.g., Belyaev & Rafikov 2010;Sterken et al. 2012), especially for the presentday heliosphere. Wallis (1987) found that during the passage through a dense cloud, dust can penetrate the heliosphere to Earth with little deflection due to the heliosphere's small size. Athanassiadou & Fields (2011);Fry et al. (2016) found that dust grains from near-Earth supernovae are typically deflected less than 1 • by the heliosphere, due to their high speeds v dust v esc (1 au) = 42 km/s far exceeding the escape speed at 1 au. The simulations we perform here may also contribute to our understanding of dust grain dynamics from near-Earth supernovae. Wallner et al. (2021) has recently found 60 Fe from a second pulse due to an earlier supernova ∼7 Myr ago, showing that nearby supernovae are relatively commonplace on geological and astrophysical timescales. These two known events are at roughly similar distances, too far to cause mass extinctions of species on Earth, although possible damage to the biosphere is an open question under study (Melott & Thomas 2019;Melott et al. 2017; Thomas et al. 2016). Closer events should occur, but less frequently. With this in mind, (Fields et al. 2020) proposed that one or more supernovae at ∼ 20 pc could have triggered extinctions at the end of the Devonian period 360 Myr ago, leading to observed global ozone depletion reflecting ionizing radiation damage from the explosion.
Motivated by these data, we consider the case of a near-Earth explosion, one of the most dramatic events the heliosphere can experience. This scenario is relevant for both the distant past and recent well-documented events. In the aftermath of a supernova, the supernova remnant (SNR) rapidly expands outwards, sweeping up the interstellar medium and eventually engulfing many surrounding stars. As the blast wave encounters these stars, it drives back their stellar winds, compressing their astrospheres. We aim to study the extent of this compression as applied to own heliosphere with a suite of numerical simulations to determine the innermost distance the supernova blast penetrates in our solar system.
To date, the only simulations of supernovae interacting with the heliosphere has been done in Fields et al. (2008), which lays the foundations for our work here. This earlier study examined the impact of supernovae out to at most 30 pc, closer than current estimates suggest for the 3 Myr event Fry et al. (2015). Our work will for the first time study blasts from supernovae out to 126 pc, in line with the results from analyses of the 60 Fe data. We also perform detailed comparisons of our solar wind model against in situ measurements from Voyager 2 and Ulysses. In addition, we study the scaling of the heliosphere dimensions with the supernova blast properties in a more detailed manner. We then use these scalings to estimate the evolution of heliosphere compression with the arrival of the supernova shock and subsequent rebound towards its present boundary.
The structure of this paper is as follows: section 2 describes the formalism and initialization of the simulations as well as expectations; section 3 presents the results of the simulations; section 4 discusses these results; and section 5 gives concluding remarks.

SIMULATION MODEL
Our goal is to examine how our heliosphere is compressed by a supernova blast wave. Since the most recent time this occurred was ∼3 Myr ago, present-day observations of this phenomenon are impossible. Therefore, we must turn to numerical simulations. We begin with the basic fluid equations, and then show that they produce a solar wind that roughly agrees with observations. Next we apply the Sedov model for a supernova remnant as the input for the blast wave. Finally we explore scaling laws for how these flows should interact.

Fluid equations
The modern heliosphere enjoys a variety of complex physics that deviates from standard hydrodynamics, such as magnetic fields, multi-fluid flows, charge exchange, and cosmic ray propagation (Pauls et al. 1995;Zank 1999). Our goal, however, is to investigate the broad changes induced by the supernova blast. Consequently, we do not attempt to compete with the sophisticated models that include these effects, such as those presented in, e.g., Pogorelov et al. (2004), Izmodenov et al. (2008), and Opher et al. (2015. Indeed, some of the relevant physics for the modern heliosphere may not be applicable for the case of an incoming supernova blast. For example, charge exchange occurs when neutral ISM atoms penetrate into the heliosphere (Baranov & Malama 1993;Pauls & Zank 1997). As a result, the solar wind's ram pressure is weakened and the boundary between the solar wind and the present-day ISM is closer than it would be for a fully ionized ISM. But as we will show, we do not expect a SNR to contain a large population of neutrals, in which case the effects of charge exchange are not at play. Out to 30 au, the one-component model of the solar wind matches the observed density very well and the velocity to a difference of less than 10% (see Fig. 4.2 in Zank 1999), a result we will confirm below.
To start, we assume that both solar wind and SNR flows can be adequately described with the basic equations of hydrodynamics, where p = nkT for an ideal gas. We use an adiabatic equation of state with γ = 5/3. We solve these equations with the Athena++ code (Stone et al. 2020). Athena++ is a grid-based magnetohydrodynamics framework, though we only make use of its hydrodynamics and neglect magnetic fields. We use the built-in HLLE Riemann solver, well-suited for our case where the kinetic energy of the flow dominates. While this solver can be diffusive, especially around contact discontinuities, it suppresses the Carbuncle instability and the "odd-even decoupling" that can appear when using other solvers (Sutherland et al. 2003;Quirk 1994). While performing our simulations, we do not see evidence for substantial diffusion around any contact discontinuities.

Solar wind initialization
The solar wind is the complex flow of gas launched from the solar corona and streaming outwards through the solar system, first predicted by Parker (1958). While the real solar wind varies with time according to solar activity (Provornikova et al. 2014;Izmodenov et al. 2008), for these simulations we adopt a constant, steady outflow. Spacecraft near Earth's orbit such as ACE and DSCOVR have taken an abundance of solar wind data. Our main interest is in how the solar wind interacts far away from the Sun, so we input the wind at 1 au using a rough average of solar wind density, speed, and thermal pressure. All grid cells within 1 au are overwritten with constant values every timestep.
The real solar wind not only varies in time, but also location: it is launched differently in the plane of the solar system than towards the poles. We use data from Ulysses's close approach in 1995 to approximate these parameters as a step function in angle, shown in Table  1. For comparison, we also show the values adopted in Fields et al. (2008), which is spherically symmetric.
While a reasonable wind initialization is made across the grid at the start of the simulation, the solar wind is allowed to relax to its steady-state profile before the supernova blast is introduced.
With our steady one-fluid model, we must check to ensure our simulated solar wind still accurately represents the observed heliosphere. To this end, we compare our relaxed solar wind in the equatorial region to Voyager 2 plasma data. Figure 1 shows density, velocity, and thermal and ram pressures of both our model and Voyager 2 data. A 240-day running average of Voyager 2 data is also shown. Our model tracks Voyager 2's density, velocity, and ram pressure well, though not accounting for temporal variations. The thermal pressure clearly has a different profile. We attribute this in part to the assumption of the adiabatic evolution of our model, whereas the real solar wind is not adiabatic. Since the thermal pressure approximately matches in the region from 5-40 au, we consider it to be sufficiently accurate in the most relevant region. Furthermore, the thermal pressure is still several orders of magnitude below the ram pressure, which will dominate the large-scale structure of the simulations.
In Figure 2 we plot the velocity as a function of heliolatitude during Ulysses' close approach. Although initially a polar step function, the abrupt change in velocity is slightly smoothed out as the wind propagates. Data from our model is taken from a distance of 2 au in order to allow the solar wind time to relax and give a better description of the distant solar wind than the input step function. The largest discrepancy here is due to temporal variability, which we do not model. The averages over time are a close match.
Currently, it is an open question as to where the supernova 3 Myr ago exploded, and therefore which direction the blast wave would come from. Proposed clusters include the Sco-Cen association (Benítez et al. 2002) and the Tuc-Hor association (Mamajek 2015;Hyde & Pecaut 2018). Sørensen et al. (2017) traced back nearby stellar clusters and statistically examined which ones were most likely to produce a nearby supernova explosion. In addition, this supernova could have been a companion to a previous star that exploded, and could have been flung outwards with orbital velocity and exploded away from its home cluster. Given these uncertainties in location, we do not assume a single place of the supernova, Figure 2. Comparison of our solar wind velocity to measurements made by Ulysses during its close approach of the solar minimum of 1995. We use the same color scheme as in Figure 1. Solar wind velocity does not change appreciably with distance ( Fig. 1) but simulations are for a distance of 2 au; data follow the Ulysses elliptical orbit and sample distances from 1.3 to 2.8 au.
but instead examine the effect of solar wind orientation with respect to the supernova. In order to increase computational efficiency, we perform our simulations in 2D cylindrical coordinates in the r-z plane, with rotational symmetry imposed about the z axis. We use two different orientations depending on whether the blast wave enters orthogonal or parallel to the plane of the solar system. Figure 3 shows a schematic geometric interpretation of these orientations. When the blast comes perpendicular to the axis of rotation as shown in panel (a), the rotational symmetry of our simulation captures the solar wind behavior, which we model with the step function in angle described above; this is the polar orientation. When the blast comes in along the plane (striking the outer planets' orbits first) as in panel (b), we use a spherically symmetric wind; we refer to this as the equatorial orientation. Both orientations have similar ram pressures, so we do not expect great differences in heliospheric structure between the orientations.
Our base mesh resolution is 256×512 cells for all but the comparison to Fields et al. (2008), which is 1024×1536. If kept uniform across the grid, several of our larger simulations would make the region within 1 au only a few cells across, resulting in a poorly-defined flow. To ensure a smooth spherical outflow, we refine the solar wind injection region for the entirety of the simulation (known as static mesh refinement, SMR). The amount of refinement depends on the outer boundaries of the mesh. Figure 4 shows the the relaxed solar wind and corresponding meshblocks in the polar orientation. Each refinement level increases resolution by a factor of 2.

Supernova blast initialization
Supernova remnants evolve over time. After the initial free expansion phase, supernova remnant morphology roughly follows a Sedov-Taylor profile (Sedov 1946;Taylor 1950) for much of its evolution. The remnant expands quasi-spherically over many parsecs until thermal emission becomes important. For SNRs several tens of parsecs across, a spherical forward shock is wellapproximated as plane wave on the scale of the solar system. Since the most important parameter in our simulations is the distance to the supernova, R SN , we invert the usual Sedov equation for distance to solve for time as where R SN is the radius of the remnant, ρ 0 is the ambient medium density, β is a numerical factor of 1.1517 for the typical γ = 5/3, and E SN is the energy of the supernova, taken to be 10 51 ergs. By applying the Rankine-Hugoniot conditions, we find the density, velocity, and thermal pressure of the gas immediately postshock, which are where the subscripts 0 and 1 indicate the ambient medium and immediate post-shock gas, respectively. γ is the adiabatic index, and v s is the shock velocity, given by v s = 2 5 Our current local interstellar environment is dominated by the Local Bubble, a region of low-density carved out by multiple supernovae (Smith & Cox 2001). While the density varies greatly with location, we use an average value of n = 0.005 cm −3 , rather than a typical Galactic ISM. The earliest supernovae to explode did so in a denser environment, and later ones encounter the swept-out low-density region from the past explosions. (Fuchs et al. 2006) estimated that 14-20 supernovae have exploded in the Local Bubble, and (Schulreich et al. 2017) and Breitschwerdt et al. (2016) performed hydrodynamic simulations showing how 16 supernovae can generate this environment. We will show that the ambient density has no effect on the distance of closest approach in our solar system, a feature of using the Sedov model. But the low density of the medium extends the duration of the Sedov phase, which should be a reasonable approximation for the supernova distances of interest here.
Furthermore, we can calculate the temperature of the post-shock material assuming an ideal gas. This yields a remarkably high temperature of For the low density observed in the Local Bubble, this model predicts that the temperature post-shock will remain over 10 6 K for 100 pc, the approximate size of the Local Bubble. In the past, it was accepted that Local Bubble consisted of a large hot component with T ∼ 10 6 K, though such assumptions have recently been challenged (Welsh & Shelton 2009;Linsky & Redfield 2021). Despite the hot environment of supernova remnants, the quick formation of dust and even molecules in SN 1987A (Matsuura et al. 2017) suggest that some component of the remnant may be neutral. We do not account for multi-fluid hydrodynamics here, instead interpreting the Sedov equations to indicate complete ionization at the forward shock.
The most important parameter for the supernova is the distance, which directly affects the strength of the blast wave, with ram pressure scaling as P ram ∼ E SN /R 3 SN . The supernova 3 Myr ago is estimated to have occurred somewhere between 60-130 pc away (Fry et al. 2015). Recently, Fields et al. (2020) proposed that a supernova at ∼ 20 pc could have contributed to biological extinctions at the end-Devonian period. To cover this range of distances, we use supernova distances of 25.3, 50.0, 63.3, 75.9, 94.9, 110.0, and 126.5 pc. Looney et al. (2006) proposed that a very nearby supernova within ∼ 1 pc exploded in the early stages of the solar nebula; we limit our simulations to the fully-formed solar system and so do not consider this case here.
From a single point in space, SNRs weaken over tens of thousands of years. Our simulations cover at most the first few years of the initial blast; this timescale is a multiple of the ∼ 0.5 yr crossing time for a flow of 100 km/s to travel 10 au. Therefore, we do not allow the remnant to weaken during our simulation. Since we are primarily interested in the closest approach of the blast wave in our solar system, we restrict ourselves to modelling the strongest part of the blast, the forward shock. See §4.2 for discussion of the behavior over longer timescales.
The supernova blast is implemented as a boundary condition in Athena++. We first let the solar wind evolve until it has reached a steady state across the entire domain. Then the boundary condition at −z max is changed to be the incoming supernova blast, flowing in the +z direction. This new boundary condition is kept constant for the duration of the simulation. All times shown in our figures define t = 0 as the introduction of the blast.

Expected heliosphere structure and stagnation distance
Our work draws upon insights from the many studies of the present-day heliosphere's interaction with the very local ISM (e.g., Opher et al. 2020;Frisch et al. 2011). In addition, several studies of stellar winds interacting with the ISM have been carried out, both observational and numerical (Meyer et al. 2021;Henney & Arthur 2019;Kobulnicky et al. 2016). These studies show that the solar wind-ISM interaction region consists of three main features: the termination shock (TS), where the solar wind slows to a subsonic velocity, the heliopause (HP), where the solar wind and in-flowing ISM meet, and the bow shock (BS), where the ISM transitions from supersonic to subsonic. In the modern day, the Voyager 1 and Voyager 2 missions have passed the HP at distances of 121.6 and 119.0 au, respectively (Burlaga et al. 2019). There is increasing evidence that there is no BS (Mc-Comas et al. 2012), indicating that the Sun's motion through the ISM is barely subsonic. This is not the case for a very rapidly-moving supernova blast.
We expect that the closest approach of the blast wave, directly on-axis, will be the point of pressure balance. (As we will show in section 4, pressure balance is instead an excellent predictor of the TS rather than the HP.) This point, also called the stagnation distance, is relevant for astrospheres and their bow shocks (Wilkin 1996;Comerón & Kaper 1998). While commonly written as a function of a star's mass loss rate as where v w is the stellar wind velocity and v * is the speed of the star relative to the ISM (Comerón & Kaper 1998), our model of the Sun has both an equatorial and polar region with slightly different mass loss rates. In order to relate the stagnation distance to our input parameters, we write it as a balance of thermal and ram pressure between the Sun and the supernova remnant. This stagnation distance is where solar wind properties are evaluated at 1 au. Fields et al. (2008) found good agreement with this (for the TS) for very nearby supernovae. Given that the SNR gas parameters depend completely on distance, r stag could equivalently be written in terms of the Sedov supernova distance as = A R SN 100 pc where A = 24.97 (23.50) au in the equatorial (polar) orientation. This final equation allows us to write the stagnation distance solely in terms of the supernova distance, given an explosion energy.
Qualitatively, different supernova distances should yield the same large-scale heliospheric shape and features among the simulations. Nearby supernovae will produce a much larger velocity than distant ones, which may affect the production of Kelvin-Helmholtz instabilities that form on the HP.

RESULTS
We run 13 simulations, numbered with respect to supernova distance. Model 1 is a comparison to model 12 of Fields et al. (2008) for a 20 pc supernova. Here, the solar wind follows their parameters, rather than the updated ones for the rest of our models. Models 2-4 are for a 25.3 pc and 50 pc supernova in the equatorial and top-down orientations. Models 5a, 5b, and 5c are all a 63.3 pc supernova in the equatorial orientation, but with higher ambient medium densities representing a more dense Local Bubble before multiple supernovae carved it out. Models 6-11 are for supernovae at larger distances in both equatorial and polar orientations. A summary of the results is given in Table 2, which lists the initial conditions for the supernova and the distances of closest approach for the TS, HP, and BS. The location of the BS is not stated for the cases where it has retreated off the grid domain.
In Figure 5 we show three density plots from the last timestep of models 6, 8, and 11. We see the qualitative expected upwind structure including the TS, HP, and BS. Any turbulence happens across the HP, seen most easily in model 8. The high-density equatorial solar wind gets bent back and is incorporated into the rest of the flow. It does not appear to be a source of turbulence. Instead, the HP drives Kelvin-Helmholtz instabilities. Figure 5 and all of our simulations show that for supernova distances consistent with recent 60 Fe measurements, the closest approach of the blast is > 10 au away. Certainly the blast can only arrive at 1 au for supernovae at extinction-level distances. But Fry et al. (2015) showed that 60 Fe abundances measured in terrestrial and lunar archive over the past 10 Myr imply a supernova distance of 50-100 pc. This reaffirms the conclusions of Fields et al. (2008) that these radioisotopes must arrive in the form of dust that can decouple from the blast at the supernova-solar wind interface. The dust must then travel ∼ 10 au, which requires large or fast grains to avoid repulsion from the solar light pressure (Athanassiadou & Fields 2011;Fry et al. 2016).

Comparison to previous work
In model 1, we use the same input parameters as Fields et al. (2008) model 12 in order to compare our two codes. We find excellent agreement with the overall structure of the heliosphere and the locations of the TS and HP upwind. The largest difference between the two is the amount of Kelvin-Helmholtz instabilities at the HP present in the previous work. We attribute this abundance to the use of adaptive mesh refinement in the previous work. This refinement was not used here in favor of static mesh refinement that resolved the solar wind at 1 au over much larger scales. While further refinement is needed to study instabilities and downwind mixing in detail, the two results are otherwise consistent.

Time-dependent features
A characteristic timescale for these simulations is the time for the blast wave to cross the simulation domain, t cross . When t ∼ t cross , the heliosphere directly upwind is in its most compressed state. Over multiple t cross , several processes continue to shape the heliosphere: Kelvin-Helmholtz instabilities along the HP, rebounding of the HP and BS, and downwind TS splitting.
The HP is a tangential discontinuity formed by solar wind and supernova fluids flowing along a surface of contact with a velocity tangential to this surface. This condition is ripe for Kelvin-Helmholtz instabilities, as we see in our simulations. Spiral features begin close to the axis of symmetry and grow as they are pushed downwind. The appearance of these instabilities is, in part, due to the resolution of our simulations, as highresolution runs promote more instabilities and greater mixing. Since our goal is primarily to determine the innermost penetration of the SNR, we do not extend additional layers of resolution to the Kelvin-Helmholtz ripples.

Discontinuity locations
Tracking the location of the discontinuities (TS, HP, and BS) over time can be used to determine their stability. The Mach number M is an excellent indicator of the location of strong shocks. Both the solar wind and blast wave are initially supersonic and must decrease below M = 1 in order to interact head-on along the axis of symmetry. We note that as the blast first enters the domain, the BS does not immediately develop. Instead, there is a smooth gradient over M until it sharpens into a single location after approximately t ∼ t cross . Figure 6 shows the Mach number along the z-axis for models 5a, 8, and 11. In this plot, the locations of the TS and BS are clear, marked by the near-vertical jump in the Mach number as it crosses M = 1. In 5a, the BS has receded out of the simulation frame, so is only marked by a small uptick at the domain boundary.
The HP is located where the two flows meet, which should be near-zero velocity along this axis. Models 5a and 11 clearly show this sharp downwards spike. The HP in model 8 is more difficult to locate. In this simulation, Kelvin-Helmholtz instabilities form close to the axis of symmetry. These ripples cause the location of the HP to change each frame, even after several t cross have passed. We attribute this to the polar orientation of both models: the polar solar wind is lower density, and the greater density contrast appears to promote the growth of Kelvin-Helmholtz instabilities close to the axis of symmetry. Nonetheless, we still use the minimum value for M as the HP with the understanding that this value has some error of ∼ 10% for models 8 and 10. We can apply this Mach number analysis over the duration of the simulation to examine the evolution of discontinuities over time. Figure 7 shows the locations of the TS, HP, and BS for model 6. Upon reaching the closest approach at t ∼ 250 days, the TS remains extremely stable. The HP has some small motion, mostly due to the effect of Kelvin-Helmholtz instabilities. The BS approaches an innermost position, then slowly retreats upstream over the duration of the simulation. This is a numerical artifact, largely due to an unintended interaction of the BS with the outer r-axis. By widening the boundary to capture the full extent of the shock so that the BS falls off the +r-direction, we find that the BS retreat is not as dramatic. However, doing so is both computationally more intensive and does not affect the other heliosphere features.

Downstream features
When the TS meets itself at a point in the downwind side of the simulation, it splits into three features labelled in Baranov & Malama (1993). These are (from the inside outwards) the Mach disk, the secondary tangential discontinuity, and the reflected shock (as shown in Figure 5. In model 1, we even see the appearance of Kelvin-Helmholtz instabilities in this secondary tangential discontinuity. The appearance of these features is a qualitative way to verify the accuracy of this code, as they are expected in the case of a fully-ionized ISM. More detailed models that include a neutral ISM component do not see evidence of these features.

Local Bubble density
For our Sedov supernova blast, we use a uniform ISM density of n = 0.005 cm −3 . Of course, the Local Bubble does not have that same density everywhere, as evidenced by the Complex of Local Interstellar Clouds. Due to self-similarity in the Sedov blast, the ram pressure does not depend on the ambient density. Therefore, the distance of closest approach should not depend on the ambient density. This relation is also reflected in eq. (15), in which ambient density does not appear in the expression for stagnation distance. To verify this relation and examine any other differences beyond stagnation distance, we ran three simulations changing only the ISM ambient density in models 5a-c. The ambient densities used here are 3.34 × 10 −26 , 1.34 × 10 −25 and 6.68 × 10 −25 g cm −3 (i.e., ρ 0 , 4ρ 0 , and 16ρ 0 ) for models 5a, 5b, and 5c, respectively.
We show these simulations after the same amount of time has passed since the introduction of the blast in Figure 8 (all are in the isotropic equatorial orientation). The most significant output of these simulations, the distance of closest approach, remains unchanged independent of the ambient ISM density, as expected. Kelvin-Helmholtz instabilities are seen in the low density simulation, but are not as apparent in the higher density cases. Larger ambient densities slow the blast, resulting in a smaller velocity difference that hinders the growth of these instabilities. Other features are largely the same, though the smaller blast wave velocity means evolution takes place on a longer timescale. Qualitatively, in Fig.  8(c), we see a clear bow shock and lack of features in the downstream termination shock, but these will continue evolving as the simulation advances.
We conclude that the ambient density does not have a large impact on our simulations. The greatest effect is on the timescale of the heliosphere compression, with a larger ambient density slowing the compression. Larger ISM densities will also increase the time for the blast wave to reach the solar system after the inciting supernova explosion.

Orientation effects
These simulations are performed in two orientations depending on the location of the supernova compared to the equatorial solar wind. The polar wind has a lower density but higher velocity, making the ram pressures similar (see Table 1). The overall structure of the heliosphere is the same for both orientations. The main difference is that, in the polar orientations, the equatorial wind is bent back through the heliosheath. It does not appear to be a source of instabilities, but reacts to those produced along the HP, as seen in Figure 5(b).
We ran three polar orientations with equatorial counterparts (numbers 3, 6, and 10; model 8 is polar but does not have a corresponding equatorial simulation at the same distance). As expected from the weaker ram pressure, the polar orientation compresses the heliosphere more. This ∼10% difference in the TS, however, is dwarfed by the different values for the supernova distance. We conclude that the supernova distance is more important than its orientation for heliosphere compression.

DISCUSSION AND ANALYSIS
These simulations of a supernova blast wave colliding with the heliosphere assume idealized hydrodynamics. Although the real situation is surely more complex, we expect that the gross features of the perturbed heliosphere are captured in our simulation. The apparent lack of a large neutral component to supernova blasts implies that our single-fluid treatment is a reasonable approximation. Thus, scaling laws and analytical arguments can be formed.
If one assumes a thin region for the heliosheath, the position of the bow shock should follow a simple analytical expression as a function of angle. As derived by (a) (b) (c) Figure 8. Density plots of models 5a, 5b, and 5c for supernovae 63.3 pc at the same time for three different ambient ISM densities: 3.34 × 10 −26 , 1.34 × 10 −25 and 6.68 × 10 −25 g cm −3 , from left to right. An animated version of (a) is available online showing the arrival of the blast wave into the solar system at time t = 0 and the subsequent compression of the heliosphere until the time of the displayed frame. Wilkin (1996), the position is Our simulations do not show a thin heliosheath, instead showing a well-separated TS, HP, and BS. Accordingly, this equation is only accurate for small values of θ directly upstream.

Discontinuity Scaling
Given the final locations of the TS, HP, and BS in each of these simulations, relations between these values and the supernova blast waves can be examined. This comparison will test the usefulness of r stag scaling.
In Figure 9, we compare the location of the TS and HP to two parameters: the stagnation distance from Equation 13 in (a) and the supernova distance in (b). When comparing to r stag , the TS has a very tight relation along the y = x line. According to this relation, r stag is an excellent predictor of the location of the TS rather than the HP. This effect is also noted in, e.g., Comerón & Kaper (1998) and Comerón & Pasquali (2007). The location of the HP appears to be dependent on orientation, approaching much closer in the polar orientations (particularly models 8 and 10) than in the equatorial one. We note that these two models also suffered from Kelvin-Helmholtz instabilities near the axis of symmetry. It is possible that the increased grid resolution of these models is responsible for this difference, though it does not affect the models in the equatorial orientation similarly.
Similarly to Figure 9(a), we can plot the TS and HP distances as a function of R SN , shown in Figure 9(b). We also plot eq. (15) for both equatorial and polar orientations, though the small difference between these two lines emphasizes how solar orientation is a less significant parameter than supernova distance. We again see the tight correlation between the TS and r stag , as well as the very apparent r stag ∝ R 3/2 SN relation. The locations of the HP shown in Figure 9(a) seem to suggest a correlation between TS and HP distances. To examine this, we plot these two distances against each other in Figure 10. We fit a line to all these points, forcing it to go through the origin. The slope of the best-fit line is 1.389 ± 0.067, in very good agreement with Fields et al. (2008) who found a slope of 1.41. The smaller slope we find is primarily due to models 8 and 10, which lie well below the line. As discussed above, this departure from the trend is likely due to the increased resolution, allowing instabilities to form close to the axis of symmetry. The other models show decent agreement with the fit, though more points would show how well this fit holds over a larger range. This relation does not hold for the present-day heliosphere, where the TS and HP are ∼100 and ∼120 au, respectively. The derived relation comes from purely fluid dynamics, and the present-day heliosphere requires a variety of more complex physical process to model correctly. Therefore, the present heliosphere is not required or even expected to fit among this trend.

Blast weakening over time
The strength of the supernova blast will weaken over time as the remnant expands. The duration of this process depends on the explosion distance and local density, but certainly takes 1 kyr. This timescale is far longer than the ∼ 1 yr duration of our simulations. Accordingly, supernova blast properties change on similarly long timescales. We can study the effect of the weakening supernova blast by considering "snapshots" of the blast properties at different times.
We have shown that r stag is an excellent predictor for the distance to the TS. By taking advantage of this relation, we can extend this analysis over the whole passage of the blast rather than just the leading edge.
The Sedov-Taylor model for supernova remnant evolution can be used to calculate the pressure balance distance for longer time scales. Normally, the Sedov solution is solved in terms of the outermost blast wave. In this instance, we wish to solve for the gas parameters for a stationary observer at a constant location from the origin.
We employ a Sedov blast wave verification code 1 to calculate the Sedov profile over the first 300 kyr after the blast. The chosen ambient density is still that of the Local Bubble, n amb = 0.005 cm −3 . At each timestep, we obtain the thermal and ram pressure and use Eqn. 13 to calculate r stag (indicating the TS position). We show the results of these calculations in Fig. 11. The vertical lines correspond to the initial arrival of the blast; their spacing in time reflects the duration of the blast travel to the Solar System for different supernova distances (eq. 4). Note that on the timescales plotted, our hydrodynamic simulations cover a thin ∼ 1 yr at the innermost region upon the blast arrival; the rest of the curves use pressure balance to find the closest approach of the supernova. Maximum heliospheric compression lasts on the order of ∼kyr, however, it takes > 100 kyr for the blast to weaken enough so that the heliosphere fully rebounds. Figure 11 has important implications for delivery of 60 Fe and other radioisotopes to the Earth and Moon. We see that the closest approach of the blast recedes quite rapidly at first. If the supernova ejecta is wellmixed into and carried by the blast, this means that the distance it must travel through the heliosphere rapidly becomes larger over time. If instead the ejecta is not well-mixed but located well behind the forward shock, it has even farther to travel in the heliosphere than in the well-mixed case. This again points to the need for the radioisotopes to arrive on dust grains having large sizes or high speeds.
After the forward shock arrival, much of the solar system will be exposed to the blast wave. Fig. 11 also shows the region of the Kuiper belt from 30 to 55 au (Sanctis et al. 2001;Chiang et al. 2003). Even for the most distant 100 pc supernova, the entirety of the Kuiper belt is exposed for ∼10 kyr.
It is interesting to note that due to the slightly sharper weakening of the closer supernovae, the outer Kuiper belt at 55 au is exposed to the supernova blast for approximately 70 kyr regardless of distance. In contrast, the inner Kuiper belt at 30 au is more sensitive to the distance, resulting in a nearly linear relation with distance.

Other solar system effects
Given the extent to which the heliosphere can be compressed, parts of the outer solar system are directly exposed to the supernova blast. This exposure may have numerous effects on the outer bodies. Stern & Shull (1988) investigated how the light from a nearby supernova could melt the outermost surface of comets. Stern (1990) calculated the erosion of these bodies due to SNRs and how small particles with radius 100 µm would be ejected from the solar system. Both of these studies can now be re-contextualized in light of the known supernovae detected by terrestrial 60 Fe: the most recent 3 Myr supernova may have "cleaned" the Oort cloud of small dust grains but left larger comet orbits unperturbed.
Due to the increased abundance of cosmic rays, the cosmic ray exposure ages of asteroid surfaces may be affected. While typical isotopes of interested for exposure ages have half-lives of less than 1 Myr (see, e.g., Michel et al. 1997), some of the longer-lived radioisotopes like 60 Fe, 53 Mn, and 26 Al could have high abundances. This effect would be more apparent in metallic meteoroids, which are exposed for longer than stony meteoroids (Ammon et al. 2009).
Potential links between SNRs and planets have so far been almost wholly unexplored. A significant challenge is to suggest not only how a supernova affects planets, but also find what observable features could remain after millions of years. We leave these investigations to future studies.

Supernova effects on astrospheres
In these simulations, we examine the supernova blast wave's effect on our own heliosphere. These simulations can be generalized to examine the effect of supernovae on astrospheres, which are driven by stellar winds from other stars.
Astrospheres are commonly observed by their bow shocks in the IR part of the spectrum. They are most commonly seen in surveys, either by the rapidly-moving runaway stars (Peri et al. 2012) or in the dusty environment of the Galactic plane (Kobulnicky et al. 2016). Red supergiants in particular, such as Betelgeuse, are expected to have enormous astrospheres nearly a parsec wide (Meyer et al. 2021). As of yet, no bow shocks have been associated with a star located within an SNR. Indeed, supernovae are prolific destroyers of ISM dust, so it is perhaps expected that astrospheres in SNRs would not be detectable in the IR due to a lack of dust. However, if discovered, they would be a novel form of stellarinterstellar interaction.
Astrospheres present an interesting complementary aspect to our own heliosphere: though we can directly probe our own heliosphere, the overall shape of the heliopause is still debated. It may have a comet-like tail 1000s of au long (e.g., Izmodenov & Alexashov 2015) or it may be truncated much closer (Opher et al. 2015). In contrast, known astrospheres cannot be probed directly, but their shape can be seen by their bow shocks. Both the solar system penetration distance and the increase in cosmic rays may alter the potential habitability of astrospheres. The proximity of stellar systems to supernovae affects the Galactic Habitable Zone (Gonzalez et al. 2001;Lineweaver et al. 2004;Morrison & Gowanlock 2015;Spinelli et al. 2021). With these simulations, more accurate calculations of cosmic ray exposure during supernova blasts can be made in order to better ascertain astrosphere viability.

Effects of solar motion
We have assumed the solar system is at rest relative to the supernova explosion. In general one expects the Sun will move relative to the supernova progenitor and blast center. A nonzero solar velocity relative to the blast would change the ram pressure seen by the heliosphere, and the impact of this change scales at δP ram /P ram ∼ v /v blast . The Sun's present motion with respect to the stellar Local Standard of Rest is 18 km/s (Schönrich et al. 2010;Zbinden & Saha 2019) and our speed relative to the very local ISM is 27 km/s; for such values v /v blast 1 when blast arrives, and the perturbation is small.
At late times, the blast speed slows and Earth's speed could become important; these effects are discussed in (Chaikin et al. 2021) in the context of 60 Fe deposition and the local environment encountered by the solar system. It remains for future work to model such effects on the heliosphere, including the possibility that the Earth's velocity is misaligned with that of the supernova blast.

CONCLUSIONS
Motivated by terrestrial detections of 60 Fe as evidence for near-Earth supernovae in the recent past, we have presented hydrodynamic simulations of the heliosphere's response to a supernova blast wave at various distances.
We match our steady solar wind to that observed by space missions and test the effect of solar wind orientation. The supernova blast is assumed to be in the Sedov phase.
The broad structure of the heliosphere is reproduced, albeit at much smaller scales than the present-day heliosphere. We verified that pressure balance gives the location of the TS over a range of distances. We applied this relation to analytically examine the heliosphere throughout the duration of the supernova remnant evolution far longer than the hydrodynamic simulations could run.
Taking advantage of rotational axisymmetry allowed for two orientations of how the blast wave strikes the heliosphere: polar, in which the blast approaches from the poles of the Sun, and equatorial, in which the blast arrives from the side. We apply a steady fast and slow solar wind originating from the poles and equator, respectively. Due to their density differences, the ram pressures of these winds are very similar. Appropriately, since the penetration distance depends chiefly on ram pressure, the orientation is found to have little influence on the global structure of the heliosphere.
This work reaffirms and builds upon the conclusions of Fields et al. (2008) that the supernova blast plasma is strongly excluded from 1 au for any plausible distance to the recent 60 Fe-depositing supernovae. The observed 60 Fe deposits on the Earth and Moon must have arrived in a form other than the plasma-namely, in dust grains. The dynamics of dust grains in the outer heliosphere are well-studied (e.g., Belyaev & Rafikov 2010), especially for the present-day heliosphere. Wallis (1987) found that during the passage through a dense cloud, dust can penetrate the heliosphere to Earth with little deflection due to the heliosphere's small size. Athanassiadou & Fields (2011);Fry et al. (2016) found that dust grains from near-Earth supernovae are typically deflected less than 1 • by the heliosphere. Our work shows provides the location of closest approach of the supernova material and thus the initial conditions for studies of the dust propagation within the compressed heliosphere to Earth.
As the SNR evolves, the blast wave will weaken and allow the heliosphere to rebound. According to our scaling laws, this process is expected to take several 100 kyr to rebound to 100 au. Our simulations do not account for supernova-formed dust dynamics, but their propagation through the heliosphere is modified by the decreased heliosphere size. This effect is especially significant for any grains that arrive within the first 100 kyr after the blast wave arrival. Grains that arrive later will need to traverse progressively more of the heliosphere in order to reach the Earth and Moon.
The principal shocks in our simulations (BS and TS) can accelerate particles through diffusive Fermi acceleration to produce anomalous cosmic rays (Zank et al. 1996;Lazarian & Opher 2009). Cosmic rays accelerated in these shocks would have an effect on the amount of radiation impinging on Earth. For very nearby supernovae, there could even be biological effects, either as a mass extinction (Gehrels et al. 2003) or lesser extinction events (Melott et al. 2017;Thomas et al. 2016). Tracking how these shocks evolve over the passage of the blast wave furthers our understanding of potential biological responses to the supernova.
The proposed Interstellar Probe mission 2 plans to launch a spacecraft several hundred au into the very local ISM in the coming decades. Such a probe would be sent out to the remains of an ancient supernova remnant, and may thus contribute to the study of how old SNRs evolve and fade into the Galactic medium.
The simulations presented here could be expanded upon in many ways, including developing a more careful treatment of non-hydrodynamic physics like magnetic fields and charge exchange. Solar activity represented as a time-varying solar wind may also be relevant for determining how instabilities affect the distance of closest approach. Such inclusions would allow for more realistic non-axisymmetric simulations, allowing us to probe how our own solar system responds to dramatic nearby events such as supernovae.