THE SPIRAL MODES OF THE STANDING ACCRETION SHOCK INSTABILITY

The physical system consists of a steady-state, standing spherical accretion shock at a radial distance r_s0 from the origin. Below the shock, the material settles subsonically onto a star of radius r_* and mass M centered at the origin of the coordinate system. This settling is mediated by a cooling rate per unit volume that mimics neutrino emission due to electron capture, $\mathscr{L} \propto p^{3/2}\,\rho$, with p being the pressure and ρ being the density (e.g., Blondin & Mezzacappa 2006; Fernández & Thompson 2009b). In all the cases studied here, cooling is relevant only in a narrow region—of the order of a pressure scale height—outside the accretor. Upstream of the shock, the fluid is supersonic, with incident Mach number $\mathcal {M}_1 = 5$ at r = r_s0, and has a vanishing energy flux. Upstream and downstream solutions are connected by the Rankine–Hugoniot jump conditions. The equation of state is that of an ideal gas of adiabatic index γ = 4/3, and the self gravity of the flow is neglected. The equations describing this steady flow are presented in Appendix A.1. The model was first developed by Houck & Chevalier (1992) to study accretion onto compact objects and subsequently used by Blondin et al. (2003), Blondin & Mezzacappa (2006), Blondin & Mezzacappa (2007), Foglizzo et al. (2007), Blondin & Shaw (2007), and Fernández & Thompson (2009b) as a minimal approximation to the postbounce stalled shock.

In contrast to Fernández & Thompson (2009b), we do not include here the effects of nuclear dissociation at the shock. This is a significant energy sink in the system, changing the Mach number and density below the shock for fixed upstream conditions, and hence affecting the linear modes of the SASI. Nevertheless, as our goal here is to characterize the basic behavior of the instability in three dimensions, we aim at keeping the calculations as simple as possible.

Neutrino heating is also neglected, in order to suppress convection in the flow, as originally done by Blondin et al. (2003). Since most of the flow below the shock is adiabatic, the entropy profile is flat and hence marginally unstable to overturns triggered by the SASI itself.

We adopt a system of units normalized to the shock radius r_s0, the free-fall velocity at this radius $v_{\rm ff0} = \sqrt{2GM/{r_{\rm s0}}}$, and the upstream density ρ₁ (Fernández & Thompson 2009b). Numerical values relevant for the stalled shock phase of core collapse supernovae are r_s0 ∼ 150 km, v_ff(r_s0) ∼ 4.8 × 10⁹M^1/2_1.3(r_s0/150 km)^−1/2 km s⁻¹, t_ff ≡ r_s0/v_ff ∼ 3.1 M^−1/2_1.3(r_s0/150 km)^3/2 ms, and $\rho _1 \sim 4.4\times 10^7 \dot{M}_\mathrm{0.3}\,M_\mathrm{1.3}^{-1/2}\,(r_{s0}/150 \,{\rm km})^{-3/2}$ g cm⁻³ (assuming a strong shock), where the mass accretion rate has been normalized to $\dot{M} = 0.3\,\dot{M}_\mathrm{0.3}\,M_\odot$ s⁻¹. We will express quantities involving angular momentum in terms of $\dot{M}{r_{\rm s0}} \simeq 4\pi \times 10^{46}(r_{\rm s0}/150 \,{\rm km})^2 \dot{M}_\mathrm{0.3}$ g cm² s⁻¹.

To obtain linear SASI eigenmodes and eigenfrequencies, we solve the differential system of Foglizzo et al. (2007), the results of which agree with time-dependent axisymmetric hydrodynamic simulations to high precision (Fernández & Thompson 2009b). The relevant equations are described in detail in Appendix A; here we summarize the main elements of the calculation. The eigenvalue problem is formulated in terms of perturbation variables that are conserved when advected. These are the perturbed mass flux h = δρ/ρ + δv_r/v_r, Bernoulli flux f = v_rδv_r + δc²/(γ − 1), entropy δS = (γ − 1)⁻¹[δp/p − γδρ/ρ], and an entropy–vortex combination δK = r²v · (∇ × δw) + ℓ(ℓ + 1)(p/ρ)δS (see also Equation (A33)), with c² = γp/ρ, v the velocity and w the vorticity. Perturbations are decomposed in Fourier modes in time and spherical harmonics in the angular direction having the general form

$$\begin{equation} \delta q = \delta \tilde{q}(r)\,e^{-i\omega t}\,Y_\ell ^m(\theta,\phi), \end{equation} \tag{ 1 }$$

where $\delta \tilde{q}(r)$ is the complex amplitude, Y^m_ℓ is a spherical harmonic, and the complex frequency satisfies ω = ω_osc + iω_grow. This choice results in all thermodynamic variables along with δv_r being proportional to Y^m_ℓ (corresponding to spheroidal modes, e.g., Andersson et al. 2001). The transverse components of the velocity involve angular derivatives of the Y^m_ℓ and have the same radial amplitude (Appendix A),

$$\begin{eqnarray} \delta v_\theta & = & \delta \tilde{v}_\Omega (r)\,e^{-i\omega t}\,\frac{\partial }{\partial \theta }Y_\ell ^m(\theta,\phi) \end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray} \hspace{17pt}\delta v_\phi & = & \delta \tilde{v}_\Omega (r)\,e^{-i\omega t}\, \frac{1}{\sin \theta }\frac{\partial }{\partial \phi }Y_\ell ^m(\theta,\phi). \end{eqnarray} \tag{ 3 }$$

The system of ordinary differential equations that determines the complex amplitudes is independent of the azimuthal number m of the mode. Boundary conditions at the shock are expressed in terms proportional to the shock displacement Δξ, or the shock velocity Δv = −iωΔξ. By imposing $\delta \tilde{v}_r(r_*) = 0$ at the inner boundary, a complex eigenvalue ω is obtained for a given set of flow parameters. For each ℓ, a discrete set of overtones is obtained, which are related to the number of nodes in the radial direction (Foglizzo et al. 2007; Fernández & Thompson 2009b). The resulting eigenmodes have both real and imaginary eigenfrequencies.

For fixed γ, $\mathcal {M}_1$, and functional form of the cooling function, the eigenfrequencies of the SASI depend only on the relative size of the shock and accreting star, quantified by the ratio r_*/r_s0. Time-dependent studies of individual SASI modes are possible by choosing model parameters such that only the fundamental mode is unstable, for a given ℓ. Based on the fact that the relevant modes in more realistic core-collapse simulations are ℓ = 1, and 2, we adopt two configurations for single-mode studies: r_*/r_s0 = 0.5 for ℓ = 1, and r_*/r_s0 = 0.6 for ℓ = 2 (see Foglizzo et al. 2007 and Fernández & Thompson 2009b for more extended parameter studies of the SASI eigenfrequencies). In addition, we explore a configuration with a larger postshock cavity, r_*/r_s0 = 0.2, to more closely resemble actual stalled supernova shocks. In this case unstable overtones of ℓ = 1 and ℓ = 2 are present.

To facilitate visualization, and unless otherwise noted, we employ the real representation of spherical harmonics. For ℓ = 1, one has

$$\begin{eqnarray} Y_1^0 & = & \sqrt{\frac{3}{4\pi }}\cos \theta \end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray} Y_1^1 & = & \sqrt{\frac{3}{4\pi }}\sin \theta \,\cos \phi \end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray} Y_1^{-1} & = & \sqrt{\frac{3}{4\pi }}\sin \theta \,\sin \phi. \end{eqnarray} \tag{ 6 }$$

This basis is entirely equivalent to the one involving complex exponentials in azimuth (e.g., Arfken & Weber 2005), but has a straightforward interpretation: Equations (4)–(6) are dipoles in the z-, x-, and y-directions, respectively. For ℓ = 2, we have

$$\begin{eqnarray} Y_2^0 & = & \sqrt{\frac{5}{4\pi }}\frac{1}{2}(3\cos ^2\theta -1)\\ \ms Y_2^1 & = & \sqrt{\frac{15}{4\pi }}\sin \theta \cos \theta \,\cos \phi\\ \ms Y_2^{-1} & = & \sqrt{\frac{15}{4\pi }}\sin \theta \cos \theta \,\sin \phi\\ \ms Y_2^2 & = & \sqrt{\frac{15}{16\pi }}\sin ^2\theta \, \cos 2\phi\\ \ms Y_2^{-2} & = & \sqrt{\frac{15}{16\pi }}\sin ^2\theta \, \sin 2\phi, \end{eqnarray} \tag{ 7 }$$

corresponding to the usual z-symmetric quadrupole (Y⁰₂), and a series of four-striped "beach balls" with alternating polarity and symmetry axis along $\hat{x}$ (Y¹₂), $\hat{y}$ (Y⁻¹₂), and $\hat{z}$ (Y^±2₂).

Whenever necessary, we will denote the traditional complex spherical harmonics as ϒ^m_ℓ.

We perform time-dependent hydrodynamic calculations to verify the evolution of linear SASI eigenmodes in three dimensions, covering the linear and nonlinear phases. To this end, we employ the publicly available code Zeus-MP (Hayes et al. 2006), which solves the Euler equations using a finite difference algorithm that includes artificial viscosity for the treatment of shocks. The default version of the code supports an ideal gas equation of state and point-mass gravity. We have extended it to account for optically thin cooling and a tensor artificial viscosity, for a better shock treatment in curvilinear coordinates (e.g., Stone & Norman 1992; Hayes et al. 2006; Iwakami et al. 2008a). Issues associated with implementing the latter are discussed in Appendix B.

The initial conditions are set by the steady-state accretion flow used in linear stability calculations. To prevent runaway cooling at the base of the flow, we impose a Gaussian cutoff in the cooling with entropy (Fernández & Thompson 2009b). The normalization of the cooling function is adjusted so that the Mach number at the surface of the accreting star is $\mathcal {M}\sim 10^{-2}$.

Spherical polar coordinates (r, θ, ϕ) are used, covering the whole sphere minus a cone of half-opening angle 5° around the polar axis. This prescription does not significantly alter the flow dynamics and ameliorates the severe Courant–Friedrichs–Lewy restriction around the polar axis (H-Th. Janka 2010, private communication). We present convergence studies in Appendix B showing that only ∼10% differences relative to using the full sphere are introduced with this prescription. Cells are uniformly spaced in the polar and azimuthal directions, while ratioed in radius. The choice of grid spacing is key to a achieve a stable background solution and to adequately capture linear growth rates, while minimizing the computational cost. The radial size Δr_min of the cells at the inner boundary determines how stable the unperturbed shock remains, because near hydrostatic equilibrium needs to be maintained. At the shock, the angular size determines how well growth rates are captured. Based on numerical experiments, we have found that Δr_min < 10⁻³r_s0 is required to obtain a stable shock within 500 dynamical times, and that 36 cells in the polar direction are the minimum needed to obtain a clean shock oscillation (ω_osc within ≲10% and ω_grow within ∼20% from the linear stability value; convergence tests are presented in Appendix B). We choose the ratio of radial spacing so that Δr = r_s0Δθ at the shock radius r_s0, and Δϕ = Δθ, resulting in a total of 56 radial (ℓ = 1), 48 polar, and 96 azimuthal cells. For ℓ = 2, we use 28 ratioed cells in radius from r_* to r_s0, and then use a uniform radial spacing until the outer boundary. This choice of grid spacing is made because the tensor artificial viscosity behaves best when applied over the longest cell dimension (e.g., Stone & Norman 1992), and hence choosing the three cell dimensions as equal as possible results in nearly isotropic dissipation except near the polar axis.² In addition, we perform one simulation at double the resolution in all dimensions (e.g., 128 × 96 × 192) to test the robustness of our results. This resolution approaches values comparable to existing two-dimensional simulations and can better capture the parasitic instabilities that mediate the saturation of the instability (Guilet et al. 2010). In Table 1, we summarize the modes evolved. With the exception of the high-resolution model (R5_L1P2_HR), all simulations were carried out until at least 200 t_ff0 well into the nonlinear phase.

The boundary conditions are reflecting at the inner radial boundary (r = r_*), periodic in the ϕ direction, reflecting at the polar boundaries (θ = {0, π} minus the cutout), and set equal to the upstream flow at the outer radial boundary (r = 4r_s0). The use of tensor artificial viscosity prevents the appearance of a numerical instability at the poles of the shock, as described by Iwakami et al. (2008b). This carbuncle instability (Quirk 1994) arises within a few dynamical times if the standard VonNeumann & Richtmyer (1950) prescription is used. However, even when including the tensor viscosity, we have found that a different numerical instability develops around the polar axis when reflecting boundary conditions and the full sphere are used. This is due to the fact that a wake is generated by the reflecting axis when the sloshing modes reach a large amplitude, resulting in a different radial velocity profile at opposite sides of the polar axis. The azimuthal velocity then develops a sawtooth instability around the axis. The cutout improves (but does not completely suppress) the development of this instability, to the point where it does not significantly alter the flow dynamics.

In order to trigger an individual spiral mode in our simulations, we drop overdense shells with angular dependence determined by real spherical harmonics (Equations (4)–(11)). We locate these shells in the upstream flow in such a way that when advected, they arrive at the shock with a time delay that corresponds to a relative phase Φ. In other words, the radial positions r₁ and r₂ of a two-shell combination satisfy

$$\begin{equation} \int _{r_1}^{r_2} \frac{{\mathrm{{\it d}}}r}{|v_r|} = \frac{\Phi }{\omega _{\rm osc}}. \end{equation} \tag{ 12 }$$

More than two modes are included straightforwardly, as only relative phases are relevant. This type of perturbation results in the excitation of an isolated spiral mode when the background accretion flow does not rotate, allowing comparison with eigenfunctions from linear stability analysis, and quantification of its contribution to the angular momentum redistribution as it grows exponentially in amplitude.

Figure 4 shows the azimuthal velocity field at the same altitudes above the equator as in Figure 3 for model R5_L11P2, for which only the fundamental mode is unstable. Despite the coarseness of the grid, the flow structure in the upper row of Figure 3 is clearly reproduced.³ The online version of the article contains an animated version of this figure.

Zoom In Zoom Out Reset image size

To test the effects of unstable overtones, we have also evolved a configuration with r_*/r_s0 = 0.2, which is closer in size to stalled supernova shocks (models R2_L11f and R2_L11h). In this case ℓ = 1 has several unstable overtones (see, e.g., Figure 13 of Fernández & Thompson 2009b), which are all excited whenever an ℓ = 1 perturbation is imposed. Figure 5 shows the resulting azimuthal velocity arising from perturbations of the type in Equation (12), employing phase differences corresponding to the fundamental and first overtones of the Y¹₁ + iY⁻¹₁ mode. In the first case, shown in the left panel, a spiral-like pattern is formed, although the presence of unstable overtones makes it irregular. In the second, on the right, an initial first overtone spiral becomes overtaken by a fundamental sloshing mode in an oblique direction, because the latter mode has a larger growth rate.

Zoom In Zoom Out Reset image size

As a mode grows in amplitude, the pressure perturbation shown in Figure 1(k) develops into a protuberance that rotates with the oscillation period of the instability. Eventually a shock forms at the interface where the two streams move toward each other, resulting in a triple point at the shock (Blondin & Shaw 2007). Figure 6 shows the morphology of the shock right before this triple point forms, for models R5_L11P2, R6_L21P2, and R6_L22P2.

Zoom In Zoom Out Reset image size

The initial phase of exponential growth of isolated modes is followed by nonlinear mode coupling and saturation. Figures 7(a) and (b) show the spherical harmonic coefficients

$$\begin{equation} a_{\ell,m}(t) = \int {\mathrm{{\it d}}}\Omega \, Y_\ell ^m\, R_s(\theta,\phi,t), \end{equation} \tag{ 13 }$$

where R_s(θ, ϕ, t) is an isobaric surface tracing the shock, as a function of time for models R5_L11P2c and R5_L11P4. In both cases, the primary sloshing mode combination remains isolated in the linear phase, along with a non-oscillatory increase in the ℓ = 2, m = 0 mode due to the increasing oblateness of the shock, likely caused by the redistribution of angular momentum (Section 4.2). Other modes are excited only after the primary mode reaches nonlinear amplitudes. Interestingly, the first modes to couple are Y^±2₂, which achieve amplitudes about 50% lower than the primary mode. The axisymmetric ℓ = 1 mode together with ℓ = 2, m = {0, ± 1} modes couple much later, and do not reach as large an amplitude as m = 2. We do not show ℓ = 3 results to keep Figure 7 legible, but a similar trend is observed in that Y^±3₃ couple first and reach the largest amplitudes of all ℓ = 3 modes. This sequence of mode coupling differs from that found by Iwakami et al. (2008a), who witnessed the ℓ = 2, m = 0, ± 1 modes to couple first and reach the largest amplitude behind ℓ = 1. However, their non-axisymmetric perturbation excites the real spherical harmonics Y⁰₁ and Y¹₁ in phase, hence it is not surprising that the mode coupling sequence differs.

Zoom In Zoom Out Reset image size

To track how the relative phase in the primary spiral mode evolves going into the nonlinear phase, we compute a normalized cross-correlation between the spherical harmonic coefficients of the corresponding sloshing modes:

$$\begin{eqnarray} && \mathrm{CC}[a_{1,1},a_{1,-1}](\Delta t) \nonumber \\ &&= \frac{\int dt\,[a_{1,1}(t)-\bar{a}_{1,1}] [a_{1,-1}(t+\Delta t)-\bar{a}_{1,-1}]}{\big(\int dt[a_{1,1}(t)-\bar{a}_{1,1}]^2\big)^{1/2}\big(\int \,dt[a_{1,-1}(t)-\bar{a}_{1,-1}]^2\big)^{1/2}},\qquad \end{eqnarray} \tag{ 14 }$$

where $\bar{a}_{\ell,m}$ is the time average of a_ℓ,m and the integrals are performed over the same time interval [t_min, t_max]. Figure 8 shows the results of applying Equation (14) to the linear and nonlinear phases of modes R5_L11P2c, R5_L11P4, and R5_L11P4, which differ only in the initial relative phase of the perturbation. In the linear regime, the cross-correlation has an initial peak at a time delay Δt₁ = Φ/ω_osc. As this amounts to shifting the second sloshing mode back into the first, Equation (14) yields values near unity.⁴ Subsequent peaks arise at Δt_n = (Φ + nπ)/ω_osc, (n = 1, 2, ...), with an amplitude (−1)ⁿexp(−nπω_grow/ω_osc), i.e., each half-period down in magnitude by the inverse exponential growth rate, with alternating sign. When Equation (14) is applied to the fully saturated nonlinear phase, the first peak shifts to a larger value of Δt relative to the linear phase, with no significant change in amplitude between subsequent peaks. Furthermore, the three nonlinear cross-correlations are nearly identical, with a period ≃12 t_ff0 and a time delay ≃2.5 t_ff0, corresponding to a phase shift ≃2π/5. This suggests that the development of spiral modes in the fully saturated stage is independent of the initial conditions, for any non-zero initial relative phase between sloshing modes. As mentioned in Section 3, a larger class of initial perturbations can therefore result in spiral modes, not just those tuned to a relative phase of π/2.

Zoom In Zoom Out Reset image size

We have also performed a simulation in which, instead of using overdense shells, we impose random cell-to-cell pressure perturbations below the shock (R5_RAN). For this, we employ the setup with r_*/r_s0 = 0.5, for which only the fundamental ℓ = 1 mode is unstable (and the most unstable for all ℓ). Figure 9(a) shows the corresponding spherical harmonic coefficients as a function of time. Initially, a sloshing mode along the z-axis is triggered due to the form of the grid. However, as this mode reaches the nonlinear phase, a sloshing mode along the y-axis emerges out-of-phase, and then another along the x-axis follows. Figure 9(b) shows the normalized cross-correlation in the nonlinear phase, for the three different ℓ = 1 mode combinations over the time interval t ∈ [200, 300] t_ff0, showing that all of them are out-of phase. Note that the phase delay between Y¹₁ and Y⁻¹₁ is nearly the same as that in Figure 8 (modulo half a period), even though Y¹₁ is still growing in amplitude. The resulting angular momentum redistribution is discussed in Section 4.2.

Zoom In Zoom Out Reset image size

To close this subsection, we address the effects of dimensionality and resolution on the saturation amplitude. Table 2 shows the rms fluctuation of the a_1,1 spherical harmonic coefficient around its mean for various models. There does not appear to be a systematic difference between the saturation amplitude in two dimension and three dimension at this resolution and with this numerical method. The changes due to the presence of the cutout around the polar axis or a resolution doubling are of the same magnitude, ∼10%. Similarly, whether the sloshing mode is isolated or part of a spiral mode seems to make a small difference, which is a weak function of the relative phase between modes. In the context of parasitic instabilities (Guilet et al. 2010), this indicates that either (1) the resolution is still too low to expose the different growth rates of parasites in two dimension and three dimension or (2) given the structure of these modes, the difference in the growth rates does not translate in sizable differences in the saturation amplitudes.

Given that the background accretion flow has not net angular momentum, any net spin-up of the protoneutron star via accreted matter involves the separation of the postshock flow into regions with angular momentum of opposite sign (e.g., Blondin & Shaw 2007). In what follows, we explore how the linear and nonlinear phases of spiral modes mediate this angular momentum redistribution and quantify the magnitude of the effect. To ease comparison with other studies, we express angular momenta in units of $\dot{M}{r_{\rm s0}}^2$, which amounts to dividing our results by 4π within the dimensionless unit system described in Section 2.1.

4.2.1. Linear Phase

Figure 10 shows the z-component of the angular momentum density integrated over a spherical surface at radius r

$$\begin{equation} l_z(r) = r^2\int \, {\mathrm{{\it d}}}\Omega \, \rho \, r\sin \theta \, v_\phi, \end{equation} \tag{ 15 }$$

for model R5_L11_HR at various times during the linear phase. The upper panel shows times covering a full oscillation cycle, and the bottom shows curves rescaled by a factor exp [ − 2ω_grow(t − t₀)] over a longer timescale, with ω_grow the growth rate measured from the spherical harmonic coefficients. We thus infer that, in the linear phase, (1) a spiral mode has a characteristic angular momentum density profile which separates the flow into two counterrotating regions, (2) this profile is non-oscillatory, and (3) it grows in amplitude at nearly twice the growth rate of the mode. The latter is consistent with the fact that l_z(r) is a second-order quantity, as the zeroth- and first-order components of Equation (15) vanish. To lowest order in Δξ/r_s0, it is made up of products of the form Δρ⁽¹⁾Δv⁽¹⁾_ϕ and ρ⁽⁰⁾Δv⁽²⁾_ϕ, where the superscripts in parentheses indicate explicitly the order of the perturbation.

Zoom In Zoom Out Reset image size

We can separate out Equation (15) into components

$$\begin{eqnarray} l_z(r) & = & l^{(0)(2)} + l^{(1)(1)} + {\rm O(3)} \end{eqnarray} \tag{ 16 }$$

$$\begin{eqnarray} l^{(0)(2)} & = & r^3\int \,{\mathrm{{\it d}}}\Omega \, \rho ^{(0)}\, \sin \theta \Delta v_\phi ^{(2)} \end{eqnarray} \tag{ 17 }$$

$$\begin{eqnarray} l^{(1)(1)} & = & r^3\int \,{\mathrm{{\it d}}}\Omega \, \Delta \rho ^{(1)}\,\sin \theta \Delta v_\phi ^{(1)}, \end{eqnarray} \tag{ 18 }$$

where O(3) denotes higher order terms. Due to the fact that ρ⁽⁰⁾ is spherically symmetric, only the ℓ = 0 component of sin θ Δv⁽²⁾_ϕ survives:

$$\begin{equation} \big[\sin \theta \Delta v_{\phi }^{(2)}\big]_{0,0} = \int \,{\mathrm{{\it d}}}\Omega \, Y_0^0\,\sin \theta \,v_\phi - {\rm O(3)}. \end{equation} \tag{ 19 }$$

The first-order components are the real parts of the complex perturbations

$$\begin{eqnarray} \Delta \rho ^{(1)} & = & \frac{1}{2}\sum _{\ell,m}\left[ \delta \tilde{\rho }_{\ell,m}\Upsilon _\ell ^m e^{-i\omega _{\ell,m} t} + \delta \tilde{\rho }^*_{\ell,m}\Upsilon _\ell ^{m*} e^{i\omega ^*_{\ell,m}t}\right] \end{eqnarray} \tag{ 20 }$$

$$\begin{eqnarray} \sin \theta \Delta v_\phi ^{(1)} & = & \frac{i}{2}\sum _{\ell,m}m\left[ \delta \tilde{v}_{\Omega,\ell,m}\Upsilon _\ell ^m e^{-i\omega _{\ell,m} t}\right.\nonumber \\ & & - \left.\delta \tilde{v}^*_{\Omega \ell,m}\Upsilon _\ell ^{m*} e^{i\omega ^*_{\ell,m}t}\right], \end{eqnarray} \tag{ 21 }$$

where ϒ^m_ℓ denotes the usual complex spherical harmonics, and the star stands for complex conjugation. Integrating over angles yields

$$\begin{eqnarray} l^{(1)(1)} &=& r^3\sum _{\ell,m} m\, e^{2{\rm Im}(\omega _{\ell,m})t}\, [ {\rm Im}(\delta \tilde{\rho }_{\ell,m}){\rm Re}(\delta \tilde{v}_{\Omega,\ell,m})\nonumber \\ & & - {\rm Re}(\delta \tilde{\rho }_{\ell,m}){\rm Im}(\delta \tilde{v}_{\Omega,\ell,m})] +{\rm O(3)}. \end{eqnarray} \tag{ 22 }$$

Figure 11 shows l_z(r) at time t = 23.5 t_ff0 for model R5_L11_HR, where the Y¹₁ + iY⁻¹₁ mode is at phase zero (cf. Figure 7). Also shown are l⁽⁰⁾⁽²⁾, with both ρ⁽⁰⁾ and sin θ Δv⁽²⁾_ϕ obtained by projection of the corresponding fields onto Y⁰₀, and the ℓ = 1, m = 1 component of l⁽¹⁾⁽¹⁾. For the latter, the real part of the complex amplitude at phase zero is the projection onto Y¹₁, while the imaginary part is minus the Y⁻¹₁ component. Most of the angular momentum density comes from l⁽⁰⁾⁽²⁾, with l⁽¹⁾⁽¹⁾ being a ∼10% correction. Together, these terms account for almost all of l_z(r), in agreement with the fact that modes with ℓ ⩾ 2, m ≠ 0 have a negligible amplitude in the linear phase.

Zoom In Zoom Out Reset image size

For comparison, Figures 12(a) and (b) show the Y^±1₁ components of the density and azimuthal velocity fields as a function of radius at t = 23.5 t_ff0, for model R5_L11_HR, along with the eigenmodes from linear stability analysis normalized by the value of their respective spherical harmonic coefficients. Reasonable agreement is found between both approaches. Thus, most of l⁽¹⁾⁽¹⁾ can be accounted for by linear theory. Figure 12(c) shows [sin θ v⁽²⁾_ϕ]_0,0 (Equation (19)). Note that the amplitude is ∼10 times smaller than the first order ℓ = 1 components in Figure 12(a).

Zoom In Zoom Out Reset image size

The angular momentum redistribution is a smooth process, mediated by exponential amplification. Thus, the formation of a triple point at the shock is not the main agent behind the spin-up of the inner region, as we have shown that angular momentum redistribution begins as soon as linear modes are excited. The triple point may, however, be related to the saturation of the primary mode, setting the magnitude of the angular momentum redistribution.

A key element for a predictive theory of the spin-up due to exponentially growing spiral modes involves calculation of the second-order velocity perturbation sin θ v⁽²⁾_ϕ, since l⁽⁰⁾⁽²⁾ dominates the spin-up. We surmise that this would require repeating the steps followed for linear perturbations (Appendix A.1), but now expanding to second order in Δξ/r_s0. Combining this result with a criterion for the saturation of the SASI (such as that from Guilet et al. 2010) would allow estimation of the maximum angular momentum redistribution achievable by this means.

4.2.2. Nonlinear Phase and Total Spin-up

We focus first on the high-resolution ℓ = 1 model, R5_L11_HR, and then present results for other parameter combinations (Table 1). Figure 13 shows the evolution of the three components of the enclosed angular momentum

$$\begin{equation} L_i(r) = \int _{r_*}^r\, l_i(s)\,{\mathrm{{\it d}}}s\, \quad (i=x,y,z) \end{equation} \tag{ 23 }$$

evaluated at r = 0.6r_s0, corresponding to the approximate radius at which l_z(r) changes sign (see Figure 10). Two clear phases can be distinguished. First, from t = 0 to 50 t_ff0, L_z grows exponentially due to the effects discussed in Section 4.2.1. Once other modes start to grow in amplitude, the torque becomes oscillatory and saturates together with the primary mode. Once full saturation has been reached, the transverse components L_x and L_y increase in magnitude and fluctuate stochastically around zero. The subsequent evolution in this fully saturated stage depends somewhat on the resolution employed. Model R5_L11_HR seems to remain in a quasi steady state until the end of the integration at t = 180 t_ff0, while model R5_L11P2c shows large amplitude and long-period oscillations in L_z, and a secular increase in the magnitude of L_x and L_y at late times.

Zoom In Zoom Out Reset image size

Figure 14 shows how the angular momentum redistribution depends on the polar and azimuthal indices of the primary spiral mode, as well as on the relative phase with which it is excited. In this case, the ratio r_*/r_s0 is chosen so that only the fundamental mode is unstable for either ℓ = 1 or ℓ = 2, hence the primary spiral mode grows isolated in the linear phase. The top panel focuses on the low-resolution runs that excite the ℓ = 1, m = ±1 modes (R5_L11P2c, R5_L11P4, and R5_L11P8). The maximum spin-up due to the primary spiral mode (along L_z) is weakly dependent on the initial relative phase and has a magnitude ${\sim }0.6\dot{M}{r_{\rm s0}}^2$. However, the time delay required to achieve this maximum can differ by a factor of two. The transverse components (L_x and L_y) do not show a significant difference, aside from the secular growth in the case Φ = π/2, which is likely due to resolution effects (Figure 13). Figure 14(b) focuses on the ℓ = 2, m = ±1 modes, which still impart a spin-up along the z-axis despite having vanishing azimuthal velocity at the equator. Note however that the magnitude of the spin-up is much smaller than ℓ = 1 spiral modes. The magnitude of the transverse components at late times for Φ = π/2 is to be taken with caution, again due to possible resolution effects. Finally, the ℓ = 2, m = ±2 modes show a behavior qualitatively similar to ℓ = 1, m = ±1, although the magnitude of the maximum angular momentum redistribution depends strongly on the relative phase. As with ℓ = 2, m = ±1, the maximum spin-up is much smaller in magnitude than ℓ = 1.

Zoom In Zoom Out Reset image size

The effects of changing the size of the postshock cavity (r_*/r_s0) on the angular momentum redistribution are shown in Figure 15 for model R2_L11fc. Here, as in Figure 5(left), the system has a size more similar to realistic core-collapse situations, but is such that many ℓ = 1 overtones are unstable. Still, a spin-up is observed, which is in fact larger than the case r_*/r_s0 = 0.5. The secular growth of the spin-up, however, is to be taken with caution due to the low resolution of the model. The angular momentum redistribution following the phase of linear growth is still $\simeq 0.6 \dot{M}{r_{\rm s0}}^2$. This initial peak, however, comes after a long time delay of 130 t_ff0 ≃ 400M^−1/2_1.3(r_s0/150 km)^3/2 ms.

Zoom In Zoom Out Reset image size

The effect of changing the initial perturbation is shown in Figure 16, which shows the angular momentum enclosed within r = 0.6r_s0 for model R5_RAN. The angular momentum amplification is directly related to the growth of spiral modes (compare with Figure 9). By the time integration stops, Y¹₁ is still growing, so L_y may increase even more in magnitude. The fact that the spin-up along the z-axis is positive (in contrast to models R5_L11_HR and R5_L11P2c) is a consequence of the random initial relative phase between the corresponding sloshing modes. The spherical symmetry of the problem implies that adding π to Φ in Equation (12) leads to angular momentum redistribution with the opposite sign.

Zoom In Zoom Out Reset image size

The order of magnitude of the maximum angular momentum redistribution, $\dot{M}{r_{\rm s0}}^2$, can be understood if at saturation the azimuthal velocity v_ϕ is close to the upstream velocity v₁. From linear theory (Appendix A.1), this follows for shock displacements of order unity, which is indeed the case for ℓ = 1 modes. Not so straightforward to explain is the size of the fluctuations in the fully saturated state, which involves interference between modes with different amplitudes and phases. The basic behavior of the spin-up during the exponential growth phase agrees with what was seen by Blondin & Shaw (2007, their Figure 7), although they did not quantify their results in terms of fundamental physical quantities. Our results agree to within a factor of a few with those of Blondin & Mezzacappa (2007). They report typical spin-ups of 2.5 × 10⁴⁷ at 250 ms (72 t_ff0 for M = 1.2 M_☉ and r_s0 = 230/1.5 = 153 km, where 1.5 is our approximate average shock radius in units of r_s0 and 230 km their quoted shock position). At similar times, we find $0.6 \dot{M}{r_{\rm s0}}^2 \simeq 10^{47}(\dot{M}/0.36\,M_\odot \,{\rm s}^{-1})({r_{\rm s0}}/153 \,{\rm km})^2$ g cm² s⁻¹. The differences could be caused by several factors: (1) their absorbing boundary condition with no cooling, which causes the shock to expand with time (e.g., Blondin et al. 2003) and (2) resolution, which we have already shown to affect the long-term behavior of the system.

Assuming that all of the angular momentum enclosed within the radius where l_z(r) changes sign is accreted onto the neutron star, and that the latter has a moment of inertia I_NS = I₄₅ × 10⁴⁵ g cm², one can write the minimum period due to spiral modes as

$$\begin{equation} P \simeq 80\, I_{45} \dot{M}^{-1}_{0.3}\left(\frac{150 \,{\rm km}}{{r_{\rm s0}}} \right)^2 \left(\frac{0.6}{f_{\rm amp}} \right) {\rm ms}, \end{equation} \tag{ 24 }$$

where f_amp is the fraction of $\dot{M}{r_{\rm s0}}^2$ achieved during the phase of exponential growth. If no exponentially growing spiral mode takes place, one can still expect to achieve f_amp ≲ 0.1 from stochastic fluctuations, in which case the spin period P ≳ 500 ms. Figure 17 shows the evolution of the spin period for a few representative modes, taking into account the three components of the angular momentum, $P = 2\pi I_{\rm NS} / \sqrt{\vphantom{A^A}\smash{\hbox{$L_x^2+L_y^2+L_z^2$}}}$, with the L_i measured at the radius at which l_z(r) changes sign.

Zoom In Zoom Out Reset image size

This simplistic picture is bound to change when more realistic physics is included. First, because it is not clear yet that spiral modes can easily develop in a convective gain region (e.g., Iwakami et al. 2008a), but also because the neutrinospheric radius, the mass accretion rate, and the shock radius change as a function of time. The onset of explosion further complicates the picture, as the mass accretion rate decreases and the postshock flow expands. Still, anisotropic cold and fast downflows piercing the convective region are expected to impart stochastic torques to the forming neutron star (Thompson 2000), an effect that we cannot capture with our current setup.

The steady accretion flow below the shock is obtained by solving the time-independent Euler equations, using an ideal gas equation of state of adiabatic index γ and the gravity of a point mass M,

$$\begin{eqnarray} \frac{\partial }{\partial r} (r^2 \rho v_r) & = & 0 \end{eqnarray} \tag{ A1 }$$

$$\begin{eqnarray} v_r \frac{\partial v_r}{\partial r} + \frac{1}{\rho }\frac{\partial p}{\partial r} + \frac{\it GM}{r^2} & = & 0 \end{eqnarray} \tag{ A2 }$$

$$\begin{eqnarray} \frac{\partial p}{\partial r} - c_s^2 \frac{\partial \rho }{\partial r} & = & (\gamma -1)\frac{\mathscr{L}}{v_r}. \end{eqnarray} \tag{ A3 }$$

The boundary conditions at r = r_s0 are given by the Rankine–Hugoniot jump conditions (e.g., Landau & Lifshitz 1987). The upstream flow is adiabatic and has Mach number $\mathcal {M}_1$ at r = r_s0. For a given ratio r_*/r_s0, the normalization of the cooling function is obtained by imposing v_r(r_*) = 0.

Foglizzo et al. (2007) introduce Eulerian perturbations of the form shown in Equation (1). To obtain a compact formulation of the system, the following perturbation variables are employed

$$\begin{eqnarray} h & = & \frac{\delta v_r}{v_r} + \frac{\delta \rho }{\rho } \end{eqnarray} \tag{ A4 }$$

$$\begin{eqnarray} f & = & v_r \delta v_r + \frac{\delta c^2}{\gamma -1} \end{eqnarray} \tag{ A5 }$$

$$\begin{eqnarray} \delta S & = & \frac{1}{\gamma -1}\left(\frac{\delta p}{p} - \gamma \frac{\delta \rho }{\rho } \right) \end{eqnarray} \tag{ A6 }$$

$$\begin{eqnarray} \delta K & = & r^2 \mathbf {v}\cdot (\nabla \times \delta \mathbf {w}) + \ell (\ell +1)\frac{p}{\rho }\delta S, \end{eqnarray} \tag{ A7 }$$

corresponding to the perturbed mass flux, energy flux, entropy, and an entropy–vortex combination, respectively. The differential system for the radial profile of linear perturbations is obtained by perturbing the time-dependent fluid equations to first order and rewriting the resulting equations in terms of (A4)–(A7), obtaining (Foglizzo et al. 2007)

$$\begin{eqnarray} \frac{\partial h}{\partial r} & = & \frac{i\delta K}{\omega r^2 v_r} + \frac{i\omega }{v_r(1-\mathcal {M}^2)} \left[\frac{\mu ^2}{c^2}f -\mathcal {M}^2h -\delta S \right] \end{eqnarray} \tag{ A8 }$$

$$\begin{eqnarray} \frac{\partial f}{\partial r} & = & \delta \left(\frac{\mathscr{L}}{\rho v_r}\right) + \frac{i\omega v_r}{(1-\mathcal {M}^2)}\left[h - \frac{f}{c^2} +\left(\gamma -1 + \frac{1}{\mathcal {M}^2}\right)\frac{\delta S}{\gamma }\right]\nonumber\\ \end{eqnarray} \tag{ A9 }$$

$$\begin{eqnarray} \frac{\partial \delta S}{\partial r} & = & \frac{i\omega }{v_r} \delta S + \delta \left(\frac{\mathscr{L}}{p v_r}\right) \end{eqnarray} \tag{ A10 }$$

$$\begin{eqnarray} \frac{\partial \delta K}{\partial r} & = & \frac{i\omega }{v_r}\delta K + \ell (\ell +1)\delta \left(\frac{\mathscr{L}}{\rho v_r}\right), \end{eqnarray} \tag{ A11 }$$

where $\mu ^2 = 1 - (1-\mathcal {M}^2)[\ell (\ell +1)c^2/r]/\omega ^2$ and $\mathcal {M}$ is the Mach number. The shock boundary conditions for Equations (A8)–(A11) are obtained from the perturbed shock jump conditions evaluated at the unperturbed shock position r_s0 (Foglizzo et al. 2007),

$$\begin{eqnarray} h_2 & = & \left(\frac{1}{v_2} - \frac{1}{v_1} \right)\Delta v \end{eqnarray} \tag{ A12 }$$

$$\begin{eqnarray} f_2 & = & \left(\frac{\mathscr{L}_1}{\rho _1 v_1} - \frac{\mathscr{L}_2}{\rho _2 v_2}\right)\Delta \xi +(v_2 - v_1)\Delta v \end{eqnarray} \tag{ A13 }$$

$$\begin{eqnarray} \delta S_2 & = & \frac{\rho }{p}\left[\frac{\mathscr{L}_1}{\rho _1 v_1} - \frac{\mathscr{L}_2}{\rho _2 v_2} - \left(\frac{{\it GM}}{{r_{\rm s0}}^2} - 2\frac{v_1 v_2}{{r_{\rm s0}}} \right)\left(1 - \frac{v_2}{v_1}\right) \right]\Delta \xi\nonumber \\ & & - \frac{v_1}{p/\rho }\left(1 - \frac{v_2}{v_1} \right)^2 \Delta v \end{eqnarray} \tag{ A14 }$$

$$\begin{eqnarray} \delta K_2 & = & \ell (\ell +1)\left[ (v_1 - v_2)\Delta v + f_2\right], \end{eqnarray} \tag{ A15 }$$

where Δξ is the shock displacement, Δv = −iωΔξ is the shock velocity, and the subscripts 1 and 2 refer to quantities above and below the shock, respectively. The reflecting boundary condition at r = r_* is

$$\begin{equation} \delta v_r = \frac{v_r}{(1 - \mathcal {M}^2)}\left[h + \delta S - \frac{f}{c^2} \right] = 0, \end{equation} \tag{ A16 }$$

resulting in a complex eigenvalue ω.

The transverse components of the perturbed Euler equation are

$$\begin{eqnarray} {-}i\omega \,\delta v_\theta + \frac{v_r}{r}\frac{\partial \,\delta v_r}{\partial \theta } + v_r\, \delta w_\phi + \frac{1}{\rho \, r}\frac{\partial \delta p}{\partial \theta } & = & 0 \end{eqnarray} \tag{ A17 }$$

$$\begin{eqnarray} &&-i\omega \,\delta v_\phi + \frac{v_r}{r\,\sin \theta }\frac{\partial \,\delta v_r}{\partial \phi } - v_r\, \delta w_\theta + \frac{1}{\rho \, r\, \sin \theta }\frac{\partial \delta p}{\partial \phi } = 0,\nonumber \\ \end{eqnarray} \tag{ A18 }$$

from which one can readily solve for δv_θ and δv_ϕ if the transverse components of the vorticity w are known. Note that, since δv_r and δp are proportional to Y^m_ℓ (Section 2.2), one has {δv_θ,  δw_ϕ} ∝ ∂Y^m_ℓ/∂θ, and {δv_ϕ,  δw_θ} ∝ ∂Y^m_ℓ/(sin θ ∂ϕ). Hence, the only additional quantity required is the radial profile of the vorticity.

Taking the curl of the Euler equation yields

$$\begin{equation} \frac{\partial \delta \mathbf {w}}{\partial t} - \nabla \times \left(\mathbf {v}\times \delta \mathbf {w}\right) - \frac{1}{\rho ^2}\left(\nabla \rho \times \nabla p\right) = 0. \end{equation} \tag{ A19 }$$

The inertial term can be expanded as

$$\begin{equation} \nabla \times (\mathbf {v}\times \delta \mathbf {w}) = v_r\left(\nabla _\Omega \cdot \delta \mathbf {w}_\Omega \right)\hat{r} -\frac{1}{r}\frac{\partial }{\partial r} \left(r\,v_r\,\delta \mathbf {w}_\Omega \right), \end{equation} \tag{ A20 }$$

where ∇_Ω · δw_Ω is the angular part of the divergence of the transverse vorticity, $\delta \mathbf {w}_\Omega = \delta w_\theta \, \hat{\theta }+ \delta w_\phi \,\hat{\phi }$. Denoting by B the baroclinic term, one obtains

$$\begin{eqnarray} B_r & = & 0 \end{eqnarray} \tag{ A21 }$$

$$\begin{eqnarray} B_\theta & = & {-}\frac{1}{\rho _0^2}\left[ \frac{1}{r\sin \theta } \frac{{\mathrm{{\it d}}}p_0}{{\mathrm{{\it d}}}r} \frac{\partial \,\delta \rho }{\partial \phi } -\frac{1}{r\sin \theta } \frac{{\mathrm{{\it d}}}\rho _0}{{\mathrm{{\it d}}}r} \frac{\partial \,\delta p}{\partial \phi }\right]\qquad \end{eqnarray} \tag{ A22 }$$

$$\begin{eqnarray} B_\phi & = & {-}\frac{1}{\rho _0^2}\left[ \frac{1}{r}\frac{{\mathrm{{\it d}}}\rho _0}{{\mathrm{{\it d}}}r} \frac{\partial \, \delta p}{\partial \theta } -\frac{1}{r} \frac{{\mathrm{{\it d}}}p_0}{{\mathrm{{\it d}}}r} \frac{\partial \,\delta \rho }{\partial \theta }\right].\qquad \end{eqnarray} \tag{ A23 }$$

The equation governing the radial profile of the transverse vorticity is then

$$\begin{equation} \frac{\partial \delta \mathbf {w}_\Omega }{\partial r} = \left[\frac{i\omega }{v_r} + \frac{1}{r} + \frac{1}{\rho _0}\frac{{\mathrm{{\it d}}}\rho _0}{{\mathrm{{\it d}}}r}\right]\delta \mathbf {w}_\Omega -\frac{1}{v_r}\mathbf {B}. \end{equation} \tag{ A24 }$$

Once the radial profile of δw_Ω is known, solving for δw_r is straightforward, as no radial derivatives are involved.

The boundary condition at the shock for Equation (A24) is obtained from Equations (A17) and (A18), once the transverse velocity below the shock is known. Imposing $\hat{t}_i \cdot (\mathbf {v}_2\, +\, \delta \mathbf {v}_2\, - \, \mathbf {v}_1) = 0$ at the shock, where $\hat{t}_i$ are the tangent vectors

$$\begin{eqnarray} \hat{t}_\theta & = & \frac{1}{r}\frac{\partial \Delta \xi }{\partial \theta }\hat{r} + \hat{\theta } \end{eqnarray} \tag{ A25 }$$

$$\begin{eqnarray} \hat{t}_\phi & = & \frac{1}{r\sin \theta }\frac{\partial \Delta \xi }{\partial \phi }\hat{r} + \hat{\phi }, \end{eqnarray} \tag{ A26 }$$

one obtains (e.g., Foglizzo et al. 2007)

$$\begin{eqnarray} \delta v_{\theta,2} & = & \frac{v_{r,1} - v_{r,2}}{{r_{\rm s0}}}\frac{\partial \Delta \xi }{\partial \theta } \end{eqnarray} \tag{ A27 }$$

$$\begin{eqnarray} \delta v_{\phi,2} & = & \frac{v_{r,1} - v_{r,2}}{{r_{\rm s0}}\sin \theta }\frac{\partial \Delta \xi }{\partial \phi }, \end{eqnarray} \tag{ A28 }$$

where the subscripts 1 and 2 denote values upstream and downstream of the shock, respectively. We then obtain

$$\begin{eqnarray} \delta \tilde{w}_{r,2} & = & 0 \end{eqnarray} \tag{ A29 }$$

$$\begin{eqnarray} \delta \tilde{w}_{\theta,2} & = & -\frac{i\omega }{v_r}\delta \tilde{v}_{\phi,2}+\frac{1}{{r_{\rm s0}}}\delta \tilde{v}_{r,2} + \frac{1}{\rho {r_{\rm s0}}}\delta \tilde{p}_2 \end{eqnarray} \tag{ A30 }$$

$$\begin{eqnarray} \delta \tilde{w}_{\phi,2} & = & \frac{i\omega }{v_r}\delta \tilde{v}_{\theta,2}-\frac{1}{{r_{\rm s0}}}\delta \tilde{v}_{r,2} - \frac{1}{\rho {r_{\rm s0}}}\delta \tilde{p}_2, \end{eqnarray} \tag{ A31 }$$

where tildes denote the radial amplitude (Equation (1)). Given that $\delta \tilde{v}_{\theta,2} = \delta \tilde{v}_{\phi,2}$, it follows that $\delta \tilde{w}_{\phi,2} = -\delta \tilde{w}_{\theta,2}$. With this boundary condition, Equations (A22)–(A24) imply that $\delta \tilde{w}_\theta (r) = -\delta \tilde{w}_\phi (r)$ for all r, and hence $\delta \tilde{v}_\theta (r) = \delta \tilde{v}_\phi (r)$ as well. For convenience, we adopt the notation

$$\begin{equation} \delta \tilde{v}_\Omega (r) \equiv \delta \tilde{v}_\theta (r) = \delta \tilde{v}_\phi (r). \end{equation} \tag{ A32 }$$

Given the radial and angular form of the transverse velocity perturbations, one obtains δw_r(r) = (∇ × v)_r = 0 everywhere.

Knowing that the transverse components of the vorticity have equal and opposite radial amplitudes, one can show from the definition of δK (Section 2.2) that

$$\begin{equation} \delta \tilde{w}_\theta = \frac{1}{r v_r}\left[\frac{\delta \tilde{K}}{\ell (\ell +1)} - \frac{p}{\rho }\delta \tilde{S} \right], \end{equation} \tag{ A33 }$$

that is, the amplitude of the vorticity can be directly obtained from the basic perturbation variables, without the need for solving Equation (A24). We have nevertheless integrated this equation as a self-consistency check.

Astrophysical hydrodynamic codes that employ finite difference algorithms to solve the Euler equations rely on artificial viscosity to broaden shocks over a few cells. This allows differencing of the flow variables without discontinuities, while obtaining the correct solutions to the jump conditions away from shocks. The standard practice is to use the prescription of VonNeumann & Richtmyer (1950), which only takes into account the spatial derivative of the velocity normal to the shock as an approximation to the divergence of the velocity field for determining compression. When curvilinear coordinates are used, however, this approach results in short wavelength oscillations behind the shock, and spurious heating of homologously contracting flows (e.g., Tscharnuter & Winkler 1979). To fix this, the artificial viscosity can be formulated as a coordinate-invariant tensor, which can correctly identify zones where the flow is compressed by becoming active only when ∇ · v < 0 (Tscharnuter & Winkler 1979). Additional constraints on the functional form come from requiring that the artificial viscosity does not act on homologously contracting or shear flows (Tscharnuter & Winkler 1979; Stone & Norman 1992; Hayes et al. 2006).

In the context of core-collapse calculations, a tensor form of the artificial viscosity in Zeus-MP has previously been employed by Iwakami et al. (2008a), Iwakami et al. (2009a), and Iwakami et al. (2009b). They found that unless the tensor formulation is used, the carbuncle instability (Quirk 1994) appears at the poles of the shock (Iwakami et al. 2008b).

We have independently implemented a tensor artificial viscosity in Zeus-MP along the lines of Stone & Norman (1992) and Iwakami et al. (2008a), with some slight modifications. First, we use the volume difference instead of metric coefficient times coordinate difference (Stone & Norman 1992). This removes singularities at the polar axis whenever the factor sin θdθ appears in the denominator. Second, we have made use of the traceless nature of the artificial viscosity tensor to rewrite the θ-component of its divergences as

$$\begin{eqnarray} \left(\nabla \cdot \mathbf {Q}\right)_{\,(2)} & = & \frac{1}{r\,\sin \theta }\frac{\partial }{\partial \theta }(\sin \theta \, Q_{22}) - \frac{Q_{33}}{r\sin \theta } \frac{\partial \sin \theta }{\partial \theta }\qquad \end{eqnarray} \tag{ B1 }$$

$$\begin{eqnarray} & = & \frac{1}{r\,\sin ^2\theta }\frac{\partial }{\partial \theta } (\sin ^2{\theta }\,Q_{22}) + \frac{Q_{11}}{r\sin \theta } \frac{\partial \sin \theta }{\partial \theta },\qquad \end{eqnarray} \tag{ B2 }$$

where the second expression is that used by Stone & Norman (1992) and Iwakami et al. (2008a). This allows to completely eliminate the metric coefficients of the denominator by implementing the volume difference described above. The net effect is to reduce the noise in the θ velocity in the supersonic regions close to the z-axis, without making any significant difference at or below the shock.

When evolving our setup with this artificial viscosity prescription, the flow remains spherically symmetric and smooth for an indefinite time as long as there are no perturbations. If, on the other hand, the VonNeumann & Richtmyer (1950) prescription is used, the carbuncle instability indeed appears out of numerical noise within ∼10 dynamical times, in agreement with Iwakami et al. (2008b).

However, when evolving sloshing SASI modes transverse to the z-axis, a runaway sawtooth instability of the ϕ component of the velocity develops around the axis whenever the amplitude (and thus the flow velocity) becomes large. This happens because the radial velocity at both sides of the axis is not the same, which in turn is caused by a wake generated by the reflecting axis. The reflection of the flow at the axis also compresses it, resulting in enhanced artificial viscosity dissipation at the axis, with irregular velocities. A similar phenomenon is observed in modes parallel to the z-axis, although the cause of this is not clear yet.

The cutout at the polar axis, described in the next subsection, suppresses the growth of this numerical instability to the point where meaningful results can be obtained. It does not, however, completely eliminate it. We proceed to the nonlinear phase based on the fact that the angular momentum contributed by regions around the shock and close to the polar axis is minor.

To investigate the reliability of results obtained with the axis cutout prescription described in Section 2.3, we performed two runs with half-opening angles 1° and 25, which fill the gap between models R5_L11P2 and R5_L11P2c. Figure 18(a) shows the spherical harmonic coefficients for modes Y¹₁ and Y⁻¹₁ as a function of time, for different values of the half-opening angle θ_min of the cone that is removed from the grid around the polar axis. Aside from a ∼10% reduction in the maximum amplitude, all curves follow nearly identical trajectories until t ≃ 80 t_ff0, where they start to diverge more noticeably. Panel (b) shows the angular momentum density integrated over a spherical surface at time t = 24 t_ff0 (Equation (15), a low-resolution counterpart to Figure 11). Again, as the opening angle increases, results decrease by about 10%, consistent with the decrease in amplitude of shock oscillations. Other than that, the qualitative behavior is identical. We conclude that our resolution-independent results are reliable, with a quantitative uncertainty of at least 10%.

Zoom In Zoom Out Reset image size

The eigenfrequencies obtained with the method described in Appendix A have been verified to high precision with time-dependent axisymmetric hydrodynamic simulations using the code FLASH2.5 (Fernández & Thompson 2009b). A basic test for the newly implemented setup is thus the reproduction of the linear SASI growth rates from linear stability analysis. This is a resolution-dependent requirement, hence convergence tests are needed. Figure 19 shows growth rates and oscillation frequencies as a function of angular and radial resolution, for axisymmetric simulations of ℓ = 1 and ℓ = 2 modes. The measurement method is the same as that described in Fernández & Thompson (2009b). With the baseline resolution employed in this study (N_θ = 48), growth rates are captured within ∼20% and oscillation frequencies within 10%.

Zoom In Zoom Out Reset image size