Heliospheric Compression Due to Recent Nearby Supernova Explosions

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

Indeed, some of the relevant physics for the modern heliosphere may not be applicable to 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 an 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 Figure 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,

$$\begin{eqnarray}&&\displaystyle \frac{\partial \rho }{\partial t}+{\boldsymbol{\nabla }}\cdot (\rho {\boldsymbol{v}})=0,\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial \rho {\boldsymbol{v}}}{\partial t}+{\boldsymbol{\nabla }}\cdot (\rho {\boldsymbol{vv}}+p)=0,\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial E}{\partial t}+{\boldsymbol{\nabla }}\cdot [(E+p){\boldsymbol{v}}]=0,\end{eqnarray} \tag{ 3 }$$

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.

The solar wind is the complex flow of gas launched from the solar corona and streaming outward 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 these chosen values every timestep.

The real solar wind not only varies in time, but with latitude: it is launched differently in the plane of the solar system than toward the poles. We use data from Ulysses's fast latitude scan 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 are 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 the 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. Voyager 2 was launched within the equatorial region of the solar wind, but has steadily drifted toward southern latitudes. It passed 25° at approximately 75 au (see, e.g., Hill et al. 2020), shortly before the termination shock (TS). Therefore, most of the data we compare to is from the equatorial solar wind.

Zoom In Zoom Out Reset image size

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, as emphasized in a recent comparison with New Horizons data (Elliott et al. 2019). 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 measured 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's fast scan. 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.

Zoom In Zoom Out Reset image size

Currently, it is an open question as to where the supernova 3 Myr ago exploded, and therefore which direction the blast wave would have 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 outward 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, 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 the angle described above; this is the polar orientation. When the blast comes in along the plane (striking the outer planets' orbits first) as shown 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.

Zoom In Zoom Out Reset image size

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 relaxed solar wind and corresponding meshblocks in the polar orientation. Each refinement level increases the resolution by a factor of 2.

Zoom In Zoom Out Reset image size

SNRs evolve over time. After the initial free expansion phase, SNR 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 well approximated as a plane wave on the scale of the solar system. Since the most important parameter in our simulations is the distance to the supernova, ${R}_{\mathrm{SN}}$, we invert the usual Sedov equation for distance to solve for time as

$$\begin{eqnarray}&&t=\sqrt{\displaystyle \frac{{R}_{\mathrm{SN}}^{5}{\rho }_{0}}{{\beta }^{5}{E}_{\mathrm{SN}}}}=27\ \mathrm{kyr}\ {\left(\displaystyle \frac{{n}_{0}}{0.01\ {\mathrm{cm}}^{-3}}\right)}^{1/2}{\left(\displaystyle \frac{{R}_{\mathrm{SN}}}{50\ \mathrm{pc}}\right)}^{5/2},\end{eqnarray} \tag{ 4 }$$

where ${R}_{\mathrm{SN}}$ is the radius of the remnant, ρ₀ is the ambient medium density, β is a numerical factor of 1.1517 for the typical γ = 5/3, and ${E}_{\mathrm{SN}}$ is the energy of the supernova, taken to be 10⁵¹ erg. By applying the Rankine–Hugoniot conditions, we find the density, velocity, and thermal pressure of the gas immediately post-shock, which are

$$\begin{eqnarray}&&{\rho }_{1}={\rho }_{0}\displaystyle \frac{\gamma +1}{\gamma -1},\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&{v}_{1}=\displaystyle \frac{2}{\gamma +1}{v}_{s},\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}&&{p}_{1}=\displaystyle \frac{2}{\gamma +1}{\rho }_{0}{v}_{s}^{2},\end{eqnarray} \tag{ 7 }$$

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

$$\begin{eqnarray}&&{v}_{s}=\displaystyle \frac{2}{5}\displaystyle \frac{{R}_{\mathrm{SN}}}{t}=\displaystyle \frac{2}{5}\sqrt{\displaystyle \frac{{E}_{\mathrm{SN}}}{{\beta }^{5}{\rho }_{0}}}{R}_{\mathrm{SN}}^{-3/2}\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&=\,730\ \mathrm{km}\,{{\rm{s}}}^{-1}\ {\left(\displaystyle \frac{n}{0.01\ {\mathrm{cm}}^{3}}\right)}^{-1/2}{\left(\displaystyle \frac{{R}_{\mathrm{SN}}}{50\ \mathrm{pc}}\right)}^{-3/2}.\end{eqnarray} \tag{ 9 }$$

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⁻³, 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 & Breitschwerdt (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

$$\begin{eqnarray}&&T=\displaystyle \frac{2{\beta }^{5}}{25{k}_{B}(\gamma +1)}\displaystyle \frac{{E}_{\mathrm{SN}}}{{r}_{s}^{3}{n}_{1}},\end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray}&&T=1.5\times {10}^{6}\ {\left(\displaystyle \frac{{R}_{\mathrm{SN}}}{100\ \mathrm{pc}}\right)}^{-3}{\left(\displaystyle \frac{{n}_{0}}{0.01\ {\mathrm{cm}}^{-3}}\right)}^{-1}\ {\rm{K}}.\end{eqnarray} \tag{ 11 }$$

For the low density observed in the Local Bubble, this model predicts that the temperature post-shock will remain over 10⁶ K for 100 pc, the approximate size of the Local Bubble. In the past, it was accepted that the Local Bubble consisted of a large hot component with T ∼ 10⁶ K, though such assumptions have recently been challenged (Welsh & Shelton 2009; Linsky & Redfield 2021). Despite the hot environment of SNRs, 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}_{\mathrm{ram}}\sim {E}_{\mathrm{SN}}/{R}_{\mathrm{SN}}^{3}$. The supernova 3 Myr ago is estimated to have occurred somewhere between 60 and 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 modeling the strongest part of the blast, the forward shock. See Section 4.2 for a 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.

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 (Kobulnicky et al. 2016; Henney & Arthur 2019; Meyer et al. 2021). These studies show that the solar wind-ISM interaction region consists of three main features: the 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 (McComas 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

$$\begin{eqnarray}&&{r}_{\mathrm{stag}}=\sqrt{\displaystyle \frac{\dot{M}{v}_{w}}{4\pi {\rho }_{0}{v}_{* }^{2}}},\end{eqnarray} \tag{ 12 }$$

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 SNR. This stagnation distance is

$$\begin{eqnarray}&&{r}_{\mathrm{stag}}=\sqrt{\displaystyle \frac{{P}_{\mathrm{sw}}+{\rho }_{\mathrm{sw}}{v}_{\mathrm{sw}}^{2}}{{P}_{\mathrm{SNR}}+{\rho }_{\mathrm{SNR}}{v}_{\mathrm{SNR}}^{2}}}\ \mathrm{au},\end{eqnarray} \tag{ 13 }$$

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

$$\begin{eqnarray}&&{r}_{\mathrm{stag}}={\left(\displaystyle \frac{25(\gamma -1)}{8{\beta }^{5}}\displaystyle \frac{{R}_{\mathrm{SN}}^{3}}{{E}_{\mathrm{SN}}}\left({P}_{\mathrm{sw}}+{\rho }_{\mathrm{sw}}{v}_{\mathrm{sw}}^{2}\right)\right)}^{1/2}\end{eqnarray} \tag{ 14 }$$

$$\begin{eqnarray}&&=A{\left(\displaystyle \frac{{R}_{\mathrm{SN}}}{100\ \mathrm{pc}}\right)}^{3/2}{\left(\displaystyle \frac{{E}_{\mathrm{SN}}}{{10}^{51}\ \mathrm{erg}}\right)}^{-1/2},\end{eqnarray} \tag{ 15 }$$

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.

In model 1, we use the same input parameters as those in 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.

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 TD 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 high-resolution 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.

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 5(a), 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 5(a), the BS has receded out of the simulation frame, so it is only marked by a small uptick at the domain boundary.

Zoom In Zoom Out Reset image size

The HP is located where the two flows meet, which should be near-zero velocity along this axis. Models 5(a) and 11 clearly show this sharp downward 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 in 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.

Zoom In Zoom Out Reset image size

When the TS meets itself at a point in the downwind side of the simulation, it splits into three features as labeled in Baranov & Malama (1993). These are (from the inside outward) the MD, the secondary TD, and the RS (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 find evidence of these features.

For our Sedov supernova blast, we use a uniform ISM density of n = 0.005 cm⁻³. 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 Equation (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 5(a)–(c). The ambient densities used here are 3.34 × 10⁻²⁶, 1.34 × 10⁻²⁵, and 6.68 × 10⁻²⁵ g cm⁻³ (i.e., ρ₀, 4ρ₀, and 16ρ₀) for models 5(a)–(c), respectively.

We show these simulations after the same amount of time has passed since the introduction of the blast shown 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 Figure 8(c), we see a clear bow shock and lack of features in the downstream TS, 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.

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

Table 2. Supernova-Heliosphere Collision Simulation Results

Label	# of SMR	Orientation a	${R}_{\mathrm{SN}}$	ρSNR	vSNR	pSNR	rstag	rTS	rHP	rBS
	Levels		(pc)	(g cm−3)	(km s−1)	(erg cm−3)	(au)	(au)	(au)	(au)
1 b	0	E	20.0	6.40e-25	688	1.00e-9	2.34	2.47	3.37	8.90
2	1	E	25.3	3.34e-26	2141	5.11e-10	3.18	3.43	4.75	14.47
3	1	P	25.3	3.34e-26	2141	5.11e-10	2.99	3.08	3.60	12.25
4	1	E	50.0	3.34e-26	771	1.99e-10	8.83	9.72	15.97	⋯
5(a)	1	E	63.3	3.34e-26	541	3.26e-11	12.56	13.72	23.88	⋯
5(b) c	1	E	63.3	1.34e-25	271	3.26e-11	12.56	13.82	23.97	⋯
5(c) c	1	E	63.3	6.68e-25	121	3.26e-11	12.56	13.82	22.61	38.18
6	1	P	63.3	3.34e-26	542	3.26e-11	11.82	12.94	18.90	39.94
7	1	E	75.9	3.34e-26	412	1.89e-11	16.51	17.53	29.98	⋯
8	2	P	94.9	3.34e-26	295	9.68e-12	21.72	22.31	23.10	64.26
9	2	E	110.0	3.34e-26	236	6.22e-12	28.80	30.76	44.63	86.33
10	2	P	110.0	3.34e-26	236	6.22e-12	27.11	26.07	27.05	80.08
11	2	E	126.5	3.34e-26	192	4.09e-12	35.53	37.21	52.15	97.27

Notes.

^a E = equatorial, P = polar orientation. Solar wind parameters for these are given in Table 1. ^b Fields et al. (2008) model 12 comparison, different solar wind parameters. ^c Models 5(b) and 5(c) use ISM densities of n = 0.02 and 0.1 cm⁻³, respectively.

Download table as:  ASCIITypeset image

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.

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.

Zoom In Zoom Out Reset image size

Similarly to Figure 9(a), we can plot the TS and HP distances as a function of ${R}_{\mathrm{SN}}$, as shown in Figure 9(b). We also plot Equation (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}_{\mathrm{stag}}\propto {R}_{\mathrm{SN}}^{3/2}$ 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.

Zoom In Zoom Out Reset image size

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.

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 of 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 timescales. 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 ⁴ 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⁻³. At each timestep, we obtain the thermal and ram pressure and use Equation (13) to calculate r_stag (indicating the TS position). We show the results of these calculations in Figure 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 (Equation (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 ∼few kyr, however, it takes >100 kyr for the blast to weaken enough so that the heliosphere fully rebounds.

Zoom In Zoom Out Reset image size

Figure 11 has important implications for the delivery of ⁶⁰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 well mixed 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. Figure 11 also shows the region of the Kuiper Belt from 30–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.

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 ⁶⁰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 ⁶⁰Fe, ⁵³Mn, and ²⁶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.

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 stellar-interstellar interaction.

Astrospheres present an interesting complementary aspect to our own heliosphere: though we can directly probe our own heliosphere, the overall shape of the HP is still debated. It may have a comet-like tail thousands of astronomical units 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.

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 the 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 ⁶⁰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 Earth's velocity is misaligned with that of the supernova blast.