# On the linear growth mechanism driving the stationary accretion shock instability

###### Abstract

During stellar core collapse, which eventually leads to a supernovae explosion, the stalled shock is unstable due to the standing accretion shock instability (SASI). This instability induces large-scale non spherical oscillations of the shock, which have crucial consequences on the dynamics and the geometry of the explosion. While the existence of this instability has been firmly established, its physical origin remains somewhat uncertain. Two mechanisms have indeed been proposed to explain its linear growth. The first is an advective-acoustic cycle, where the instability results from the interplay between advected perturbations (entropy and vorticity) and an acoustic wave. The second mechanism is purely acoustic and assumes that the shock is able to amplify trapped acoustic waves. Several arguments favouring the advective-acoustic cycle have already been proposed, however none was entirely conclusive for realistic flow parameters. In this article we give two new arguments which unambiguously show that the instability is not purely acoustic, and should be attributed to the advective-acoustic cycle. First, we extract a radial propagation timescale by comparing the frequencies of several unstable harmonics that differ only by their radial structure. The extracted time matches the advective-acoustic time but strongly disagrees with a purely acoustic interpretation. Second, we present a method to compute purely acoustic modes, by artificially removing advected perturbations below the shock. All these purely acoustic modes are found to be stable, showing that the advected wave is essential to the instability mechanism.

###### keywords:

instabilities — waves — shock waves — supernovae: general — methods: analytical## 1 Introduction

Several observational as well as theoretical arguments show that most core collapse supernovae are asymmetric explosions. The observational evidence include polarisation measures (Leonard et al., 2006), double peaked oxygen lines (Maeda et al., 2008), and the large velocities of neutron stars (Lai, 2001). On the theoretical side, the failure of realistic spherically symmetric calculations to produce an explosion (Liebendörfer et al., 2001) suggests that the breaking of the initial spherical symmetry is essential to the explosion mechanism. Indeed most proposed mechanisms are fundamentally asymmetric: the modern version of the neutrino-driven mechanism (Marek & Janka, 2009), the acoustic mechanism (Burrows et al., 2006), and the magneto-centrifugal mechanism (Wheeler et al., 2002).

The spherical symmetry can be broken by several phenomena. In the magneto-centrifugal mechanism, the asymmetry is due to a fast rotation in conjunction with a strong magnetic field. In the case of moderate rotation, two non radial instabilities can develop during the stalled shock phase: the neutrino-driven convection, and the stationary accretion shock instability usually called SASI which induces non spherical oscillations of the shock (Blondin et al., 2003). Of these two instabilities, SASI is the most efficient in creating a global asymmetry (Foglizzo et al., 2006) and helping the explosion (Marek & Janka, 2009; Burrows et al., 2006). In addition to its crucial role in the explosion mechanism, SASI has a number of potentially important consequences: it affects the kick and the spin of the neutron star (Scheck et al., 2006; Blondin & Mezzacappa, 2007), causes an emission of gravitational waves (Kotake et al., 2007), and creates a time dependence of the neutrino signal (Marek et al., 2009).

But what is SASI ? Which physical phenomenon drives these shock oscillations ? Surprisingly, the answer to this basic question remains somewhat uncertain. Indeed, two mechanisms have been proposed – the advective-acoustic cycle (Foglizzo et al., 2007) and a purely acoustic mechanism (Blondin & Mezzacappa, 2006) – and distinguishing them remains a controversial issue. The aim of this article is to clarify this question.

In Section 2, we review previous works by explaining the two instability mechanisms and the arguments proposed to distinguish them. On the basis of these works, the advective-acoustic cycle is generally favoured but the purely acoustic mechanism cannot be ruled out for realistic flow parameters. In Section 3, we propose a new method to extract a radial propagation time from the frequency of the unstable modes. For this purpose we compare the frequencies of several harmonics that differ only by their radial structure. As the radial advective-acoustic and purely acoustic times are significantly different, this allows us to distinguish between the two cycles and discards the purely acoustic mechanism. In Section 4, we describe a method to compute purely acoustic modes, which are shown to be stable. In Section 5, we give a detailed analysis of the frequency spectra of advective-acoustic and purely acoustic cycles, which confirms the validity of the method used in Section 3. Finally in Section 6, we summarise our work and conclude that the instability mechanism behind SASI is the advective-acoustic cycle.

## 2 The two proposed mechanisms: summary of previous works

### 2.1 Different models of SASI

SASI has been observed to grow in many different models of stellar collapse, which have very diverse degrees of complexity. These include state of the art realistic numerical simulations (Marek & Janka, 2009; Burrows et al., 2006), as well as models where successive simplifications have been made: approximate neutrino transport (Scheck et al., 2008), neglect of the stellar structure (Ohnishi et al., 2006; Yamasaki & Yamada, 2007; Iwakami et al., 2008), simple equation of state with the effect of neutrinos parameterised by cooling/heating functions (Blondin et al., 2003; Blondin & Mezzacappa, 2006; Foglizzo et al., 2007; Yamasaki & Foglizzo, 2008; Fernández & Thompson, 2009b, a), adiabatic approximation (Blondin & Mezzacappa, 2007), and finally the most simplified of all these models is the planar toy model (Foglizzo, 2009; Sato et al., 2009). Very often, the simpler the model, the deeper the physical understanding. For example, models neglecting the heating have clearly shown that SASI is distinct from neutrino-driven convection (Blondin et al., 2003), a fact missed by most realistic models. The neglect of the stellar structure allows to set up an initial steady state, which can be studied with a perturbative analysis (Foglizzo et al., 2007; Yamasaki & Yamada, 2007). Finally, the simplest of these models, the planar toy model, allows for an analytical treatment where the instability mechanism can be clearly demonstrated.

It is widely accepted that the same instability is at work in all these models. Therefore, although we base our investigation on the model of Blondin & Mezzacappa (2006), we expect our conclusions to hold more generally. We choose this model because it is simple enough to perform a perturbative analysis with clean boundary conditions. This idealised flow was first described by Houck & Chevalier (1992) and later studied by Blondin & Mezzacappa (2006) to describe SASI. It has then become a classical model of SASI and has been studied also by Foglizzo et al. (2007); Fernández & Thompson (2009b); Fernández (2010); Guilet et al. (2010). In this flow, the fluid is modelled as a non self-gravitating perfect gas with an adiabatic index of 4/3. The effects of neutrinos are simply modelled by a cooling function, while neutrino heating is neglected. The proto-neutron star is treated as a hard surface on which the cooling fluid settles down. The shock is stationary at a radius which can be freely chosen by adjusting the normalisation of the cooling function. This model and its linear analysis will be extensively used in the remainder of the paper. We refer the reader to Foglizzo et al. (2007) for a description of the linear eigenmodes calculation, and to Appendix A for a slightly different (but equivalent) formulation of the equations governing it, which has been used in this paper.

### 2.2 The two mechanisms

The two mechanisms proposed to explain the growth of SASI are schematically illustrated by Figure 1. The advective-acoustic cycle (left panel) is a cycle between two waves : an advected wave composed of entropy and vorticity and an acoustic wave propagating outward (Foglizzo & Tagger, 2000; Foglizzo, 2002; Blondin et al., 2003; Foglizzo et al., 2007). Note that the acoustic wave also propagates in the transverse direction. These waves are coupled via two coupling processes taking place at the shock and at a smaller radius in the deceleration region. The acoustic wave reaching the shock makes it oscillate and creates the entropy-vorticity wave. This wave is advected towards the proto-neutron star (hereafter called PNS) and its deceleration creates an acoustic feedback thus closing the loop. The cycle is unstable if the amplitude of the wave at the end of a cycle is larger than it was at the beginning.

The purely acoustic mechanism has been described in Blondin & Mezzacappa (2006) as an acoustic wave trapped in the postshock cavity. This acoustic mode is assumed to be ampliÞed through its interaction with the shock. Even when the acoustic propagation is mostly azimuthal some radial propagation is also present between the shock wave and the inner turning point (Fig. 6 of Blondin & Mezzacappa (2006)). Such an acoustic path can be interpreted as a purely acoustic cycle between the shock and the turning point (Figure 1, right panel). We stress that this description is quite general and does not assume a purely radial propagation. The transverse propagation is not restricted or disregarded neither here nor anywhere in this paper.

Any of these two cycles could in principle be unstable depending on the efficiency of the different coupling processes. For example, the purely acoustic cycle might be unstable if the acoustic wave is amplified when it is reflected at the shock, as assumed by Blondin & Mezzacappa (2006). Alternatively, the advective-acoustic cycle can be unstable if gradients below the shock allow for an efficient acoustic feedback, and if the shock is strong enough to efficiently create an entropy-vorticity wave (Foglizzo & Tagger, 2000; Foglizzo, 2001, 2002, 2009). By studying a simplified model of the advective-acoustic cycle, Foglizzo (2009) could explain salient features of SASI, in particular that it is dominated by low frequency and large scale non-radial modes. On the other hand, the purely acoustic cycle lacks a reference model where the properties of such acoustic instability could be unambiguously determined. As a consequence, it has not so far explained the basic features of SASI.

### 2.3 Timescales

A large part of this article deals with timescales: in particular how they determine the oscillation frequency, and how this can be used to discriminate between the two mechanisms. It is therefore useful to recall the characteristic timescales of the two cycles. The duration of a radial advective acoustic cycle can be written as:

(1) |

where is the effective coupling radius where the acoustic feedback is produced^{1}^{1}1Note that this is a simplified description of the instability: in reality the feedback is distributed radially rather than produced at a single radius.. Foglizzo et al. (2007) and Scheck et al. (2008) have suggested that the coupling radius is close to the location where the deceleration is strongest, however there remains some uncertainty in the precise location of this coupling. Similarly, the duration of a radial purely acoustic cycle is:

(2) |

where is the radius where the acoustic wave propagating down is reflected (i.e. its turning point). Finally, another acoustic timescale which plays a role potentially for both of the cycles is the azimuthal acoustic time:

(3) |

The two radial timescales suffer some uncertainty due to the lack of knowledge about and , both of which can in principle depend non trivially on the mode considered. In the model studied in this paper and briefly described in Section 2.1, it is however possible to obtain a useful upper limit by supposing that these two radii coincide with the PNS surface^{2}^{2}2With the cooling function considered here (, ), the advection time until the PNS surface is finite. However with a less violent cooling law (, ) this time is infinite, so that no useful upper limit can be obtained. It is then more relevant to estimate the coupling radius to be where the velocity gradient is strongest (Foglizzo et al., 2007).. In the following, we choose this prescription, keeping in mind that a shorter time (hence a higher frequency) can be obtained by changing and . The azimuthal acoustic time also suffers an uncertainty due to the radius at which it should be estimated. Here again a useful upper limit can be obtained by choosing the shock radius .

### 2.4 WKB description of the two cycles (following Foglizzo et al. (2007))

In order to quantitatively describe the two cycles, one needs to separate perturbations below the shock into the three types of waves: two acoustic waves (propagating up and down) and the advected wave. Below the shock the flow is close to adiabatic (Foglizzo et al. (2007), and see Section 4), thus the entropy and vorticity are simply advected. The main difficulty comes from the acoustic perturbations, which can be separated between waves propagating up and down only in the WKB approximation. The WKB approximation is valid if the wavelengths of the acoustic and advected waves are small compared to the length scale of the gradients. It can therefore be applied only to the modes with a high enough frequency, i.e. satisfying:

(4) |

Foglizzo et al. (2007) performed such a WKB analysis of the flow considered by Blondin & Mezzacappa (2006). In the following, we summarise the principle of this analysis, which could in principle be applied more generally to other flows.

Constants describing the coupling efficiency can be defined in the following way. When an acoustic wave propagating up with an amplitude hits the shock, it creates an advected wave with an amplitude as well as an acoustic wave propagating down with an amplitude ^{3}^{3}3the definition of the variable is recalled in Appendix A. Choosing another variable to measure the amplitude would change the coupling constants at the shock and in the gradients, but not the global cycle constant determining the stability of the cycles.. The ratio of the amplitude of the waves created to that of the incoming wave defines the coupling constants: and . These can be computed analytically using the boundary conditions at the shock and by separating the perturbations into the three types of waves owing to the WKB method (Foglizzo et al., 2005; Foglizzo, 2009).

When an entropy-vorticity wave is advected towards the PNS, it creates an acoustic wave propagating up with an efficiency characterised by the ratio of the amplitudes (measured at the shock): . Similarly the reflection in the gradients of the acoustic wave propagating down is characterised by: . These two constants can be computed numerically by integrating the differential system describing the perturbations, imposing modified boundary conditions at the shock ( to obtain , to obtain ).

Finally, using these four coupling constants, one can define the global constants that describe the efficiency of the two cycles (the amplitude ratio between the beginning and the end of the cycles): and . The cycle constants are complex numbers that depend on the complex frequency . Since the two cycles take place simultaneously, an eigenmode must satisfy the following relation (Foglizzo & Tagger, 2000):

(5) |

If these constants are calculated with a real frequency then the instability criterion for a single cycle can be written: or (depending on the cycle considered), which is equivalent to the criterion stated in Section 2.2. The growth rate can then be approximated by:

(6) |

where is the duration of the cycle considered (here the advective-acoustic one). When a second cycle is present, it can interfere constructively or destructively with the first cycle depending on their relative phase. This causes oscillations of the growth rate, which are indeed observed in the eigenspectrum of SASI (Figure 7 of Foglizzo et al. (2007)). These oscillations contain information about the efficiency and the timescale of the two cycles at work. Foglizzo et al. (2007) could thus extract directly from the eigenspectrum the values of the constants and (note that this is independent from the WKB analysis, and only assumes the presence of two cycles). The value of the cycle constants obtained from the WKB method and from the analysis of the eigenspectrum are in very good agreement (their Figures 8 and 9). They show that the advective-acoustic can be unstable with a cycle constant reaching , while the purely acoustic cycle is always stable with typically .

### 2.5 Which cycle drives SASI?

The WKB analysis of Foglizzo et al. (2007) summarised in the last subsection has proven that the advective-acoustic cycle is unstable and can explain SASI, while the purely acoustic cycle is stable. However this result is valid only for high frequency modes, satisfying Equation 4. The first few harmonics do not satisfy this inequality as for the fundamental mode. One may consider that the WKB method can be safely applied for the tenth harmonics and higher (satisfying ). Unfortunately these high frequency modes are unstable only for large shock radii (depending on the cooling function chosen). Thus strictly speaking in the more realistic range of shock radius , the WKB analysis cannot conclude on the instability mechanism. Foglizzo et al. (2007) have however argued by a continuity argument that it is not necessary to invoke a different instability mechanism for the low frequency modes. Indeed, the eigenspectrum vary smoothly both with the shock radius and mode frequency, and the eigenfunction of the low frequency modes resembles that of the high frequency modes. It may thus seem natural to assume that all SASI modes are due to the advective-acoustic cycle. Another argument in this direction is that the WKB approximation usually gives quite good results even when the small parameter (here ) is of order unity. This is however not a formal proof, and it cannot be excluded that the low frequency modes at moderate shock radius originate from a different instability mechanism.

Laming (2007) has studied analytically the stability of a spherical accretion shock by establishing an approximate dispersion relation, using a method inspired from Vishniac & Ryu (1989). This equation contains terms due to advection, which Laming (2007) proposes to remove in order to assess the importance of the advective-acoustic cycle in the instability mechanism. The conclusion he initially reached by this mean depended on the radius of the shock. At large shock radii, the omission of the advective term stabilised the flow, thus favouring the advective-acoustic cycle. On the contrary, at small shock radii advection affected very little the growth rate, favouring a purely acoustic interpretation of the instability. However, this analysis has been criticised by Yamasaki & Foglizzo (2008) because the dispersion relation was established with the use of dubious approximations. Furthermore an erratum published later on (Laming, 2008) corrects a few mistakes and changes the conclusion: the advective-acoustic cycle is then favoured for all shock radii.

A number of authors have tried to distinguish the two mechanisms by using the oscillation frequency observed in the simulations. These studies have led to diverging conclusions: Ohnishi et al. (2006) concluded that the advective-acoustic cycle is responsible for the instability, while Blondin & Mezzacappa (2006) concluded in the favour of the purely acoustic mechanism. This confusion probably comes from the fact that both of the mechanisms can explain the oscillation frequency if one makes the right assumptions. Indeed, Ohnishi et al. (2006) excluded the purely acoustic mechanism because the oscillation period is much longer than the radial acoustic time, while Blondin & Mezzacappa (2006) have shown that the oscillation period could be explained if one considers a non radial acoustic path. On the other hand, Blondin & Mezzacappa (2006) excluded the advective-acoustic cycle without considering the possibility that the acoustic feedback be created at a larger radius than the PNS surface. The determination of the frequency of each cycle thus suffers some uncertainty due to our lack of knowledge of the coupling radius (for the advective-acoustic cycle), and of the non-radial path to be considered (for a purely acoustic cycle). Later, Scheck et al. (2008) have shown that the oscillation period observed in their simulations was compatible with both mechanisms if the advective-acoustic coupling radius lies where the deceleration is maximum, or if an ad hoc non radial acoustic path is chosen (their Figure 16).

Fernández & Thompson (2009b) have proposed an argument that can be applied to the regime in which the WKB method is inapplicable. They followed the variation of a SASI mode growth rate as the shock radius is varied, which changes the ratio of the azimuthal acoustic time to the advection time from the shock to the PNS. They observed that the mode has a maximum growth rate when the two times coincide, and interpreted this as an indication that the advective-acoustic cycle is at work (because both advection and acoustic times play a role in the instability). While this observation is indeed suggestive that the advection plays a role in determining the growth rate, it does not rule out the possibility of the acoustic cycle playing a dominant role with only a minor help from the advection.

Some more indirect information about the instability mechanism can be found by considering the effect of a magnetic field on SASI. Endeve et al. (2010) have performed numerical simulations of SASI in the presence of a radial magnetic field (a split monopole). They obtained the surprising result that a magnetic field, which magnetic pressure is negligible compared the thermal one (), is capable of entirely stabilising SASI (in their model 2DB14Am). This seems to contradict a purely acoustic mechanism because the acoustic waves should be almost unaffected by such a weak magnetic field. By contrast, this result can be naturally interpreted in the framework of the advective-acoustic cycle. Indeed Guilet & Foglizzo (2010) have shown that the advective-acoustic cycle is significantly affected by a magnetic field if the Alfvén speed is comparable to the advection velocity, which can happen even if the magnetic pressure is negligibly small. In their model 2DB14Am, the Alfvén speed is indeed comparable to the advection velocity in the vicinity of the PNS.

To summarise this subsection, the advective-acoustic mechanism is favoured by most arguments. However none of these arguments can rule out the purely acoustic mechanism for realistic flow parameters. The goal of this paper is to find a method to definitely distinguish the two mechanisms for realistic flow parameters.

### 2.6 This article

As discussed in the last subsection, previous studies were not successful in determining the instability mechanism from SASI oscillation frequency, because the fundamental mode frequency can be explained by both mechanisms. This is due to the fact, that the (radial) advective-acoustic time and the azimuthal acoustic time are both comparable to the oscillation period, as is illustrated by Figure 2. In Section 3 we present a method to avoid this difficulty by extracting a radial time from the frequencies of several unstable radial harmonics. The large difference between acoustic and advective-acoustic radial times allows us to distinguish clearly the two mechanisms (the radial acoustic time is shorter by a factor 3 to 10, as shown in Figure 2).

Another reason for the remaining ambiguity on the instability mechanism is that, for realistic flow parameters, the unstable modes cannot be described by the WKB method. In Section 4 we present a new method to compute purely acoustic modes, which partially avoids this difficulty. This method makes use of the fact that the wavelength of advected perturbations is shorter than the acoustic wavelength. As a consequence, the identification of advected perturbations through a WKB approximation is valid at lower frequency.

Finally, in Section 5 we develop an analytical description of the frequency spectrum of the two cycles, in order to confirm the validity of the argument presented in Section 3. This analytical treatment is checked against the acoustic modes computed in Section 4.

In the following we apply our two arguments to the classical model of SASI briefly described in Section 2.1. The reader is referred to Blondin & Mezzacappa (2006) and Foglizzo et al. (2007) for a more detailed description. We focus on a shock radius , which is typical of supernovae and for which the unstable modes cannot be described by the WKB method. We also choose the following parameters: , , .

## 3 Extraction of a radial propagation timescale from the eigenspectrum

Previous studies that tried to determine the instability mechanism using the oscillation frequency (Blondin & Mezzacappa, 2006; Ohnishi et al., 2006; Scheck et al., 2008) have based their discussion entirely on the mode that grew fastest in their simulation. This most unstable mode is usually the lowest frequency mode, and its frequency is consistent with a radial advective-acoustic time as well as with a non-radial acoustic time. The specificity of our study is that we consider other (less) unstable modes, which can provide us with additional constraints on the instability mechanism. These modes are accurately computed by a linear analysis (Foglizzo et al., 2007), but may be hard to distinguish clearly in numerical simulations (although they should be present) because they grow more slowly and thus have a smaller amplitude than the most unstable mode. We focus in particular on higher frequency harmonics having the same transverse structure as the most unstable mode. For most astrophysically relevant values of the shock radius, several of these modes are unstable. For example, for our fiducial value , three modes are unstable (Figure 10). We emphasise that these modes have a higher frequency than the fundamental mode, but not high enough to allow for a WKB analysis.

What timescale determines the frequency change of the higher harmonics as compared to the fundamental mode ? We argue that this timescale cannot correspond to an azimuthal propagation because all these modes share the same transverse structure (determined by the spherical harmonics index ). As illustrated by Figure 3, these modes differ only by their radial structure suggesting that the frequency change should instead be controlled by a radial time, either advective-acoustic or purely acoustic. As a consequence, we suggest that a radial propagation timescale can be extracted from the frequency difference between two successive modes that share the same transverse structure (i.e. same spherical harmonics). We define such a time in the following way:

(7) |

where and are the frequencies of the n-th and n+1-th mode of SASI for a given spherical harmonics index (ranking the modes from the lowest to the highest frequency). From the physical argument stated above, we expect this time to be of the order of the radial acoustic time if the instability mechanism were purely acoustic, or of the order of the radial advective-acoustic time if SASI were due to the advective-acoustic cycle. Given the qualitative nature of the argument, we may not however expect a precise agreement. Indeed the analysis of the frequency spectrum detailed in Section 5 confirms that this is only approximately true. However thanks to the large difference between radial times (Figure 2), this approximate estimate is enough to be conclusive. We note that this method bears some similarity with the extraction of the radial structure of the Sun using the ”large frequency spacing” in helioseismology (Deubner & Gough, 1984).

Some light can be shed on the above argument by drawing an analogy with the very classical case of acoustic modes in a rectangular box (illustrated by Figure 4). Let us assume that the box is filled with a uniform gas at rest and delimited by reflecting rigid walls. The frequency eigenspectrum of acoustic modes is then:

(8) | |||||

(9) |

where and are the back and forth acoustic propagation times in the and directions, while and are the number of half wavelengths in the width of the box (in and direction). To draw the analogy, we consider a box elongated in the -direction, with representing the radial direction and representing the transverse direction ( or ). The correspondences are thus , and . distinguishes the different modes that share the same spherical harmonics index. corresponds to the lowest frequency fundamental mode, with a frequency that is dictated by the transverse acoustic time. Higher integer values of correspond to the higher frequency modes, which are the focus of this section. As is increased the frequency increases due to the first term in Equation (8), which contains the acoustic time in the -direction (analogous to the radial acoustic time). This confirms the qualitative argument given in the last paragraph. As the frequency is not simply proportional to , the frequency difference between two successive modes is not exactly the ”radial” frequency , suggesting that the Equation (7) is only an approximate estimate even in this simplest case. However the extracted timescale does tend toward the radial time in the limit where . As the azimuthal acoustic time is significantly larger the radial time, we can expect the method of extraction of a radial timescale to give reasonable results.

The radial timescales extracted by using the mode duplets , , and (all of them with ) are compared to the advective-acoustic and purely acoustic radial times in Figure 5. The timescales extracted from mode duplets are represented only on the radius interval in which both of the modes are unstable. The three extracted timescales are close to each other, which gives some confidence on their physical relevance. Furthermore, they roughly match the advective-acoustic time: within for moderate shock radius in the range (the focus of this article), and within for larger shock radii. By contrast, the extracted timescales are significantly longer than the radial acoustic time: the ratio varies between 3.5 for small shock radii and 6 for larger shock radii. This result clearly excludes the purely acoustic mechanism and argues in favour of the advective-acoustic cycle as the origin of SASI.

The difference at large shock radii between the advective-acoustic and the extracted timescale does not prevent the argument to be conclusive for several reasons. First, as compared to the large discrepancy with the purely acoustic time, this disagreement is quite slight. Second, it can easily be interpreted in the context of an advective-acoustic cycle by adjusting the effective coupling radius. Indeed the advective-acoustic time has been computed between the shock and the PNS surface, and as such should be considered as an upper limit. By contrast, it is impossible to reconcile in the same way the acoustic time with the measured timescales, because the acoustic time is too short. Finally we note that the instability mechanism at large shock radii was anyway not much debated, because the WKB analysis of Foglizzo et al. (2007) can be conclusively applied to higher frequency harmonics.

## 4 Computation of purely acoustic eigenmodes

In this section, a method to compute purely acoustic modes is described. In this method, one separates the advected perturbations just below the shock, and subtract them from the total perturbations at the shock. Let us emphasise that a WKB description of the acoustic waves is not needed, since we do not need to distinguish acoustic waves propagating up and down. As a result, the method is applicable to lower frequency modes. We stress that we do not make any assumption regarding the transverse propagation of acoustic waves. The WKB method is applicable when the wavelength of acoustic waves is shorter than the scale of the gradients: , while the separation of advected perturbations only requires that the advection wavelength be shorter than the gradient scale:

(10) |

For example for , the first three SASI modes have . We may thus expect that the advected wave can be extracted reasonably well for the second and third modes, but only marginally for the first mode.

Another assumption implicitly made here is that, just below the shock, the entropy is simply advected, which is true if the flow is adiabatic. The flow just below the shock is indeed usually rather close to adiabaticity. The quasi-adiabaticity of the postshock flow can be quantified by the dimensionless parameter (where is a dimensionless entropy defined in Foglizzo et al. (2007)), which ranges from for a shock radius at , to for .

In order to be able to subtract the entropy and vorticity perturbations (described by the variables and defined in Appendix A), one needs to determine the values of the other variables and associated to this advected wave. This can be done as in the WKB analysis of Foglizzo et al. (2007); Yamasaki & Foglizzo (2008), by requiring that these perturbations are advected:

(11) | |||||

(12) |

Using the differential system recalled in Appendix A to remove the radial derivatives (neglecting the non adiabatic terms in Equations (62)-(63) due the quasi-adiabatic assumption), Equations (11) and (12) can be expressed by the following linear system of two equations:

(13) | |||||

(14) |

Note that Equation (13) (signifying that is advected) is equivalent to requiring that the pressure perturbation vanishes (compare Equations (13) and (56)). Thus, when the advected wave is subtracted from the perturbations below the shock, the pressure perturbation is unchanged. This system of equations can be solved to obtain the following expression of the advected perturbations as a function of the amplitude of the entropy and vorticity perturbations (described by the variables and ):

(15) | |||||

(16) | |||||

In order to compute purely acoustic modes, the advected wave is then subtracted from the perturbations at the shock to define new boundary values and :

(17) | |||||

(18) |

where and are determined by the boundary conditions at the shock (Equations (65) and (65)). and are computed using Equations (15) and (16), and the values of and at the shock (Equations (67) and (68)). Using these modified boundary conditions, the acoustic modes can then be computed numerically as described for example in Foglizzo et al. (2007).

Figure 6 compares the growth rate of the most unstable acoustic modes (left panel) and SASI modes (right panel). The difference is striking: the purely acoustic modes are all stable for all shock radii in the explored range , while SASI modes have significant growth rates. This shows that the purely acoustic cycle cannot account for SASI and that the advected wave is essential to the instability. The most ”unstable” acoustic mode (actually the least damped) for each is the lowest frequency one. The discussion of the purely acoustic frequency spectrum is delayed to Section 5.2.2.

In order to quantify the effect of non-adiabaticity at the shock on the acoustic modes, we have used boundary conditions which either take into account or neglect the cooling processes in the computation of and . The associated change in the frequency and growth rate of the modes is barely noticeable (less than for any shock radius), confirming the validity of our quasi-adiabatic approximation.

The error due to the high frequency approximation can in turn be quantified in the following way. Noting that there is some arbitrariness in the choice of the variables that should be advected, one can require another variable to be advected and compare the result with the previous calculation. For example requiring the perturbation of radial velocity to be advected instead of , we get a different definition of the advected perturbations:

(19) | |||||

(20) |

where and are two parameters vanishing in the limit , and defined by:

(21) | |||||

(22) |

As expected, the two prescriptions for the advected perturbations are equivalent in the limit of vanishing and . These parameters are very small for the second and third acoustic modes: and respectively for (note that these two modes have a higher frequency than the second and third SASI modes as discussed in Section 5.2.2). The high frequency approximation is thus well justified and the two prescriptions are in reasonable agreement: the damping rates agree within for the second mode, and for the third mode. These two modes are enough to conclude that SASI is not a purely acoustic instability since the second and third SASI modes are unstable. The approximation is slightly less good for the fundamental acoustic mode: in this case . The damping rates in the two prescriptions differ by , but the modes are still definitely stable whatever the prescription.

## 5 Detailed analysis of the advective-acoustic and purely acoustic frequency spectra

In this section, we study in detail the frequency spectrum expected from advective-acoustic and purely acoustic cycles. In the subsection 5.1, we study a simple toy model (Foglizzo, 2009), which allows an analytical study of the frequency spectrum. In the following subsection 5.2.1, we then give an approximate generalisation to the more realistic model of SASI used in Sections 3 and 4. This analysis confirms the validity of the argument presented in Section 3. Furthermore it allows a direct comparison between analytical expectations and the frequency spectrum of SASI. This direct comparison leads to the same conclusion about the instability mechanism, in favour of the advective-acoustic cycle.

### 5.1 Planar toy model

#### 5.1.1 Presentation

We consider the toy model introduced by Foglizzo (2009). It consists of a planar shock lying above a zone of deceleration imposed by an external potential step. The potential step is localised within a length scale of , and is separated from the shock by a region of size where the flow is uniform. The gas is assumed to be perfect, adiabatic (with an adiabatic index ), and non-self-gravitating. The direction perpendicular to the shock is called , and referred to as vertical. Modes have a transverse structure along (parallel to the shock), which is referred to as the horizontal direction. The domain is a box of horizontal size . After a Fourier transform in time and horizontal direction , the eigenmodes are characterised by their complex frequency and the number of horizontal wavelengths in the width of the box. The horizontal wave number is then:

(23) |

The reader is referred to Foglizzo (2009) for more details about this toy model, the equations governing it, and the numerical method used to solve them. In the remaining of this subsection, the numerical calculations have used the parameters: , , , and .

Cycle constants and describing the advective-acoustic and purely acoustic cycle can be defined as in Section 2 (but here without any approximation). An eigenmode with a complex eigenfrequency must then verify Equation (5), and contains both cycles in some proportion. In order to separate the two cycles and study their properties independently from each other, we follow Foglizzo (2009) in considering advective-acoustic modes verifying:

(24) |

which corresponds to a situation where the acoustic wave propagating downward would be artificially removed. Similarly, purely acoustic modes follow:

(25) |

which signifies that the advected entropy-vorticity wave is ignored. This is equivalent to the method described in Section 4 to compute purely acoustic modes.

As in Foglizzo (2009), we define new cycle constants and that do not include the propagation of the waves between the shock and the potential step. These are defined as:

(26) | |||||

(27) |

where and are the vertical numbers of the advected wave, and of the acoustic waves propagating down (+) and up (-):

(28) | |||||

(29) | |||||

(30) |

#### 5.1.2 Eigenfrequencies of the advective-acoustic cycle

Equations (24) and (26) are combined to obtain the condition verified by an advective-acoustic mode:

(31) |

where is the duration of a one dimensional advective-acoustic cycle (, ):

(32) |

Equation 31 can be solved approximately in the regime where the growth rate is low () and where the acoustic wave is propagating (). In this case, the modulus of this complex relation gives access to the growth rate (Foglizzo, 2009):

(33) |

In a similar way, the frequency may be obtained from the argument of the same relation, which gives an information on the wave phase:

(34) |

where is the phase of , and describes the phase shift taking place during the two coupling processes. Equation (34) gives a second order polynomial in , which physical solution is:

(35) | |||||

where we defined two frequencies typical of the advective-acoustic cycle and of the horizontal propagation:

(36) | |||||

(37) |

Note that in principle depends on the frequency, therefore Equation (35) gives explicitly the frequency only if we know the value of . In the following of this subsection, we suppose which turns out to be a good approximation in this toy model. This is not always true in the spherical model as will be discussed in Section 5.2.2.

Figure 7 compares this analytical result with eigenmodes computed numerically. The agreement is very good both with advective-acoustic modes and with the full modes (i.e. the modes that verify Equation (5)), in the regime of propagative acoustic waves. The left panel shows that the approximations made are valid, and that there is no significant phase shift. The right panel shows that the frequency of the full modes is dictated by the advective-acoustic cycle, as could be expected from the results of Foglizzo (2009). One can nonetheless notice that the purely acoustic cycle has an influence on the growth rate of the modes (especially those close to marginal stability, (Foglizzo, 2009)), since the unstable modes are not exactly the same on the two panels. The origin of the stable mode at very low frequency, which does not coincide with any curve, is rather uncertain (it lies in the regime of evanescent acoustic waves, which is not described by Equation (35)).

From Figure 7 and Equation (35), one may notice that two successive modes sharing the same transverse structure are separated by a frequency difference that is close to the frequency of a vertical advective-acoustic cycle . This is exact in the case of one-dimensional modes, but only approximative for . Indeed a transverse advective-acoustic mode has a slightly higher frequency than the corresponding 1D mode (by a factor at most for a propagative mode). This higher frequency may seem counterintuitive at first sight, because the geometrical path along which the waves propagate during one cycle is longer. This longer path is the reason why the duration of the cycle increases with inclination in Equation (33) giving the growth rate. However this duration corresponds to the group velocity of acoustic waves, while the frequency is determined by a phase relation. The cycle time relevant to the frequency is thus given by the phase speed. The vertical phase speed of an acoustic wave increases with inclination, because the wave fronts get closer to the vertical direction. This explains the higher frequency of non radial modes.

#### 5.1.3 Eigenfrequencies of the purely acoustic cycle

The frequency of the purely acoustic modes can be obtained by a very similar method to that used for the advective-acoustic modes. For this purpose, one combines Equations (25) and (27) to obtain:

(38) |

where is the duration of a one-dimensional purely acoustic cycle ():

(39) |

By taking the argument of this complex equation, we obtain:

(40) |

from which one can deduce the following expression for the frequency of purely acoustic modes:

(41) |

where is the frequency associated to a vertical acoustic propagation:

(42) |

Equation (41) is analogous to the frequency of the acoustic modes of a rectangular box with periodic boundary conditions in and reflexive boundary conditions in (it is only slightly modified by advection):

(43) |

Thus as argued in Section 3, if the number of wavelengths in the transverse direction () increases, the frequency increases typically by corresponding to a propagation time in the transverse direction. Similarly if the number of wavelength in the vertical direction () increases, the frequency increases typically by corresponding to an acoustic propagation along the vertical direction.

Equation (41) (with ) is compared with the numerically computed eigenmodes in Figure 8. It agrees with the purely acoustic modes (all stable), showing here again that the approximation of small growth rate holds and that the phase shift is not very significant. By contrast, the frequency of the full modes do not match the purely acoustic prediction. This confirms that the unstable modes are not dominated by the purely acoustic cycle, and that the study of the frequency spectrum enables to discriminate between the two mechanisms.

Although the purely acoustic cycle does not determine the mode frequencies, it can be observed from Figure 8 that the acoustic timescale influences the modes stability. Indeed, the full modes in the vicinity of an acoustic branch are more likely to be unstable than the other modes. Furthermore, the most unstable modes (for a given ) are all close to the branch . These observations can be interpreted by two phenomena. First, a constructive interference between the purely acoustic and advective-acoustic cycles takes place when the mode frequency (determined by the advective-acoustic cycle) is close to that of a purely acoustic cycle (Foglizzo, 2002). This destabilising effect is more pronounced for modes close to marginal stability (Foglizzo, 2009). Second, the growth rate associated with the advective-acoustic cycle alone is maximum when the acoustic propagation is close to horizontality (Foglizzo, 2009), which happens at the frequency of the fundamental acoustic branch .

### 5.2 Spherical model

#### 5.2.1 Approximate generalisation of the analytical results

The description of waves propagation given in the previous subsection is exact in a uniform flow, and could be generalised to a non uniform flow in spherical geometry with the use of the WKB approximation. In this framework, the terms of the form (in Equations (26) and 27) should be replaced by terms of the following form: . However, as our goal is to discriminate the mechanism at work in the regime where the WKB approximation does not hold, we will not develop rigorously such a description. We simply give an approximate generalisation, from which one may expect the right order of magnitude but not an exact value. This is however sufficient to distinguish the two mechanisms as noted in Section 3.

The timescales associated to a radial advective-acoustic cycle, a radial purely acoustic cycle, and a transverse acoustic propagation have been defined for a spherical model in Section 2.3 (Equations (1)-(3)). One can define from these timescales, typical frequencies which are a direct generalisation of the frequencies of the planar model , , and :

(44) | |||||

(45) | |||||

(46) |

We propose to estimate the frequency spectra of advective-acoustic and purely acoustic cycles by directly introducing these typical frequencies in Equations (35) and (41) (for this purpose should be replaced by ). This procedure suffers from different sources of uncertainty: the formula itself is approximate, the phase shift is a priori unknown, and the typical frequencies , , and are uncertain respectively due the lack of knowledge of the coupling radius, the turning point of acoustic waves, and the radius at which it should be evaluated. When comparing to SASI modes in Section 5.2.2, we will use a range of values for the radius where and are evaluated in order to illustrate this uncertainty.

The above formulation can be used to assess the validity of the argument developed in Section 3. An inspection of Figures 7, 8, 9 and 10 confirms visually that the difference between two successive modes that share the same transverse structure corresponds approximately to a radial propagation frequency: for the advective-acoustic cycle, and for the purely acoustic cycle. This visual observation can be recovered from Equation (35) for the advective-acoustic cycle, if one assumes :

(47) | |||||

where we have defined , , and . Similarly for the purely acoustic cycle, Equation (41) with the assumption that leads to the following relation:

(48) | |||||

Thus we expect the timescale extracted from the eigenfrequencies in Section 3 to match one of the radial times (depending on the instability mechanism), for sufficiently large values of and if the two modes have the same phase shift. The deviation at low can be estimated from the second term in the right hand side of Equations (47) and (48). For example for the parameters , , , it is already less than for both the advective-acoustic and the purely acoustic modes (note that for the purely acoustic modes is the second mode). This confirms the validity of the argument given in Section 3, and is roughly consistent with the fact that the three extracted times differ by at most . We note that a difference in the phase shift between two successive modes could pollute the measure of the radial time, but it is not expected to consistently change this measure for different mode duplets.

#### 5.2.2 Comparison with SASI and purely acoustic eigenfrequencies

Figure 9 compares the frequency spectrum of purely acoustic modes with the analytical estimate. The method described in Section 4 gives two kinds of modes: the acoustic modes in which we are interested, as well as modes due to a secondary advective-acoustic cycle with the entropy being created by cooling processes acting on the acoustic wave. This last category of modes are clearly distinguished and significantly more stable than the purely acoustic modes, and therefore we do not show them in this discussion of the purely acoustic modes. Contrary to the case of the planar toy model, it is necessary to invoke a phase shift so that the analytical estimate agrees with the numerically obtained eigenfrequencies. We find that a phase shift of for all modes with , and for all modes with gives a reasonable match to the eigenspectrum (Figure 9). The origin of this phase shift remains obscure. Nevertheless, the comparison confirms that the analytical estimate accounts for the frequency of acoustic modes within the uncertainties discussed in the previous subsection and illustrated by the dashed lines in Figure 9.

Figure 10 compares the frequency of SASI eigenmodes with the predictions from the advective-acoustic cycle (left panel) and the purely acoustic cycle (right panel). The advective-acoustic cycle prediction agrees with the frequency of all the non radial modes, with a precision of ten to a few tens of percents. On the other hand, the purely acoustic cycle approximately predicts the frequency of the most unstable modes for each , but significantly overestimates that of other modes. This is consistent with the result of Section 3 that the acoustic radial time is significantly shorter than that extracted from the eigenspectrum. It is also consistent with the fact that the SASI eigenspectrum is much more densely populated than the acoustic eigenspectrum (compare Figure 9 and Figure 10, remembering that the frequency scale is different).

Similarly to the toy model, the most unstable SASI modes for each spherical harmonics index are close to the purely acoustic branch . This can again be interpreted by a higher efficiency of the advective-acoustic cycle, as well as a constructive interference with the purely acoustic mode. It is also probably linked with the observation of Fernández & Thompson (2009b) that the growth rate of a mode is maximum when the advection time coincides with the azimuthal acoustic time.

Finally, we observe that the (stable) radial modes are systematically shifted with respect the theoretical curve. This can be explained by a phase shift taking place during the advective-acoustic cycle, as shown by Fernández & Thompson (2009a). Indeed, a phase shift of then gives the following frequency spectrum for the modes:

(49) |

which is in much better agreement with the numerical eigenmodes.

## 6 Conclusion

Two new arguments have been presented which allow us to conclude that SASI cannot be explained by a purely acoustic mechanism, and should thus be interpreted as an advective-acoustic cycle. Contrary to previous works, these two arguments do not rely on a WKB analysis of the acoustic waves, and are therefore valid for realistic values of the shock radius.

First we proposed a method to extract a radial propagation timescale from the frequency difference between two successive SASI modes. This timescale agrees fairly well with the advective-acoustic time, but strongly disagrees with the radial acoustic time. This discards the purely acoustic cycle as the driving mechanism of SASI, and argues in favour of the advective-acoustic cycle. A detailed analysis of the frequency spectrum provided in Section 5 confirmed the validity of this argument.

Second, we described a new method to compute purely acoustic modes by artificially subtracting the advected wave created by the shock. All these acoustic modes were found to be stable, thus showing that the advected entropy-vorticity wave created by the oscillation of the shock plays a central role in the instability mechanism.

As a side product, we provided an analytical description of the frequency spectra of the advective-acoustic and purely acoustic cycles. This description was checked in the planar toy model of Foglizzo (2009) and the model of Blondin & Mezzacappa (2006). The similarity between the frequency spectra in these two models, supports the relevance of the toy model for core collapse supernovae. Finally, we note that unexplained phase shifts during the coupling processes were observed in the second model. This calls for a better understanding of the coupling process in a cooling layer.

## Acknowledgments

The authors are grateful to John Blondin for stimulating discussions. This work has been partially funded by the Vortexplosion project ANR-06-JCJC-0119 and the SN2NS project ANR-10-BLAN-0503. J.G. acknowledges support by STFC.

## Appendix A Differential system and boundary conditions with the new variables of Yamasaki & Foglizzo

Following Yamasaki & Foglizzo (2008), we adopt the variables defined by:

(50) | |||||

(51) | |||||

(52) | |||||

(53) |

where denotes the Eulerian perturbation. As compared to the variables used by Foglizzo et al. (2007), and are unchanged, the new variable is related to by , and is slightly changed to incorporate the perturbed heating . The usual physical quantities can be expressed as a function of these new variables as:

(54) | |||||

(55) | |||||

(56) | |||||

(57) |

The velocity in non radial directions has a different spatial structure and is expressed as a function of angular derivative of the spherical harmonics in the following way:

(58) | |||||

(59) |

The advantage of this new formulation is that the differential system describing the perturbations becomes particularly compact (note that it is strictly equivalent to that of Foglizzo et al. (2007)):

(60) | |||||

(61) | |||||

(62) | |||||

(63) |

where is defined as in Foglizzo et al. (2007) by:

(64) |

and the explicit expressions of and are given in Foglizzo et al. (2007) (Equations (A7) and (A8)).

The boundary conditions at the shock are obtained by imposing the Rankine-Hugoniot relations for the perturbed quantities, which are written as:

(65) | |||||

(66) | |||||

(67) | |||||

(68) |

where the subscripts ’’ and ’1’ refer to the values just below and above the shock, respectively; is the radial displacement of the shock surface.

## References

- Blondin & Mezzacappa (2006) Blondin J. M., Mezzacappa A., 2006, ApJ, 642, 401
- Blondin & Mezzacappa (2007) Blondin J. M., Mezzacappa A., 2007, Nature, 445, 58
- Blondin et al. (2003) Blondin J. M., Mezzacappa A., DeMarino C., 2003, ApJ, 584, 971
- Burrows et al. (2006) Burrows A., Livne E., Dessart L., Ott C. D., Murphy J., 2006, ApJ, 640, 878
- Deubner & Gough (1984) Deubner F.-L., Gough D., 1984, ARA&A, 22, 593
- Endeve et al. (2010) Endeve E., Cardall C. Y., Budiardja R. D., Mezzacappa A., 2010, ApJ, 713, 1219
- Fernández (2010) Fernández R., 2010, ApJ, 725, 1563
- Fernández & Thompson (2009a) Fernández R., Thompson C., 2009a, ApJ, 703, 1464
- Fernández & Thompson (2009b) Fernández R., Thompson C., 2009b, ApJ, 697, 1827
- Foglizzo (2001) Foglizzo T., 2001, A&A, 368, 311
- Foglizzo (2002) Foglizzo T., 2002, A&A, 392, 353
- Foglizzo (2009) Foglizzo T., 2009, ApJ, 694, 820
- Foglizzo et al. (2005) Foglizzo T., Galletti P., Ruffert M., 2005, A&A, 435, 397
- Foglizzo et al. (2007) Foglizzo T., Galletti P., Scheck L., Janka H.-T., 2007, ApJ, 654, 1006
- Foglizzo et al. (2009) Foglizzo T., Guilet J., Sato J., 2009, in M. Heydari-Malayeri, C. Reyl’E, & R. Samadi ed., SF2A-2009: Proceedings of the Annual meeting of the French Society of Astronomy and Astrophysics Asymmetric explosion of core collapse supernovae. pp 143–+
- Foglizzo et al. (2006) Foglizzo T., Scheck L., Janka H.-T., 2006, ApJ, 652, 1436
- Foglizzo & Tagger (2000) Foglizzo T., Tagger M., 2000, A&A, 363, 174
- Guilet & Foglizzo (2010) Guilet J., Foglizzo T., 2010, ApJ, 711, 99
- Guilet et al. (2010) Guilet J., Sato J., Foglizzo T., 2010, ApJ, 713, 1350
- Houck & Chevalier (1992) Houck J. C., Chevalier R. A., 1992, ApJ, 395, 592
- Iwakami et al. (2008) Iwakami W., Kotake K., Ohnishi N., Yamada S., Sawada K., 2008, ApJ, 678, 1207
- Kotake et al. (2007) Kotake K., Ohnishi N., Yamada S., 2007, ApJ, 655, 406
- Lai (2001) Lai D., 2001, in D. Blaschke, N. K. Glendenning, & A. Sedrakian ed., Physics of Neutron Star Interiors Vol. 578 of Lecture Notes in Physics, Berlin Springer Verlag, Neutron Star Kicks and Asymmetric Supernovae. pp 424–+
- Laming (2007) Laming J. M., 2007, ApJ, 659, 1449
- Laming (2008) Laming J. M., 2008, ApJ, 687, 1461
- Leonard et al. (2006) Leonard D. C., Filippenko A. V., Ganeshalingam M., Serduke F. J. D., Li W., Swift B. J., Gal-Yam A., Foley R. J., Fox D. B., Park S., Hoffman J. L., Wong D. S., 2006, Nature, 440, 505
- Liebendörfer et al. (2001) Liebendörfer M., Mezzacappa A., Thielemann F.-K., Messer O. E., Hix W. R., Bruenn S. W., 2001, Phys. Rev. D, 63, 103004
- Maeda et al. (2008) Maeda K., Kawabata K., Mazzali P. A., Tanaka M., Valenti S., Nomoto K., Hattori T., Deng J., Pian E., Taubenberger S., Iye M., Matheson T., Filippenko A. V., Aoki K., Kosugi G., Ohyama Y., Sasaki T., Takata T., 2008, Science, 319, 1220
- Marek & Janka (2009) Marek A., Janka H.-T., 2009, ApJ, 694, 664
- Marek et al. (2009) Marek A., Janka H.-T., Müller E., 2009, A&A, 496, 475
- Ohnishi et al. (2006) Ohnishi N., Kotake K., Yamada S., 2006, ApJ, 641, 1018
- Sato et al. (2009) Sato J., Foglizzo T., Fromang S., 2009, ApJ, 694, 833
- Scheck et al. (2008) Scheck L., Janka H.-T., Foglizzo T., Kifonidis K., 2008, A&A, 477, 931
- Scheck et al. (2006) Scheck L., Kifonidis K., Janka H.-T., Müller E., 2006, A&A, 457, 963
- Vishniac & Ryu (1989) Vishniac E. T., Ryu D., 1989, ApJ, 337, 917
- Wheeler et al. (2002) Wheeler J. C., Meier D. L., Wilson J. R., 2002, ApJ, 568, 807
- Yamasaki & Foglizzo (2008) Yamasaki T., Foglizzo T., 2008, ApJ, 679, 607
- Yamasaki & Yamada (2007) Yamasaki T., Yamada S., 2007, ApJ, 656, 1019