Oblique Shock Breakout in Supernovae and GammaRay Bursts: II. Numerical Solutions For NonRelativistic Pattern Speeds
Abstract
Nonspherical explosions develop nonradial 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, accretioninduced collapses, and propagating detonations. We present twodimensional, nestedgrid Athena simulations of nonradial shock emergence in a frame comoving with the breakout pattern, focusing on the adiabatic, nonrelativistic limit in a plane stratified envelope. We set boundary conditions using a known selfsimilar 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 selfsimilar. Waves and vortices, which are not predominantly due to grid effects, emanate from this region. The postshock 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: general1. 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 whiplike motion described by the similarity solution of Gandel’Man & FrankKamenetskii (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 nonradial 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 nonradial. 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 nonradial 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 selfsimilar in the comoving frame: it could be oscillatory, or the outflow could interact with the star to spoil the apparent selfsimilarity. 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 comoving with the breakout pattern, and we can focus our attention on the twodimensional, stationary flow of an adiabatic gas, separated by a shock front from planestratified upstream fluid. Furthermore, we specialize to the case where the postshock flow is radiation dominated (adiabatic index ) and consider a cold polytropic atmosphere (index or ) for the preshock matter. Finally, we assume that the motions are purely nonrelativistic, 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, selfsilmilar 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 boxsize 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 timeaveraged sense) by the postshock 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 magnetohydrodynamics 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 onedimensional problem in which a shock front moves through cold gas away from a highpressure 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 selfsimilar theory of Sakurai (1960).
For our twodimensional simulations Athena was compiled with Message Passing Interface (MPI) and Static Mesh Refinement (SMR) enabled. We used a HartenLaxvan LeerContact (HLLC) Riemann solver with secondorder 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, selfsimilar 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 selfsimilarly, 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 timedependent, onedimensional flow into two dimensions. This involves several steps. First, we match the upstream density in the onedimensional solution to the depthdependent 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 seventhorder polynomial fit to the selfsimilar 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 preshock 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 ‘selfsimilar’ 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 postshock 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 selfconsistent 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 selfconsistent.
In the project’s initial stages we used a small wedge of very highdensity 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 selfsimilar 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 selfsimilar 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. Postshock 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 shortperiod 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 highestresolution 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 postshock 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.
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 nonradial, 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 selfsimilarly 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
(1) 
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 selfsimilar offset to the point of shock emergence in a definite way.) This deep selfsimilar 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 selfsimilar 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 selfsimilar (or both), it nevertheless suggests that, if the selfsimilar solution guides the flow in any way, its shock normal will be as close as possible to the selfsimilar 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 largeamplitude 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.
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 highestresolution subgrid, shown in Figure 8, reveals that KelvinHelmholz rolls develop along this vortex sheet.
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 entropyrelated 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 selfsimilar boundary conditions and the dynamics of smooth twodimensional 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.
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 timeaveraging 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, selfsimilar 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%,
(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 stellarframe 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 selfsimilar 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  

(A1)  
(A2)  
(A3)  
(A4)  
(A5) 
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
(3) 
The outflow speed in the star’s rest frame is
(4) 
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 boxsize 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 highpitched acoustic oscillations, those of periods , which would be expected from interactions between the flow and the discretized grid.
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 timeaveraged 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 breakup 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, selfsimilar 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. Selfsimilarity 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 surfaceskimming 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 finiteamplitude 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
(5) 
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 sixsevenths 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 finitevolume 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 point^{a}^{a}Location, 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 outflowing ejecta.  Kink altitude^{b}^{b}Altitude, in the final steady state, of a kink in the shock angle.  Match altitude^{c}^{c}Altitude, in the final steady state, of the bestfit 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 
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.
The structure of the shock, the shapes of the postshock streamlines, and the final entropy distribution are all very insensitive to . Oscillation frequencies are not strongly dependent, except that the lowerresolution runs show a suppression of some of the shortwavelength motions. It appears that there is very little difference between the highestresolution 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 nonradial 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 powerlaw 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 higherresolution 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 ejectaenvelope 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
Model  

RSG^{b}^{b}Model s15s7b2 of Woosley & Weaver (1995), provided by Stan Woosley.  14  490  1.20  0.566  0.520  
BSG^{c}^{c}Shigeyama & Nomoto (1990) model for the progenitor of SN 1987A, provided by Ken’ichi Nomoto.  15  50  3.88  0.0254  0.409  
Ic^{d}^{d} 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 
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 pressuredominated 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
(6) 
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 twodimensional constantvelocity flow, and must be truncated at because threedimensional 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
(7) 
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 nonradial 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 ,
and
therefore, with equation (6),
(8) 
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 deepflow limit (), we find in the Appendix that
(9) 
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 corecollapse 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 (), nonrelativistic 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 postshock 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 ejectastar shear flow, for later investigation.
6. Energy in nonradial motions
Some of the energy diverted by the oblique shock into nonradial 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 nonradial component of the energy flux is therefore . Computing the nonradial 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
(10) 
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, nonrelativistic 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), accretioninduced 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
(11) 
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 outflowing and returning matter is guaranteed to be complicated. Spherical and nonsteady 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 nonradial 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, twodimensional simulations of oblique supernova shock breakout, focusing on the limit of adiabatic, nonrelativistic 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 postshock flow around the breakout point is neither steady nor selfsimilar: 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 stronglydeflected) 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 nonradial 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 twodimensional simulations are undoubtedly more pronounced than they would be three dimensions, thanks to the wellknown inverse nature of 2D turbulence cascades (Kraichnan, 1967). As these vortices are superimposed on a rapidly expanding flow, we see little to suggest that threedimensional results would be qualitatively different; however this must be checked with future simulations.
While the simulations presented here exploit the local, essentially twodimensional 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 nonradial 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 nonradial 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, nonradial flows eliminate the shock breakout emission, which was expected in spherical theory, over much of the stellar surface. By deflecting ejecta, nonradial 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.
Appendix A Deep flow, smalldeflection 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, selfsimilar solution to shock acceleration and postshock 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 shallowangle flow ().
In the deep flow where corrections due to nonradial 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 smallangle limit.
In this limit each flow quantity will tend far downstream () toward a powerlaw 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
(A1) 
Mass conservation requires , so
(A2) 
The postshock entropy is , where and ; therefore
(A3) 
The pressure is so
(A4) 
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 ,
(A5) 
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
(A6) 
where and ; and
(A7) 
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.
References
 Fryer et al. (1999) Fryer, C., Benz, W., Herant, M., & Colgate, S. A. 1999, ApJ, 516, 892
 Gandel’Man & FrankKamenetskii (1956) Gandel’Man, G. M. & FrankKamenetskii, 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., RamirezRuiz, 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