Oblique Shock Breakout in SNe

Oblique Shock Breakout in Supernovae and Gamma-Ray Bursts: II. Numerical Solutions For Non-Relativistic Pattern Speeds

Pegah Salbi,11affiliation: Department of Astronomy and Astrophysics, University of Toronto, 50 St. George Street, Toronto, Ontario, M5S 3H4, Canada Christopher D. Matzner,11affiliation: Department of Astronomy and Astrophysics, University of Toronto, 50 St. George Street, Toronto, Ontario, M5S 3H4, Canada Stephen Ro,11affiliation: Department of Astronomy and Astrophysics, University of Toronto, 50 St. George Street, Toronto, Ontario, M5S 3H4, Canada and Yuri Levin22affiliation: Monash Center for Astrophysics, Monash University, Clayton, VIC 3800, Australia salbi@astro.utoronto.ca

Non-spherical explosions develop non-radial flows as the pattern of shock emergence progresses across the stellar surface. In supernovae these flows can limit ejecta speeds, stifle shock breakout emission, and cause collisions outside the star. Similar phenomena occur in stellar and planetary collisions, tidal disruption events, accretion-induced collapses, and propagating detonations. We present two-dimensional, nested-grid Athena simulations of non-radial shock emergence in a frame comoving with the breakout pattern, focusing on the adiabatic, non-relativistic limit in a plane stratified envelope. We set boundary conditions using a known self-similar solution and explore the role of box size and resolution on the result. The shock front curves toward the stellar surface, and exhibits a kink from which weak discontinuities originate. Flow around the point of shock emergence is neither perfectly steady nor self-similar. Waves and vortices, which are not predominantly due to grid effects, emanate from this region. The post-shock flow is deflected along the stellar surface, and its pressure disturbs the stellar atmosphere upstream of the emerging shock. We use the numerical results and their analytical limits to predict the effects of radiation transfer and gravity, which are not included in our simulations.

Subject headings:
gamma rays: bursts – hydrodynamics – shock waves – (stars): supernovae: general
slugcomment: Submitted to ApJ

1. Introduction

The arrival of a normal shock at the surface of an exploding star is associated with a flash of radiation and the release of the fastest stellar ejecta, which can then interact with circumstellar material in the earliest phase of a supernova remnant. All of these phenomena result from shock acceleration in the steeply declining density profile of the outer stellar envelope: a whip-like motion described by the similarity solution of Gandel’Man & Frank-Kamenetskii (1956), Sakurai (1960), and Matzner & McKee (1999). The outcome is very different, however, when the explosion shock is not aligned with the density gradient, so that non-radial flows develop (Matzner et al. 2013, hereafter Paper 1). In this case, the pattern of shock emergence moves across the star’s surface at some speed . A natural length scale is the ‘obliquity’ depth at which the shock front is expected to move outward at . So long as radiation is trapped at this depth, motions become strongly non-radial. The shock and ejecta velocity are both limited (to and , respectively, in the star’s rest frame) and the ejecta spray is expected to suppress the escape of photons which would otherwise contribute to the photon flash. In the aftermath, it is possible that non-radial ejecta will collide outside the star.

Fully characterizing this behavior requires numerical simulations, as the analytical arguments presented in Paper 1 leave many questions unresolved. The detailed shape of the shock front, the shape and terminal angle of the outflowing stream lines, and the motions of matter and energy in the emerging flow are all to be determined. Moreover, Paper 1 showed that the flow immediately around the point of breakout cannot be both steady and self-similar in the comoving frame: it could be oscillatory, or the outflow could interact with the star to spoil the apparent self-similarity. We construct numerical experiments to address these questions and glean additional insights into the nature of oblique breakouts.

2. Physical problem and numerical implementation

We simulate only a single, asymptotic version of the oblique breakout problem, which was also a focus of attention in Paper 1. We ignore the diffusion of radiation, which is valid if is below the depth at which diffusion becomes important. We also ignore the curvature of the stellar atmosphere and any curvature of the shock front: a necessary condition for this to be valid is that is much less than the stellar radius. We can thus consider the problem in a frame of reference co-moving with the breakout pattern, and we can focus our attention on the two-dimensional, stationary flow of an adiabatic gas, separated by a shock front from plane-stratified upstream fluid. Furthermore, we specialize to the case where the post-shock flow is radiation dominated (adiabatic index ) and consider a cold polytropic atmosphere (index or ) for the pre-shock matter. Finally, we assume that the motions are purely non-relativistic, and that stellar gravity and hydrostatic pressure are both negligible. This combined limit is most applicable to aspherical explosions of compact stars, as was demonstrated in section 5 of Paper 1, but it remains a useful analogue in all cases where obliquity affects the flow. We shall return to consider the physics we have limited in later sections (radiation transfer in §5; gravity in §7; relativity and finite depth in Table 1).

These choices have several advantages. First, as discussed in Paper 1, the shock front is effectively horizontal far below a depth of . Because we ignore stellar curvature, all the streamlines which originate from this zone must approach the known planar, self-silmilar solution. We use this fact when determining both our outer boundary conditions and our initial conditions, as we describe more fully below. The pattern speed, obliquity scale, and density at the obliquity scale set our code units: .

Second, by ignoring the length scales associated with curvature of the star or the shock front, as well as the length scale on which radiation diffusion becomes important, we are left with only one physical length scale, , against which to judge our numerical parameters. If we conduct our simulation in a square box of width whose square grid has a finest scale , its outcome can be characterized by the box-size parameter

and the resolution parameter

Third, we expect that once the outcome is normalized to its natural units of length, time, and density (, , and , respectively), it is uniquely determined (at least in a time-averaged sense) by the post-shock adiabatic index and the polytropic index . (We address only the case in our simulations.) So long as our boundary conditions are consistent with this flow in the limit of an infinitely large box, this expectation gives us confidence that the results (apart from an initial transient) must limit to the definite physical solution as .

2.1. Code and code test

We employ the Athena magneto-hydrodynamics code (Stone et al., 2008) with magnetic fields turned off. Although this code is thoroughly tested on other problems, we wished to check its performance in the context of an accelerating shock. We therefore set up a one-dimensional problem in which a shock front moves through cold gas away from a high-pressure region (a region of high initial temperature, next to a reflecting boundary condition), descending an initial density distribution (plus a floor density times lower than the minimum for ), where, as in Paper 1, the coordinate is an altitude. The adiabatic index is . The results of this test are presented in Figure 1, where we plot the time to shock emergence against the initial depth. The results adhere within 0.4% to the theoretical power law derived by Ro & Matzner (2013) using the self-similar theory of Sakurai (1960).

Figure 1.— Shock velocity as a function of fluid density in a one-dimensional test problem (§2.1). A stellar atmosphere of depth two arbitrary units is resolved by grid cells, of which the outer 75 are plotted as dots, and a shock is launched from the innermost region (right-hand side in this figure). Plotted is the time to shock emergence, calculated using the time of greatest compression in each cell, against the distance to the surface. Finite resolution affects shock propagation within about ten cells from the surface, and our initialization affects it for the deepest matter, but for depths between and 0.7 in these units, the shock acceleration law agrees with the self-similar theory of Sakurai (1960; dashed line) within 0.4%.

For our two-dimensional simulations Athena was compiled with Message Passing Interface (MPI) and Static Mesh Refinement (SMR) enabled. We used a Harten-Lax-van Leer-Contact (HLLC) Riemann solver with second-order reconstruction, with a Corner Transport Upwind (CTU) unsplit integrator, and Courant number of 0.8. The computations were performed on 256 cores of the GPC cluster at the SciNet facility located at the University of Toronto (Loken et al., 2010).

We employ a sequence of three nested grids, in which one grid (of resolution up to ) is nested within an identical grid scaled up by a factor of two in size, and this in turn is nested within a third identical grid, another factor of two larger, for an effective linear dynamical range . (A code error prevented us from using even more than three resolution levels.)

2.2. Initial and boundary conditions

We appeal to the known planar, self-similar solution, and the fact that the flow limits to this solution for large initial depths () to set our initial and boundary conditions. In this planar solution, profiles of velocity, density, and pressure behave self-similarly, adhering to fixed functional forms which scale in amplitude and length scale as the shock approaches the stellar surface. To use it, we translate this time-dependent, one-dimensional flow into two dimensions. This involves several steps. First, we match the upstream density in the one-dimensional solution to the depth-dependent initial density in the 2D grid. Then we match the shock strength, ensuring where ; this enforces the definition of from Paper 1. Third, we associate each time in the 1D solution with a location in the 2D flow: where , and copy all of the fluid variables from the 1D solution to the 2D simulation volume. (This step is performed using a seventh-order polynomial fit to the self-similar solution, accurate to a few tenths of a percent in each variable, for computational ease.) Matter ahead of the shock front is assigned its initial density, for ). For we set the density equal to a low value (a lower density than for any ), and pressure to a floor value of in code units. In other words, we ignore the portion of the 1D solution in which ejecta flies away (, ); this material does not contribute to the final solution, and introduces very strong gradients which lead to numerical problems. The pressure floor is applied also to the pre-shock fluid (). Finally, we set the horizontal velocity: everywhere: this enforces the inflow and outflow of matter across the left and right grid boundaries, respectively. These initial conditions are depicted in the first panel of Figure 2.

Except for our treatment of the material, these ‘self-similar’ initial conditions correspond to a scenario in which the fluid responds only to the vertical component of its pressure gradient, as it does in the 1D solution. For matter which originates deep within the star () this approximation is valid, as its shock normal and the post-shock pressure gradients are indeed vertical. For matter which originates in the oblique zone () it is far from correct, and we expect a transient period of readjustment, lasting a few flow times , before it can establish a self-consistent stationary state. Because our runs are limited to finite volume (finite ), we cannot perfectly reproduce the ideal solution.

These initial conditions also provide boundary conditions for the remainder of the simulation. The bottom boundary, and the inflow region () of the left boundary are set to enforce inflow: the initial fluid variables are fixed (in a ghost region) for all time. All other boundaries are assigned outflow conditions, in which accelerations at the boundary are calculated using a linear extrapolation of the fluid variables. The final stationary solution has supersonic outflow in all the outflow regions of the grid boundary, so our choice is self-consistent.

In the project’s initial stages we used a small wedge of very high-density zones to launch an oblique shock upward through the density gradient. Although the flow very close to the stellar surface was similar in these runs to what we saw later with self-similar boundary conditions, we deemed this procedure unsatisfactory. A sequence of runs with different wedge parameters (height, angle, overdensity factor) revealed that a large portion of the simulation volume was affected by these parameter choices. Mixing of matter between the wedge and the stellar envelope, and the dynamical reaction of the wedge, complicated the analysis. Furthermore high density contrasts between the wedge and the envelope induced frequent crashes of the code. For these reasons we report only on runs with the self-similar boundary conditions described above, which have the benefit that they become exact in the limit .

3. Results - Fiducial simulation

Figure 2 depicts the progression of the density distribution in our fiducial run, for which , , the largest grid extends over , and the finest grid extends over , . Because this run represents the largest box and the highest resolution achieved in our study, we wish to describe its features before examining the effects of resolution (§4.1) and box size (§4.2) on the results.

As expected, we see a period of transient evolution from the initial state. Post-shock matter flies upward, establishing an outflow above the stellar surface, and is deflected upstream (to negative ) by its internal pressure gradients. The final panel of Figure 2 represents the statistically stationary final state, in which only short-period oscillations persist. Figures 3, 4, 5, 6, and 7 depict the shock structure, compression rate, specific vorticity, pressure, and entropy, respectively, of the final stationary state, and figure 8 provides a magnified view of the entropy structure in the highest-resolution subgrid, immediately around the breakout region. We do not plot velocity vectors, but we note that, in the shock’s frame, they are almost constant in magnitude (except in the immediate post-shock region where the shock is vertical) and they are parallel to contours of entropy (except where vortices develop). They are therefore very similar to the velocity field predicted in Figure 1 of Paper 1.

In the following subsections we comment on several aspects of the final stationary state.

(a) Initial conditions and refinement zones
(b) Transient evolution
(c) Final stationary state
Figure 2.— Snapshots of density distribution in the fiducial run. The simulation is performed in a frame comoving with the pattern of shock emergence as it traverses the stellar surface. Matter flows from the left and across the shock front before accelerating outward. Note: in all figures except Figure 8, only the coarsest grid level is rendered. The boxes in the top panel denote the location of the finer grids.

3.1. Shock structure

During the emergence of an aspherical supernova shock, part of the stellar surface has been shocked while another part has not. As a consequence of its acceleration, the explosion shock (which divides these regions) curves outward, with its normal vector becoming more non-radial, as it approaches the stellar surface. Several observations made in Paper 1 form our expectations on the shape of the shock front.

First, because we limit ourselves to adiabatic flow and neglect stellar curvature, the shock front will evolve self-similarly while it is normal, i.e., while it is much deeper than the obliquity scale. This fact implies a limiting form for the locus of points which define the shock surface, which indeed we have already imposed in our boundary conditions. In particular, our definitions imply


where for and , and our initial and boundary conditions correspond to the choice . (This choice differs from Paper 1, where we enforced at the point where the shock meets the surface. The two definitions are not equivalent, because the flow connects the self-similar offset to the point of shock emergence in a definite way.) This deep self-similar limit is plotted as a blue dashed line in Figure 3.

Second, the theory of Ishizuka et al. (1964) predicts a definite functional form for the shock front (up to a scaling factor) as its normal vector curves away from the vertical. A key feature of this prediction is that the shock angle cannot evolve continuously past a specific value (in the theory, this is when the shock tangent is 20.7 from the vertical). If the shock is to reach the surface, it must experience a kink, then travel to the surface at a set angle. However, this theory is only approximate. Its predictions are therefore tentative, as is clear from the fact that it is not consistent with the self-similar evolution in the deep limit (). Paper 1 hypothesized that the shock transitions smoothly from its deep limit to the form predicted by Ishizuka et al. at some reference depth.

Third, the similarity analysis presented in Paper 1 makes a very weak prediction that the shock should reach the surface vertically (with it a horizontal normal vector). This analysis provides an analytical form for the flow structure about the point of breakout, but it is not physical (for astrophysical values of and ) because it requires the sine of the angle between the shock normal and the vertical to exceed unity by a small amount (). While this primarily implies that the surface flow is either not steady or not self-similar (or both), it nevertheless suggests that, if the self-similar solution guides the flow in any way, its shock normal will be as close as possible to the self-similar condition. The closest physical solution is a vertical shock front ().

The compression rate in the fiducial run is depicted in Figure 4 at a late time () for which the flow has entered its final, stationary state; the color scale is chosen to highlight shocks while also making visible large-amplitude sound waves. The form of the shock locus corresponds well to our expectations, as depicted in Figure 3. In particular, the shock curves toward the surface in a manner similar to the Ishizuka et al. theory, experiences a kink at finite depth, and is nearly vertical from there to the breakout point. We measure a terminal shock angle of approximately from the vertical, rather than , just below the kink.

Figure 3.— Primary shock locus (red line) at in the fiducial run, compared with the asymptotic power-law form used to set boundary conditions (eq. 1, blue dash-dot line) and with the approximate model of Ishizuka et al. (1964, green dash-dot line). The two free parameters of this model are set by matching the value and slope to eq. (1) at , a choice which optimizes its match to the simulation. The simulation shock experiences a kink at , whereas the model’s kink is at . As in Paper I, we extend the model vertically to the stellar surface from this kink location. Comparing the original stellar surface (black dashed line) and the lower boundary of the region disturbed by outflow (black solid line) illuminates the ejecta-envelope interaction discussed in §3.4. The purple dotted line depicts equation (5), which approximates this interaction. The comoving-frame and rest-frame deflection angles ( and , respectively) are depicted for clarity.

We see in Figure 4 that the kink in the primary shock radiates a weak shock downstream, and in Figures 5 and 7 we see that a vortex sheet and mild entropy discontinuity also emanate from this point, before traveling downstream at a terminal angle . All of these features are required at a kink in the shock (Landau & Lifshitz, 1959, § 102, fig. 83). A close examination of the highest-resolution subgrid, shown in Figure 8, reveals that Kelvin-Helmholz rolls develop along this vortex sheet.

Figure 4.— Rate of compression () in the final, stationary state of the fiducial run.
Figure 5.— Specific vorticity () in the final, stationary state of the fiducial run.

Another weak shock emerges from the point of breakout, which we take to be the topmost location at which the primary shock meets matter of zero entropy, ; the distribution of the entropy-related quantity is plotted in Figure 7.

Several additional shock features are apparent in Figure 4, including two very weak shocks emanating from the curving primary shock, and a third weak shock which emanates from the point at which the primary shock meets the lower grid boundary. We suspect that all of these are related to the effect of finite box size, and are caused by discrepancies between the self-similar boundary conditions and the dynamics of smooth two-dimensional flow. We examine this question again in §4.2, where we vary to assess the influence of such discrepancies. The secondary shocks are hardly visible in the pressure distribution (Figure 6), indicating they are much weaker than the primary shock.

Figure 6.— Pressure in the final, stationary state of the fiducial run. The points A, B, and C at which we trace the time histories of various quantities, are overplotted.
Figure 7.— The entropy-related quantity in the final, stationary state of the fiducial run; see Figure 8 for details around the breakout region.
Figure 8.— The highest-resolution subgrid: entropy in the final state () of the fiducial run. The compression of the stellar surface upstream of the main shock is easily visible, as is the Kelvin-Helmholz instability which develops along the vortex sheet which emerges from the shock kink. Note: the image resolution does not approach the native resolution of this subgrid.

3.2. Distribution of ejecta

Observational implications of oblique shock breakout depend critically on the distribution of fluid quantities in the ejecta: that is, the distribution of mass flux and entropy relative to angle in the far field. A related question is how the final angle of ejection (, if measured in the comoving frame) relates to the initial depth ().

To gauge these distributions we must assign an approximate to each streamline, as follows. First we obtain an averaged final stationary state by time-averaging the conserved quantities over the period in the fiducial run. The velocity vector is not perfectly aligned with , thanks to acceleration outside the simulation volume. However, we can correct for these extra accelerations by appealing to the behavior of the planar, self-similar flow: a mass element’s velocity , expressed in the star’s rest frame, is related to its final value and the sound speed by a definite relationship. To an accuracy better than 0.2%,


Although this strictly applies only to the planar limit which holds deep within the star, we apply the corresponding velocity correction at every point in our solutions. We do this first by evaluating the local stellar-frame Mach number (), then computing , then boosting the fluid velocity by this amount in the direction opposite to the local pressure gradient. We take the terminal angle to be the angle of the resulting velocity.

This procedure should be valid in the very deep flow, where it captures the behavior of the self-similar solution, and in the shallow, highly supersonic region, where the velocity correction tends to zero. Any error incurred on intermediate streamlines should decline for large values of the box size parameter , because simulations with larger capture regions of higher Mach number. Fortunately, the angular correction in our fiducial run is at most a few degrees, and our estimate for hardly varies along each streamline.

The distribution of final quantities is plotted against angle in Figure 9. We list fits to several angular ranges within the fiducial run in Table 1, where we compare to the limit of deep planar flow (derived in the Appendix). This comparison, indicated by the dashed curves in Figure 9, indicates that our numerical estimate for is not especially accurate in the deep flow.

Quantity Fit: all Analytical limit () Equation
Table 1Ejecta distribution: fiducial runaaQuantities measured in the time-averaged final steady state of the fiducial run at a distance from the breakout point.

As viewed in the shock’s frame, the stellar envelope moves azimuthally (horizontally) through the shock front and is deflected radially (vertically) by an angle (see Figure 3), the terminal value of which is . Viewed in the star’s frame, matter at rest is struck by an inclined shock and ejected from the star at a terminal angle from the radial (vertical) direction. These two angles are related by the fact that the terminal speed matches the inflow speed in the comoving frame (a consequence of energy conservation), i.e. . This implies


The outflow speed in the star’s rest frame is

Figure 9.— Distributions of flow quantities as functions of the estimated terminal angle in the time-averaged final state of our fiducial run. The quantity is the distance to the breakout point, so all three quantities plotted are constant along streamlines at large . All values are normalized to natural units: . Dashed curves represent the known planar, self-similar solution, extrapolated as a power law in from deep planar limit (), as derived in the Appendix. The feature at corresponds to the shear layer visible in Figure 5, and emanates from the kink in the main shock.

3.3. Unsteady behavior

Our final state is not steady, as it exhibits acoustic oscillations, entropy fluctuations, and vortices. Unsteadiness is evident in the time history of pressure, entropy, and specific vorticity at the locations we plot in Figure 10. Our box-size and resolution studies show that the dominant oscillation frequencies are independent of and only weakly dependent on , so we infer that they are representative of the physical solution. There is no sign of any high-pitched acoustic oscillations, those of periods , which would be expected from interactions between the flow and the discretized grid.

Figure 10.— Time dependence of entropy (), specific vorticity () and pressure () at the points A, B, and C labeled in Figure 6. In each case we plot the mean-subtracted quantity normalized to unit standard deviation over the time interval shown: std. For clarity we show only .

The oscillation periods in vortical, entropy, and pressure variations are similar, and range from a fraction of to several . Because the shock structures are quite steady, we interpret the entropy fluctuations as being due to vortices moving fluid across the time-averaged streamlines. Unlike vorticity and entropy, pressure oscillations are restricted to certain flow zones: in particular, they are confined to be downstream of the shock which emanates from the breakout point. A plausible source for all these oscillations is the unstable break-up of the shocked slab discussed below.

3.4. Interaction of ejecta with the stellar surface

Paper 1 raised the possibility that interaction between the ejecta and the stellar atmosphere might be the reason that flow around the breakout point cannot be described with a steady, self-similar solution. In this case the pressure of the ejecta, enhanced perhaps by shocks in the region where it interacts with the star, would drive a weak shock and a downward flow in the upstream (otherwise undisturbed) stellar envelope. Self-similarity is then destroyed, or at least dramatically changed, because the density no longer varies as the th power of distance from the breakout point.

This scenario plays out precisely as described within our simulations. Pressure in the surface-skimming ejecta () compresses the outermost regions of the stellar envelope, depressing its interface with the stellar surface from to in the fiducial run. This effect is most visible in Figure 3. The primary shock front shows a visible feature where it intersects this layer, and we strongly suspect that some of the oscillations and vortices in the downstream flow are the consequence of this layer’s evolution behind the shock. In Figure 4, the local compression rate shows finite-amplitude sound waves or weak shocks in those ejecta that skim the stellar surface, and these radiate as sound waves into adjacent streamlines.

Let us compare the depression of the stellar surface layer with a simple calculation based on the observed ejecta pressure, which in the fiducial run is . Except for a minor correction due to the acceleration of the shell, this matches the pressure behind the atmosphere shock (subscript ), which we idealize as moving vertically downward. The instantaneous speed of the envelope shock, , equals because of the motion of the breakout point relative to the stellar surface. Using , this implies . Setting and integrating from at (as in the upstream boundary of our simulations), we find


As depicted in Figure 3, this is an excellent approximation to the shape of the inward shock in the fiducial run (); it predicts this shock will be found at a depth of near the breakout point. Furthermore, given the compression factor of 7, it predicts that the effective breakout point will be at six-sevenths of this depth, or . This is very close to the measured value of we list in Table 2.

It is important to note that the depression of the stellar surface diverges, albeit very slowly, with distance from the breakout point. If we associate the initial distance with the stellar radius, equation (5) indicates that the breakout depth will be about . This raises the possibility that the flow on scales of will be qualitatively changed if relative to what we can simulate in a finite box. We consider this and other finite-volume effects in §4.2.

4. Numerical effects

To disentangle physical phenomena from those imposed by the numerical implementation, we independently vary both the physical resolution and the size of the simulation volume. The simulations we use for this are listed in Table 2. In the subsections below we concentrate on these parameters’ effects, some of which we have already mentioned.

Run range range Levels Breakout pointaaLocation, in the final steady state, where the main shock breaches undisturbed stellar matter; its because of the downward reaction of the upstream stellar surface to out-flowing ejecta. Kink altitudebbAltitude, in the final steady state, of a kink in the shock angle. Match altitudeccAltitude, in the final steady state, of the best-fit matching between the deep asymptotic solution and the Ishizuka et al. (1964) approximate theory, as in Figure 3.
Fiducial 683 24 -17 to 7 -7 to 17 80 3 -1.32 -6.48
R341 341 24 -17 to 7 -7 to 17 80 3 -1.30 -6.4
R171 171 24 -17 to 7 -7 to 17 73.6 3 -1.28 -6.26
B12 683 12 -8 to 4 -4 to 8 80 3 -1.23 -4.60
Table 2Simulations used to explore resolution and box size effects

4.1. Finite resolution

Our resolution study consists of simulations with identical grid geometries and identical , in which varies between 683 and 171, corresponding to individual grids (at each resolution level) which vary between and . Pressure, entropy, and specific vorticity distributions for the lowest and highest resolution runs are depicted in Figure 12. Figure 11 shows the effects of on the oscillation modes.

Figure 11.— Time dependence of specific vorticity (), entropy (), and pressure () at point B (Figure 6) for runs with different numerical parameters. In each plot the fiducial run is a solid blue curve; the quarter-resolution run R171 is a black dashed curve; and the half box-size run B12 is a red dotted line. For clarity we show only .

The structure of the shock, the shapes of the post-shock streamlines, and the final entropy distribution are all very insensitive to . Oscillation frequencies are not strongly -dependent, except that the lower-resolution runs show a suppression of some of the short-wavelength motions. It appears that there is very little difference between the highest-resolution runs; suggesting that these faithfully represent the continuum limit . However it is always possible for very slow dependence, or phenomena which depend on threshold Reynolds numbers, to spoil this fidelity.

We have not conducted any runs of sufficiently low resolution, , such that the region of non-radial flow is poorly resolved. However we speculate in §8 about what our results might imply for global simulations which lack resolution of the surface layers.

4.2. Finite volume

Our simulation volume cannot approximate the ideal infinite case, especially as we were unable to use more than three levels of refinement. What makes this limitation especially severe for our physical problem is the fact that our boundary conditions are derived from the planar limit, which is only valid for matter which originates deep within the star (streamlines with ). Because the shock locus curves rapidly upstream (), it is difficult to reach these deep streamlines in a simulation of finite volume which must also capture the obliquity zone and the outflow region. We therefore expect to see -dependent features associated with a readjustment from the flow we impose at the boundaries toward something more representative of the ideal solution.

Prime candidates for these features are the weak shocks which emanate from the lower grid boundary: one at , and another where it meets the primary shock. A couple weak shocks are also launched from the curving primary shock, where there is no clear interaction with the stellar surface, and these may also depend on the simulation volume. As for the geometrical shape of the primary shock, the approximate theory of Ishizuka et al. (1964) would suggest that this is quite insensitive, in a power-law atmosphere, to the details of the simulation – at least, up to an overall scaling (the value of ) and the absolute location of shock breakout ( for ).

All of these expectations are validated in the first and third rows of Figure 12, which compare the pressure in simulations of different and identical . Their primary shock structure is very similar, except that going from the lower to the higher-resolution run, the depth of the shock kink is greater by and the point of breakout is shifted to the left by . The shapes of the streamlines are very similar, and they are almost identical once we apply a translation and rescaling to bring the breakout point and shock kink together. Those features which genuinely differ between the two runs are those which originate from the lower boundary. The weak shocks radiated by the primary shock also depend weakly on the distance to the box boundary.

However, as we note in §3.4, there is another feature of the flow which clearly depends on its finite : the ejecta-envelope interaction upstream of the breakout region. Our calculation there suggests the breakout point descends as in runs with progressively larger .

5. Radiation Diffusion and Appearance

RSGbbModel s15s7b2 of Woosley & Weaver (1995), provided by Stan Woosley. 14 490 1.20 0.566 0.520
BSGccShigeyama & Nomoto (1990) model for the progenitor of SN 1987A, provided by Ken’ichi Nomoto. 15 50 3.88 0.0254 0.409
Icdd Model CO6 of Woosley et al. (1999), for the progenitor of SN 1998bw, provided by Stan Woosley. Coefficients listed here are modified slightly by relativistic effects for energies .
5 0.2 5.12 0.350 0.435
Table 3Characteristics of model core-collapse supernovaeaaPolytropic parameters and are fit to hydrostatic regions in the outer 20% of the stellar radius. Shock parameters and are derived from fits by Ro & Matzner (2013); coefficient is adjusted to the stellar profile as described by Tan et al. (2001). The Tan et al. shock velocity model is used to calculate . The symbol indicates the ratio ; in the bipolar explosion model of Paper 1, where is an elongation factor.

In the diffusion approximation the radiative flux defines a diffusivity and a photon diffusion speed , where is the local scale length of radiation pressure. In a radiation pressure-dominated flow applies to the diffusion of pressure. In the context of a breakout flow, therefore, the dynamical effects of radiation diffusion are determined by the local Péclet number : diffusion is negligible where and strong where . In the context of real explosions, what does this mean for the validity of our adiabatic simulations and for observations of oblique breakouts?

Within our flow becomes constant with distance along each outflow streamline, because , while and . Fitting the angular dependence over all angles within our fiducial run we find that with and (see Table 1 and the Appendix for more information) so that


For comparison, an extrapolation from the deep, planar flow would give and .

In the numerical fit, the prefactor takes a minimum value of 0.13 in the turbulent boundary layer between ejecta and star (), where fluctuates. Inspecting this layer within the numerical solution, we see that dips to a minimum of where . Purely adiabatic calculations are appropriate (and breakout emission is entirely suppressed) where , so our neglect of diffusion along the stellar surface requires . If this condition is not satisfied, the deep radial flow will be unaffected but radiation diffusion will limit the hydrodynamic deflection to those streamlines that have . (To compare to the parameter used in Paper 1, note that . For the case, this implies .)

A similar criterion, discussed in Paper 1, compares the characteristic flow speed to a global radial diffusion speed , where is the optical depth to infinity at angle from the breakout point. If along a streamline, we expect its radiation to be trapped, but if then it may escape to an external observer. (Transfer across streamlines eases this criterion somewhat, by allowing photons from regions with to diffuse into streamlines on which ; however this is only likely where .) The integral diverges logarithmically at large in two-dimensional constant-velocity flow, and must be truncated at because three-dimensional effects set in on the scale of the stellar radius. We conduct the radial integration outward from the breakout point using our numerical results, and then extrapolate to assuming in each direction. If we assume the resulting is virtually identical to , so


and because the angle at which is only 1.2 times smaller than the angle at which (provided there is such an angle). Practically speaking, this means that the condition for photon diffusion to affect the flow () is only slightly less restrictive than the condition for photons to escape the system entirely ().

To evaluate the importance of diffusion in real stars, let us start by assuming (as we have in the simulations) that , so that the zone of oblique flow is in a thin outer layer where . The shock velocity in this zone, neglecting non-radial motions, is , if and are parameters in the theory developed by Matzner & McKee (1999), which can both be refined to account for details of the stellar structure (Tan et al., 2001; Ro & Matzner, 2013, see also eqs. (A6) and (A7)). The characteristic shock velocity is independent of location in a spherical explosion and nearly constant in the cases we consider here. Identifying the depth where ,


therefore, with equation (6),


where is a characteristic optical depth of the envelope, and is the polytropic index.

Our simulations only provide information about the case , but we can extrapolate to other polytropic indices by appealing to the fact that our numerical fit to is quite close to what we would have obtained by extrapolating from the deep, planar flow. In this deep-flow limit (), we find in the Appendix that


and , , and are all functions of . It is therefore reasonable to scale the numerical fits for and by the factors and , respectively. Although we do not yet know the accuracy of this extrapolation for significantly different , we adopt this approach in Table 3.

In Table 3 we use the properties of the model core-collapse explosions introduced in Paper 1 to identify , the critical angle at which , as well as the depth ratio and the relativity factor , in terms . We use this particular combination because, in the toy model of an aspherical explosion introduced by Paper 1, is a known function of latitude if the explosion is elongated by the factor at breakout: .

The last three columns of Table 1 allow us to estimate the range of over which the assumptions of our numerical models hold: initial plane symmetry (), non-relativistic flow , and adiabatic flow (). In all cases the lower limit, comes the requirement of plane symmetry. In extended stars, diffusion sets the upper limit ( and in red and blue supergiants, respectively), but in the compact Ic progenitor the upper limit is set by relativity ( and for and , respectively), because shock propagation continues to relativistic speeds in such stars. While our results are barely applicable to explosions in diffuse stars like red supergiants (for which diffusion sets in at a significant depth), they are relevant for more compact progenitors so long as the explosion is not too relativistic.

Finally, we note that the maximum post-shock pressure is reached at the location of the kink, where the shock velocity is and therefore . This corresponds to a maximum blackbody temperature in local thermodynamic equilibrium:  K for our model {RSG, BSG, Ic} progenitors. Observed photons may be significantly more energetic if their population falls short of LTE (Katz et al., 2010), or they may be nearly an order of magnitude lower if they are degraded adiabatically (to the pressures found along the stellar surface) before escaping. We leave this and other questions, such as the possibility of bulk Comptonization in the zone of ejecta-star shear flow, for later investigation.

6. Energy in non-radial motions

Some of the energy diverted by the oblique shock into non-radial flow will be available to power transients from collisions outside the stellar surface. To calculate this, we note that the component of velocity parallel to the stellar surface in the star’s frame (, in our terminology) equals at large distances. The non-radial component of the energy flux is therefore . Computing the non-radial energy per unit stellar surface area for all streamlines diverted more than in the shock’s frame (or in the star’s frame), we find


This integral diverges in the limit , because whereas . However, the energy in significantly deflected motions – deflected by at least in the star’s frame () – is about per unit area. This tends to strongly weight those regions of the stellar surface where the pattern speed is lowest, because . Integrating over the surface of the bipolar explosion model introduced in Paper 1, we find that a fraction of the explosion energy is channeled into such motions in our {BSG, Ic} model progenitors, so long as the conditions for planar, non-relativistic flow are all met.

7. The effect of gravity

In many of the astrophysical circumstances where oblique shock breakouts should occur, gravity is not guaranteed to be a negligible perturbation to the dynamics. Examples include the ejection of planetary atmospheres during planetary collisions and giant impacts (Genda & Abe, 2003); stellar collisions (Hawley et al., 2012); eruptions of luminous blue variable stars, such as Carinae (Smith, 2013), accretion-induced collapses of white dwarfs (Fryer et al., 1999; Tan et al., 2001), detonations on neutron stars (Zingale et al., 2001; Townsley et al., 2012), and compression shocks in stellar tidal disruption events (Guillochon et al., 2009).

To be specific, let us consider shock ejection from a spherical star with radius , escape speed , and surface gravity . The ratio between gravity and the characteristic acceleration of an oblique shock breakout is . Therefore, so long as the explosion is nearly spherical (), there exists a broad range of for which gravity is a small perturbation on the scales of our simulation (), but gravity is nevertheless important on scales of order because . Under these conditions we can extrapolate the influence of gravity from our current results, much as we did for radiation diffusion in § 5.

First, what is the condition for matter to be ejected? Matter be accelerated to in the rest frame of the star; since the maximum speed in that frame is , a necessary condition for escape is . However, while this condition guarantees that matter will escape for a brief period, it is probably not, in fact, sufficient for the ejection of matter in steady state. The reason is that any streamline along which matter is bound to the star must fall back, and its trajectory can intersect that of matter ejected at in the star’s frame. In particular, matter cast forward of the breakout region in the shock’s frame () may rain back on that region. Furthermore, since the mass flow per unit angle decreases rapidly with , the ram pressure of this returning matter may overwhelm that of the otherwise escaping streamlines. Matter with achieves a speed in the star’s frame, so we estimate the criterion for mass ejection to be


i.e., the patten speed must exceed the surface Kepler velocity .

This planar result is little more than an educated guess, however, for several reasons. The interaction between out-flowing and returning matter is guaranteed to be complicated. Spherical and non-steady effects will always be important for matter ejected at speeds of order . Moreover, itself could change if the explosion lifts enough material to change the gravitational potential (e.g. Wyman et al., 2004).

If condition (11) is not met, matter cast vertically at speed (in the shock’s frame) reaches a height before raining back. This and the nearby streamlines will return to collide with the matter emerging from the breakout region.

Second, how does gravity affect the deep flow when ? Equation (4) indicates that matter deflected by a sufficiently small amount, , moves below the escape velocity in the star’s frame. However, as spherical effects may be important at this depth, and as this matter is not in a region of strongly non-radial flow, the division between capture and escape is better described within a spherical theory (e.g., §2.6 of Tan et al. 2001).

8. Conclusions

We have presented high resolution, two-dimensional simulations of oblique supernova shock breakout, focusing on the limit of adiabatic, non-relativistic flow in a thin layer near the stellar surface in order to compare our results against the analytical predictions and suggestions of Paper 1. As expected, we find that the primary shock curves toward the stellar surface, experiences a kink, and reaches the surface almost vertically. The arrival of the shock at the surface is accompanied by a spray of matter ejected at a range of angles, including along the stellar surface. The stellar surface is compressed by ejecta pressure, so that the breakout point is inside the original stellar radius. The post-shock flow around the breakout point is neither steady nor self-similar: it exhibits acoustic oscillations, entropy fluctuations, and vortices. We study the effects of grid resolution and box size on the ejecta behavior, and show that the waves and vortices around the breakout point are independent of box size and only weakly dependent on resolution.

Our simulations assume radiation diffusion, relativistic effects, gravity, and stellar curvature are all negligible, yet provide some guidance for scenarios in which they are not. Diffusion, in particular, is more important for the shallowest (most strongly-deflected) material than for the deepest matter; we are therefore able to define a deflection angle above which diffusion is important and radiation may escape to be observed. Radiation often is completely trapped () in compact stars, as predicted by Paper 1. Gravitational effects are negligible for the strongly non-radial portions of the flow when the lateral pattern speed is well above the stellar escape velocity, but outflow is quenched by returning matter when is approximately the stellar Kepler velocity.

The vortices we observe in our strictly two-dimensional simulations are undoubtedly more pronounced than they would be three dimensions, thanks to the well-known inverse nature of 2D turbulence cascades (Kraichnan, 1967). As these vortices are superimposed on a rapidly expanding flow, we see little to suggest that three-dimensional results would be qualitatively different; however this must be checked with future simulations.

While the simulations presented here exploit the local, essentially two-dimensional nature of the problem to obtain very high effective resolution (up to 16,384), this is usually not possible in global simulations. In such simulations, the effect of limited resolution on the structure of non-radial flows will be critically dependent on our resolution parameter , i.e., the comparison between the local obliquity scale and the minimum grid spacing . If this ratio is of order unity or less, we anticipate that the maximum possible deflection is limited to the small value appropriate for matter at the minimum resolvable depth: in the star’s frame, this is an angle of approximately if the polytropic index is (equations [3] and [A1], taking and ). In other words, insufficient resolution should suppress the development of non-radial flows much as radiation diffusion does within extended stellar progenitors.

One might question whether simulations which assume perfectly trapped radiation can be relevant to observations, so we end by reviewing why they are. When radiation is truly trapped, as in aspherical type I supernovae, non-radial flows eliminate the shock breakout emission, which was expected in spherical theory, over much of the stellar surface. By deflecting ejecta, non-radial motions allow collisions outside the star which may give rise to a novel form of transient. For these reasons, and because the ejecta speeds are limited, circumstellar interactions and the early supernova light curve are altered. Our local simulations provide quantitative estimates for the output of each patch of the stellar surface, which can be integrated to give global estimates as we have done in §6.

When radiation is only partially trapped on the scales of interest, as in the explosions of blue supergiants, our results provide a way to estimate the critical angle above which radiation will escape ( in § 5), which we also expect to be the limiting deflection angle. We anticipate that the emerging luminosity will match the kinetic luminosity above in our adiabatic simulations, but radiation hydrodynamic simulations will be required to test this.

Finally, when radiation is poorly trapped (), as in the explosions of extended red supergiant stars, our simulations are not relevant: the dynamics of shock breakout are described, at least locally, by the theory for spherical explosions.

PS, CDM, and SR are supported by an NSERC Discovery Grant; YL’s research is supported by an ARC Future Fellowship. We thank Ian Parrish and Shane Davis for scientific discussions and extremely useful help with Athena, Kristen Menou for bringing to our attention the relevance of this theory to planetary collisions, and Nathan Smith and John Bally for pointing out the possible relevance to Car. We are also very grateful to James Guillochon for insightful comments and to Chris McKee and the referee for stimulating questions. Most of the computations were performed on the GPC supercomputer at the SciNet HPC Consortium. SciNet is funded by the Canada Foundation for Innovation under the auspices of Compute Canada; the Government of Ontario; Ontario Research Fund - Research Excellence; and the University of Toronto.
(a) Pressure , fiducial run
(b) Entropy , fiducial run
(c) Specific vorticity, fiducial run
(d) , quarter-resolution run R171
(e) , quarter-resolution run R171
(f) Specific vorticity, run R171
(g) , half-box-size run B12
(h) , half-box-size run B12
(i) Specific voriticity, run B12
Figure 12.— Comparison of the pressure, entropy, and specific vorticity between runs with different numerical parameters. All snapshots are taken at , and the region plotted is covered by all three runs. The effect of lower resolution can be seen in the spatial scale of oscillations in the middle panels, and the effect of a smaller box size can be seen in the breakout location and the appearance of a weak shock discontinuity on the lower panels.

Appendix A Deep flow, small-deflection limit of ejecta distribution

Within our idealization that is very much less than the stellar radius (and any other scale over which shock properties vary), the very deep flow () is governed by the planar, self-similar solution to shock acceleration and post-shock flow. As a result, we can use the solutions described by Sakurai (1960) and Matzner & McKee (1999) to determine the angular dependence of all the fluid variables which emerge from this deep region, i.e., the shallow-angle flow ().

In the deep flow where corrections due to non-radial motions are negligible, our definitions and the shock acceleration law imply . In planar flow each fluid element reaches a terminal vertical (radial) velocity ; given a total speed in the shock frame (a consequence of energy conservation), this implies a final angle in the small-angle limit.

In this limit each flow quantity will tend far downstream () toward a power-law form

where is the combination of , , and of the same dimensions as , and , , and are dimensionless constants. So, for instance, the relation can be expressed


Mass conservation requires , so


The post-shock entropy is , where and ; therefore


The pressure is so


To construct the pressure scale length consider that (where are unit vectors). The first term dominates for , so

The local diffusion parameter from § 5 is defined as . Because ,


Written out with , this yields the expressions given in equation (6) with the constants in equation (9).

Finally we wish to provide fits for and which, being specialized to the case , are simpler than, and comparably accurate to, those given by Ro & Matzner (2013): these are


where and ; and


where and . Equation (A6) has an r.m.s. absolute error of , and equation (A7) has an r.m.s. relative error of 0.35%, when compared against the numerical solutions obtained by Ro & Matzner.


  • Fryer et al. (1999) Fryer, C., Benz, W., Herant, M., & Colgate, S. A. 1999, ApJ, 516, 892
  • Gandel’Man & Frank-Kamenetskii (1956) Gandel’Man, G. M. & Frank-Kamenetskii, D. A. 1956, Soviet Physics Doklady, 1, 223
  • Genda & Abe (2003) Genda, H. & Abe, Y. 2003, Icarus, 164, 149
  • Guillochon et al. (2009) Guillochon, J., Ramirez-Ruiz, E., Rosswog, S., & Kasen, D. 2009, ApJ, 705, 844
  • Hawley et al. (2012) Hawley, W. P., Athanassiadou, T., & Timmes, F. X. 2012, ApJ, 759, 39
  • Ishizuka et al. (1964) Ishizuka, T., Hashimoto, Y., & Ōno, Y. 1964, Progress of Theoretical Physics, 32, 207
  • Katz et al. (2010) Katz, B., Budnik, R., & Waxman, E. 2010, ApJ, 716, 781
  • Kraichnan (1967) Kraichnan, R. H. 1967, Physics of Fluids, 10, 1417
  • Landau & Lifshitz (1959) Landau, L. D. & Lifshitz, E. M. 1959, Fluid mechanics
  • Loken et al. (2010) Loken, C., Gruner, D., Groer, L., Peltier, R., Bunn, N., Craig, M., Henriques, T., Dempsey, J., Yu, C.-H., Chen, J., Dursi, L. J., Chong, J., Northrup, S., Pinto, J., Knecht, N., & Van Zon, R. 2010, Journal of Physics Conference Series, 256, 012026
  • Matzner et al. (2013) Matzner, C. D., Levin, Y., & Ro, S. 2013, ApJ, 779, 60
  • Matzner & McKee (1999) Matzner, C. D. & McKee, C. F. 1999, ApJ, 510, 379
  • Ro & Matzner (2013) Ro, S. & Matzner, C. D. 2013, ApJ, 773, 79
  • Sakurai (1960) Sakurai, A. 1960, Comm. Pure Appl. Math, 13
  • Shigeyama & Nomoto (1990) Shigeyama, T. & Nomoto, K. 1990, ApJ, 360, 242
  • Smith (2013) Smith, N. 2013, MNRAS, 429, 2366
  • Stone et al. (2008) Stone, J. M., Gardiner, T. A., Teuben, P., Hawley, J. F., & Simon, J. B. 2008, ApJS, 178, 137
  • Tan et al. (2001) Tan, J. C., Matzner, C. D., & McKee, C. F. 2001, ApJ, 551, 946
  • Townsley et al. (2012) Townsley, D. M., Moore, K., & Bildsten, L. 2012, ApJ, 755, 4
  • Woosley et al. (1999) Woosley, S. E., Eastman, R. G., & Schmidt, B. P. 1999, ApJ, 516, 788
  • Woosley & Weaver (1995) Woosley, S. E. & Weaver, T. A. 1995, ApJS, 101, 181
  • Wyman et al. (2004) Wyman, M. C., Chernoff, D. F., & Wasserman, I. 2004, MNRAS, 354, 1053
  • Zingale et al. (2001) Zingale, M., Timmes, F. X., Fryxell, B., Lamb, D. Q., Olson, K., Calder, A. C., Dursi, L. J., Ricker, P., Rosner, R., MacNeice, P., & Tufo, H. M. 2001, ApJS, 133, 195
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description