A Parametric Study of the SASI Comparing General Relativistic and Nonrelativistic Treatments

We present numerical results from a parameter study of the standing accretion shock instability (SASI), investigating the impact of general relativity (GR) on the dynamics. Using GR hydrodynamics with GR gravity, and nonrelativistic (NR) hydrodynamics with Newtonian gravity, in an idealized model setting, we vary the initial radius of the shock, and by varying its mass and radius in concert, the proto-neutron star compactness. We investigate four compactnesses expected in a post-bounce core-collapse supernova (CCSN). We find that GR leads to a longer SASI oscillation period, with ratios between the GR and NR cases as large as 1.29 for the highest-compactness suite. We also find that GR leads to a slower SASI growth rate, with ratios between the GR and NR cases as low as 0.47 for the highest-compactness suite. We discuss implications of our results for CCSN simulations.


INTRODUCTION
Since the discovery of the standing accretion shock instability (SASI; Blondin et al. 2003), which many twoand three-dimensional simulations performed to date demonstrate becomes manifest during a core-collapse supernova (CCSN) in the post-shock accretion flow onto the proto-neutron star (PNS), groups have made efforts to understand its physical origin and its effects on the supernova itself.The SASI is characterized in two dimensions (2D) by a large-scale "sloshing" of the shocked fluid, and in three dimensions (3D) by additional spiral modes (Blondin & Mezzacappa 2007).It is now generally accepted that turbulent neutrinodriven convection plays a major role in re-energizing the stalled shock (e.g., see Burrows et al. 2012;Hanke et al. 2013;Murphy et al. 2013;Couch & Ott 2015;Melson et al. 2015b;Lentz et al. 2015;Abdikamalov et al. 2015;Radice et al. 2015;Melson et al. 2015b,a;Roberts et al. 2016;Müller et al. 2017;Summa et al. 2018;Radice et al. 2018;O'Connor & Couch 2018a;Vartanyan et al. 2019;Müller et al. 2019;Burrows et al. 2019;Yoshida et al. 2019;Powell & Müller 2020;Stockinger et al. 2020;Müller & Varma 2020;Vartanyan et al. 2022;Nakamura et al. 2022;Matsumoto et al. 2022).The same simulations that led to the above conclusion also generally exhibit the SASI, with outcomes ranging from convectiondominated flows to SASI-dominated flows, and flows where neither dominates.Strong SASI activity and, in some cases, SASI-aided explosions have been reported in, for example, the three-dimensional simulations of Summa et al. (2018), O'Connor & Couch (2018a), and Matsumoto et al. (2022).A more precise determination of the relative role played by these two instabilities in the explosion mechanism, on a case-by-case basis (i.e., for different progenitor characteristics; e.g., see Scheck et al. 2008;Hanke et al. 2012Hanke et al. , 2013;;Couch & O'Connor 2014;Fernández et al. 2014;Melson et al. 2015a;Abdikamalov et al. 2015;Fernández 2015), will require advances in current three-dimensional models to include full general relativity, rotation, magnetic fields, and the requisite neutrino interaction physics with realistic spectral neutrino kinetics, all at high spatial resolution.It is also important to note that, while convection-dominated and SASI-dominated scenarios may lie at the extremes of what is possible, it is not necessary for one or the other instability to be dominant to play an important role -specifically, for the more complex cases where neither dominates, it would be very difficult to determine precisely the relative contribution from these two instabilities.
Several studies have concluded that the SASI is an advective-acoustic instability, in which vortical waves generated at the shock advect to the surface of the PNS, which in turn generate acoustic waves that propagate back to the shock and further perturb it (Foglizzo et al. 2006(Foglizzo et al. , 2007;;Yamasaki & Yamada 2007;Laming 2007Laming , 2008;;Foglizzo 2009;Guilet & Foglizzo 2012).This perturbation generates more, stronger, vortical waves, which advect to the PNS surface, thus creating a feedback loop that drives the instability.An alternative explanation for the SASI is the purely acoustic mechanism, in which acoustic perturbations just below the shock travel around the post-shock region and constructively interfere with each other, generating stronger acoustic perturbations and thereby feeding the instability (Blondin & Mezzacappa 2006).A recent study (Walk et al. 2023) suggests that the acoustic mechanism may play a particularly important role in the SASI when rotation is included, implying that the origins of the SASI excitation may depend on conditions between the shock and the PNS.
Other numerical studies focus on particular aspects of this instability, such as the hydrodynamics of the SASI (Ohnishi et al. 2006;Sato et al. 2009;Iwakami et al. 2014), spiral modes (Blondin & Shaw 2007;Iwakami et al. 2008;Fernández 2010), the spin-up of the possible remnant pulsar (Blondin & Mezzacappa 2007), the effect of nuclear dissociation (Fernández & Thompson 2009), saturation of the instability (Guilet et al. 2010), the generation and amplification of magnetic fields (En-deve et al. 2010(En-deve et al. , 2012)), the relative importance of the SASI and convection in CCSNe (Cardall & Budiardja 2015), the generation of, and impact on, gravitational waves by the SASI (Kotake et al. 2007(Kotake et al. , 2009;;Kuroda et al. 2016;Andresen 2017;Kuroda et al. 2017;Andresen et al. 2017;O'Connor & Couch 2018a;Hayama et al. 2018;Andresen et al. 2019;Mezzacappa et al. 2020Mezzacappa et al. , 2023;;Drago et al. 2023), and the effects of rotation (Yamasaki & Yamada 2005;Yamasaki & Foglizzo 2008;Walk et al. 2023;Buellet et al. 2023).Some of these studies included sophisticated microphysics, such as realistic equations of state (EoSs) and neutrino transport; however, with the exception of Kuroda et al. (2017), none of these studies solved the general relativistic hydrodynamics (GRHD) equations, instead solving their nonrelativistic (NRHD) counterparts, some with an approximate relativistic gravitational potential.It has been demonstrated that GR effects are crucial to include in CCSN simulations (Bruenn et al. 2001;Müller et al. 2012;Lentz et al. 2012;O'Connor & Couch 2018b), yet the SASI itself has not been fully investigated in the GR regime.A recent paper (Kundu & Coughlin 2022) does analyze steady state accretion through a stationary shock onto compact objects in a Schwarzchild geometry and compares with Newtonian solutions, and posits that GR may have a non-negligible impact on the SASI.They find that, for conditions expected in exploding CCSNe, the freefall speed is of order v ∼ 0.2 c (with c the speed of light), and the differences between the GR and NR solutions are of order 10%.For conditions expected in failed CCSNe (i.e., supernovae where the shock is not revived, in which case the freefall speed can be v ∼ 0.66 c), the differences can be larger.
The timescales that likely influence the SASI depend on signal speeds associated with advective and acoustic modes in the region between the shock and the PNS surface (Blondin & Mezzacappa 2006;Foglizzo et al. 2007;Müller 2020).Motivated in part by Dunham et al. (2020) and Kundu & Coughlin (2022), we expect SASI simulations to behave differently depending on whether or not the treatment of the hydrodynamics and gravity are general relativistic.Indeed, Dunham et al. (2020) presented GR steady-state solutions and presented preliminary results from GR SASI simulations, but did not compare results from NR and GR simulations.Specifically, we expect both advective and acoustic modes to be influenced by the different post-shock structure in the GR case relative to the NR case.
This leads to our main science question: How does a general relativistic treatment of hydrodynamics and gravity affect the oscillation period and growth rate of the SASI?To begin to address this, we present the first comparison of the SASI in both a nonrelativistic and a general relativistic framework, using an idealized model with four compactnesses, meant to span the range of conditions expected in CCSN simulations.We focus our attention on the linear regime and characterize the SASI by its growth rate and oscillation period, as was done in Blondin & Mezzacappa (2006).To capture both the linear regime of the SASI and its transition to the nonlinear regime, we perform our assessment via seven axisymmetric numerical simulations using GRHD and GR gravity, with the PNS represented by a point mass and gravity encoded in a Schwarzchild spacetime metric.To better assess the impact of GR, we also perform simulations using the same parameter sets but with NRHD and Newtonian gravity, again with the PNS represented by a point mass, but in this case gravity is encoded in the Newtonian potential.We will often use "NR" to refer to the case of Newtonian gravity and NRHD.
We use a system of units in which c = G = 1 and also make use of the Einstein summation convention, with Greek indices running from 0 to 3 and Latin indices running from 1 to 3.

Relativistic Gravity: Conformally-Flat Condition
We use the 3+1 decomposition of spacetime (see, e.g., Banyuls et al. 1997;Gourgoulhon 2012;Rezzolla & Zanotti 2013, for details), which, in the coordinate system x = t, x i , introduces four degrees of freedom: the lapse function, α (x), and the three components of the shift vector, β i (x).We further specialize to the conformallyflat condition (CFC, Wilson et al. 1996), effectively neglecting the impact of gravitational waves on the dynamics.This is a valid approximation when the CCSN progenitor is nonrotating (Dimmelmeier et al. 2005), as is the case for our simulations.The CFC forces the components of the spatial three-metric, γ ij (x), to take the form where ψ (x) is the conformal factor and the γ ij are the components of a time-independent, flat-space metric.
We choose an isotropic spherical-polar coordinate system, as it is appropriate to our problem and is consistent with the CFC; the flat-space metric is and the lapse function, conformal factor, and shift vector take the form given in Baumgarte & Shapiro (2010), where r > R Sc is the isotropic radial coordinate measured from the origin and R Sc := M/2 is the Schwarzchild radius in isotropic coordinates for an object of mass M .The line element under a 3+1 decomposition in isotropic coordinates takes the form We note here that the proper radius, R (r), corresponding to the coordinate radius, r > R Sc , is defined by where we used Eqs.(1-2) with ψ given by Eq. ( 4).Under the CFC, the square root of the determinant of the spatial three-metric is Our choice of isotropic coordinates is consistent with our implementation of general relativistic hydrodynamics in the CFC approximation.Our comparison of the Newtonian and general relativistic results is conducted within this gauge.(In practice, the hydrodynamics equations in the Newtonian limit in the isotropic gauge are identical in form to the standard Newtonian hydrodynamics equations, and we were able to run our NR simulations with our Newtonian hydrodynamics code.)

Relativistic Hydrodynamics
We solve the relativistic hydrodynamics equations of a perfect fluid (i.e., no viscosity or heat transfer) in the Valencia formulation (Banyuls et al. 1997;Rezzolla & Zanotti 2013), in which they take the form of a system of hyperbolic conservation laws with sources.Under our assumption of a stationary spacetime, the equations can be written as where U := U (t, r, θ, φ) is the vector of evolved fluid fields, is the vector of fluxes of those fields in the i-th spatial dimension, and S := S (U ) is the vector of sources, where D is the conserved rest-mass density, S j is the component of the Eulerian momentum density in the j-th spatial dimension, and τ := E − D, with E the Eulerian energy density.The component of the fluid three-velocity in the j-th spatial dimension is denoted by v j , and is the Lorentz factor of the fluid, both as measured by an Eulerian observer.The relativistic specific enthalpy as measured by an observer comoving with the fluid; i.e., a comoving observer, is h := 1 + (e + p) /ρ, where ρ is the baryon mass density, e is the internal energy density, and p is the thermal pressure, all measured by a comoving observer.Finally, i.e., γ ik γ kj = δ i j .See Rezzolla & Zanotti (2013) for more details.
We close the hydrodynamics equations with an ideal EoS, p (e) = (Γ − 1) e, where Γ ∈ (1, 2] is the ratio of specific heats.For this study, we set Γ = 4/3.We further assume the EoS is that of a polytrope; i.e., where K > 0 is the polytropic constant, whose logarithm can be considered a proxy for the entropy, S; i.e., S ∝ log p/ρ Γ .The constant K takes different values on either side of a shock, in accordance with physically admissible solutions.Eq. ( 14) is consistent with Eq. ( 13) through the first law of thermodynamics for an isentropic fluid.

Nonrelativistic Hydrodynamics
Under the 3+1 formalism of GR and the CFC, the effect of gravity is encoded in the metric via the lapse function, the conformal factor, and the shift vector, whereas with Newtonian gravity, the metric is that of flat space and the effect of gravity is encoded in the Newtonian gravitational potential, Φ, Of course, the NRHD equations can be recovered from the GRHD equations by taking appropriate limits; i.e., v i v i ≪ 1, p, e ≪ ρ, and Φ ≪ 1, and setting α ≈ α nr := 1 + Φ and ψ = 1 − Φ/2.In the case of Newtonian gravity and NRHD, we solve where and where P ij := ρ v i v j + p γ ij and we assume Φ is due only to the point source PNS,

STEADY-STATE ACCRETION SHOCKS
We take initial conditions for our simulations from steady-state solutions to Eq. (9) (GR) and Eq. ( 16) (NR).To determine the steady-state solutions, we assume the fluid distribution is spherically symmetric and time-independent.Following Blondin et al. (2003), we consider a stationary accretion shock located at r = R sh with a PNS mass M pns , PNS radius R pns , and a constant mass accretion rate Ṁ > 0. We assume a polytropic constant ahead of the shock, K = 2 × 10 14 erg cm −3 / g cm −3 4/3 , chosen so that the pre-shock flow is highly supersonic (all of our models have a pre-shock Mach number greater than 15).Given that our steady-state solutions have constant entropy between the PNS surface and the shock, they are convectively stable.This enables us to isolate the SASI and study its development.

Relativistic Steady-State Solutions
Focusing on the equation for D, we find (temporarily Manipulation of the equations for D and τ in Eq. ( 9) yields the relativistic Bernoulli equation, where B > 0 is the relativistic Bernoulli constant.At spatial infinity, the fluid is assumed to be at rest and the spacetime curvature negligible, so α ∞ = W ∞ = 1.Further, at spatial infinity, we assume the fluid to be cold; i.e., (e + p) /ρ ≪ 1, so that h ∞ = 1 and B ∞ = 1.Since B is a constant, B = 1 everywhere.Given K, Eqs. ( 14), ( 21), and ( 22) (with B = 1) form a system of three equations in the three unknowns, ρ, v, and p. From initial guesses v 0 = − 2 M pns /r, ρ 0 = − Ṁ / 4π r 2 v 0 , and p 0 = K ρ Γ 0 , the first two of which are obtained from the Newtonian approximation at a distance r > R sh for highly supersonic flow, we define dimensionless variables u 1 := ρ/ρ 0 , u 2 := v/v 0 , and u 3 := p/p 0 .These are substituted into the system of equations, which are then solved with a Newton-Raphson algorithm to determine the state of the fluid everywhere ahead of the shock.To join the pre-and post-shock states of the fluid at r = R sh , we apply the relativistic Rankine-Hugoniot jump conditions (i.e., the Taub jump conditions, Taub 1948) to obtain ρ, v, and p just below the shock.Once the state of the fluid just below the shock is found, the polytropic constant for the post-shock fluid is computed with Eq. ( 14) and the same system of equations is solved for the state of the fluid everywhere below the shock.

Nonrelativistic Steady-State Solutions
The steady-state solution method for the nonrelativistic case (taken from Blondin et al. (2003)) follows a similar procedure as the relativistic case, except we begin from the NR equations for mass density and energy density, where From these, and making the same assumptions as in the relativistic case, we arrive at a system of two equations for the three unknowns, ρ, v ≡ v r , and p, with h nr = h − 1 = (e + p) /ρ the nonrelativistic specific enthalpy and B nr the nonrelativistic Bernoulli constant.Following Blondin et al. (2003), we set B nr = 0.As in the GR case, we close this system with Eq. ( 14).

Comparison of NR and GR Steady-State Solutions
For the two lower-compactness cases (see Eq. ( 33)) (e.g., see Bruenn et al. 2013;Melson et al. 2015b;Burrows et al. 2020), Figure 1 shows steady-state accretion shock solutions as functions of coordinate distance r. (Note that in Figure 1 and Figure 2 below, for additional information, the steady-state solutions are plotted from the shock all the way down to r = 3 km.For our numerical simulations, the inner boundary is placed at R pns ∈ {20, 40} km, which then defines the compactness of our models.) In general, the magnitudes of the density, velocity, and pressure just below the shock decrease as the shock radius increases (e.g., see Eqs. (1-3) in Blondin et al. 2003).From the top-right panel of Figure 1, it can be seen that the velocities in the GR and NR cases agree well near the shock and deviate from each other for smaller radii, with the velocity being smaller in the GR case than in the NR case.The top-left and bottom-left panels show that the densities and pressures in the GR case are larger than their NR counterparts at smaller radii.The slope of the NR density profile matches expectations of ρ (r) ∝ r −3 (Blondin et al. 2003), but the GR density profile deviates noticeably from this as the inner-boundary is approached.From the bottom-right panel, it can be seen that the lapse function and its Newtonian approximation begin to deviate from each other near r = 40 km, with the degree of deviation increasing for smaller radii.
For the two higher-compactness cases (e.g., see Liebendörfer et al. 2001;Walk et al. 2020), Figure 2 shows steady-state accretion shock solutions as functions of coordinate distance r.The profiles show the same trends as those in Figure 1, although in this case the trends are more pronounced.One notable difference is that the location of the largest deviation in the velocity between the NR and GR cases occurs further in, near r = 12 km.Another notable feature of both Figures 1 and 2 is that the fluid velocity in the GR case is consistently slower than that in the NR case, in agreement with Kundu & Coughlin (2022).
We compare our numerical results with an estimate provided by Müller (2020), which provides an analytic estimate of the oscillation period of the SASI, T aa , based on the assumption that the SASI is an advectiveacoustic cycle, in which a fluid parcel advects from the shock to the PNS surface in time τ ad , which generates acoustic waves that propagate from the PNS surface to the shock in time τ ac .We modify that formula to include the effects of GR by including the metric factor, which involves the conformal factor and which converts the radial coordinate increment to the proper radial distance increment, and by replacing the nonrelativistic signal speeds with their relativistic counterparts, where λ r 0 and λ r + are the radial signal speeds of matter and acoustic waves, respectively.Using our metric, the signal speeds are (Rezzolla & Zanotti 2013) where c s is the sound-speed and where the second equality in each expression is the nonrelativistic limit.For 1D problems, this expression depends only on the steadystate values of c s and v r .We also compare our results with an estimate based on the assumption that the SASI is a purely acoustic phenomenon.We define a time, T ac , as the time taken by an acoustic perturbation to circum-navigate the post-shock cavity at a characteristic radius R ac , assuming v θ ≪ c s γ 22 , is the acoustic wave-speed in the θ-dimension (Rezzolla & Zanotti 2013) and h θ is the scale factor in the θdimension.The expression for T ac is inspired by the Lamb frequency, which relates the frequency of an acoustic wave to the radius at which it turns around (Hansen et al. 2004).We delay the specification of R ac until Section 5.
In Figure 3, we plot λ r,gr + /λ r,nr + , λ θ,gr + /λ θ,nr + , and λ r,gr 0 /λ r,nr 0 as functions of η, defined as for all of our models.In all cases, the signal speeds are slower in the GR case.This difference is accentuated in the higher-compactness models and for smaller radii.
The growth rate of the SASI depends on the steadystate conditions below the shock; however, obtaining an analytic estimate for this is difficult, and although there have been efforts to explain the physics governing the growth rate assuming nonrelativistic models (e.g., Blondin & Mezzacappa 2006;Foglizzo et al. 2007;Laming 2007Laming , 2008;;Foglizzo 2009;Guilet & Foglizzo 2012), an analytic estimate remains an open question and no estimate exists for a GR model.Here, we aim to compare the NR and GR growth rates with numerical simulations.
We emphasize that our goal is to characterize the SASI in terms of its period and growth rate, and to compare their GR and NR values.Determining its physical origin, whether it be advective-acoustic or purely acoustic, is beyond the scope of this study.The rough estimates provided by Eq. ( 27) and Eq. ( 30) are merely intended as points of reference for our numerically determined values.

SIMULATION CODE AND SETUP
We perform our simulations with thornado1 , an opensource code under development, aiming to simulate CC-SNe.thornado uses high-order discontinuous Galerkin (DG) methods to discretize space and strong-stabilitypreserving Runge-Kutta (SSP-RK) methods to evolve in time.For details on the implementation in the nonrelativistic case, see Endeve et al. (2019) and Pochik et al. (2021); for the relativistic case, see Dunham et al. (2020) and Dunham et al. (2024) (in prep.).All of our simulations use the HLL Riemann solver (Harten et al. 1983), a quadratic polynomial representation (per dimension) of the solution in each element, and third-order SSP-RK methods for time integration (Gottlieb et al. 2001) with a timestep ∆t := C cfl 1 d(2 k+1) min i∈{1,...,d} ∆x i / λ i , where C cfl = 0.5 is the CFL number, k the degree of the approximating polynomial (in our case, k = 2), ∆x i the mesh width in the i-th dimension, λ i the fastest signal speed in the i-th dimension, and d the number of spatial dimensions (Cockburn & Shu 2001).
Two important aspects of successful implementations of the DG method are mitigating spurious oscillations and enforcing physical realizability of the polynomial approximation of the solution.To mitigate oscillations, thornado uses the minmod limiter and applies it to the characteristic fields (see, e.g., Shu 1987;Pochik et al. 2021).For the interested reader, we set the β tvd parameter for the limiter, defined in Pochik et al. (2021), to 1.75 for all runs.thornado also uses the troubledcell indicator described in Fu & Shu (2017) to decide on which elements to apply the minmod limiter; for the threshold of that indicator we use the value 5 × 10 −3 .To enforce physical realizability of the solution in the NR case, thornado uses the positivity limiter described in Zhang & Shu (2010), and in the GR case, uses the limiter described in Qin et al. (2016); for the thresholds of both limiters we use the value 10 −13 .
The hydrodynamics in thornado has recently been coupled to AMReX2 , an open-source software library for block-structured adaptive mesh refinement and parallel computation with MPI (Zhang et al. 2019); however, our simulations are all performed on a uni-level mesh.
Our computational domain, D, is defined for all models as D := [R pns , 1.5 R sh (t = 0)] × [0, π].The radial extent allows us to determine whether or not the SASI has become nonlinear, which we define to be when any radial coordinate of the shock exceeds 10% of the initial shock radius.All simulations are evolved sufficiently long to achieve ten full cycles of the SASI.In some cases, the shock exceeds our threshold of nonlinearity before completing ten full cycles; in those cases, we only use data from the linear regime.
The PNS is treated as a fixed, spherically symmetric mass in order to maintain a steady-state, and we ignore the self-gravity of the fluid.Because the largest accretion rate we consider is 0.5 M ⊙ s −1 and because this lasts for a maximum of 550 ms in our 2D models, the most mass that would accrete onto the PNS is 0.275 M ⊙ and therefore would provide a sub-dominant contribution to M pns in our simulations.
We consider models with three free parameters: the mass of the PNS, M pns , the radius of the PNS, R pns , and the initial radius of the shock, R sh .We also varied the mass accretion rate, Ṁ , but found the oscillation periods and growth rates to be insensitive to this parameter, and we do not discuss these models further; all following discussion is for models with an accretion rate Ṁ = 0.3 M ⊙ s −1 .Our choice of parameters is motivated by the physical scale of CCSNe; the ranges of our parameter space are informed by models from Liebendörfer et al. ( 2001 2020), and can be found in Table 1.We also classify our simulations by their compactness, ξ, which we define as (O'Connor & Ott 2011), In our model naming convention, we first list whether the model used GR or NR along with the dimensionality (1D or 2D), followed by the mass of the PNS in Solar masses, the radius of the PNS in kilometers, and lastly the shock radius in kilometers; e.g., the 2D GR model with M pns = 1.4 M ⊙ , R pns = 40 km, and R sh = 150 km is named GR2D M1.4 Rpns040 Rsh1.50e2.If we compare an NR model with a GR model having the same parameters we may drop that specification from the model name.If no confusion will arise, we may also drop the dimensionality.
The inner radial boundary corresponds to the surface of the PNS.To determine appropriate inner-boundary conditions in the GR case, we assume D and τ follow power laws in radius; from the initial conditions, we extrapolate D and τ in radius using a least-squares method with data from the innermost five elements on the grid to determine the appropriate exponents.The radial momentum density interior to R pns is kept fixed to its initial value.We leave the outer radial boundary values fixed to their initial values for all fields.In the polar direction, we use reflecting boundary conditions at both poles.For the inner-boundary conditions in the nonrelativistic case, we assume ρ (r) ∝ r −3 and E (r) ∝ r −4 , the latter of which follows from our assumption of ρ ∝ r −3 , Eq. ( 13), Eq. ( 14), the assumption of Γ = 4/3, and the assumption of a small velocity at the PNS surface.
For the ξ ∈ {0.4,0.7} models we enforce a radial resolution of 0.5 km per element for all runs, which we found necessary for the shock in an unperturbed model to not deviate by more than 1% over 100 advection times.This is shown in the top panel of Figure 4, which plots the relative deviation of the shock radius from its initial position as a function of t/τ ad for runs with different radial resolutions for model GR1D M1.4 Rpns040 Rsh1.20e2.These results suggest that the steady-state is not maintained if the radial resolution is too coarse; e.g., greater than about 1 km.
For the ξ ∈ {1.8, 2.8} models we enforce a radial resolution of 0.25 km per element for all runs in order to maintain the same radial resolution of the pressure scale height, p |dp/dr| −1 , as in the ξ ∈ {0.4,0.7} models, while also ensuring that the shock does not deviate from its initial location by more than 1%.This can be seen in the bottom panel of Figure 4, which plots the same quantity as the top panel, but for model GR1D M2.8 Rpns020 Rsh6.00e1.
To verify that our chosen angular resolution of 64 elements (∼2.8 • ) is sufficient to resolve the angular variations of the fluid, we run two additional simulations of model NR2D M2.8 Rpns040 Rsh1.20e2, one with 128 angular elements and one with 256 angular elements.From those runs, we extract the best-fit values for the growth rates and oscillation periods (see Section 5) and find them to not significantly differ from those of the 64-angular-element run.
Our simulations are initialized with the steady-state solutions discussed in Section 3, and we take extra care to minimize initial transients.The initial conditions we obtain come from solving Eqs. ( 14), (21), and ( 22) in the GR case, and Eqs. ( 14), (25), and (26) in the NR case, which are not exact solutions of our discretized equations, so transients will be present when simulations are initialized with these solutions.To mitigate the effects of the transients, the fields are set up in 1D with the method described above and then evolved for 100 advection times, which was experimentally determined to be of sufficient duration to quell any transients.We verify that the system has achieved a steadystate by plotting, in Figure 5, the maximum, at each snapshot, of the absolute value of the normalized time derivative of the mass density versus t/τ ad for model GR1D M1.4 Rpns040 Rsh1.20e2 (top panel) and model GR1D M2.8 Rpns020 Rsh6.00e1 (bottom panel); other  models exhibit similar behavior.For example, in the top panel, it can be seen that the model settles down after approximately 35 advection times, followed by two slight increases near 45 and 50 advection times, then settles down until we end the simulation after 100 advection times.We attribute the two slight increases to limiters activating when the shock crosses an element boundary.
The relaxed 1D data is mapped to 2D, a perturbation to the pressure is applied (see below), and the system is evolved.The initial relaxation in 1D removes numerical noise and allows for a smaller perturbation amplitude, leading to a longer-lasting linear regime, which makes for a cleaner signal.We seed the instability by imposing a pressure perturbation onto the steady-state flow below the shock, with p (r) the steady-state pressure at radius r, and where δp (r, θ) /p (r) is defined in a scale-independent manner as with η defined in Eq. ( 32), and where η c = 0.75 and σ = 0.05.The perturbation is not allowed to extend into the pre-shock flow.We opted to perturb the postshock flow as opposed to the pre-shock flow after having tried various pre-shock and post-shock perturbations and finding that pre-shock perturbations generate noise when they cross the shock front, thus creating a noisy signal and making it more difficult to extract quantities of interest.Similarly, we choose a Gaussian profile because the smoothness of the profile generates less noise than, e.g., a top-hat profile.The cos θ factor is meant to excite an ℓ = 1 Legendre mode of the SASI.While this perturbation method does not exactly mimic the hydrodynamics inside a CCSN, it is sufficient to study the SASI in the linear regime.

RESULTS AND DISCUSSION
Here we discuss our analysis methods and compare the SASI growth rates and oscillation periods for our 7 models in GR and NR.

Analysis Methods
To extract SASI growth rates from our simulations, we follow Blondin & Mezzacappa (2006) and expand a quantity affected by the perturbed flow, A (t, r, θ), in Legendre polynomials, where we normalize the P ℓ such that with δ ℓℓ ′ the Kronecker delta function.Then, After experimenting with several quantities, we decided to use the quantity proposed by Scheck et al. (2008), where v θ is the fluid velocity in the polar direction as measured by an Eulerian observer, having units of rad s −1 .With G ℓ , we compute the power in the ℓ-th Legendre mode, H ℓ (t), by integrating over a shell below the shock, bounded from below by r a = 0.8 R sh and from above by r b = 0.9 R sh , For the Newtonian runs, ψ = 1.We experimented with different values of r a and r b and found that a thin shell just below the shock gave the cleanest signal.
To extract the SASI growth rate and oscillation period from H 1 , we begin by fitting the simulation data to the function (Blondin & Mezzacappa 2006) where ω is the growth rate of the SASI, T is the oscillation period of the SASI, F is an amplitude, and δ is a phase offset.We fit the data to the model using the Levenberg-Marquardt nonlinear least squares method (e.g., see Moré 1978), provided by scipy's curve fit function, which also provides an estimate on the uncertainty of the fit via the diagonal entries of the covariance matrix; we use this to define the uncertainty in the growth rate.The temporal extent over which we perform the fit is defined to begin after one SASI oscillation and to end after seven SASI oscillations, where, for simplicity, we use Eq. ( 27) to define the period of one SASI oscillation.
The period we report is obtained by performing a Fourier analysis: We integrate the lateral flux in the radial direction, from r a to r b (defined as in the computation of the H ℓ ), integrate the result over ) , and then take the Fourier transform of that result using the fft tool from scipy.From the result, F r θ ( T ), we define the period of the SASI as the value of T corresponding to the peak of the Fourier amplitudes, and we define the uncertainty in the period as the full width at half maximum (FWHM) of the Fourier amplitudes.The FFT is computed over the same time interval as the aforementioned fit.(We do not use the period returned from curve fit because, for some models, pollution from higher-order modes spoils the ability of our fitting function to capture the ℓ = 1 mode.)

Overall Trends
Here we discuss trends that appear across our models.We summarize our results in Table 2, which lists the model, the best-fit oscillation period and uncertainty, the best-fit growth rate and uncertainty, the product of the best-fit growth rate and best-fit oscillation period, the period assuming an advective-acoustic mechanism, T aa (estimated using Eq. ( 27)), and the period assuming a purely acoustic mechanism, T ac (estimated using Eq. ( 30)).When using Eq. ( 30), we choose the characteristic radius R ac to be 0.85 R sh , which is the midpoint of the shell in which we compute the power.
As a first general trend, we see that, for a given PNS radius, the oscillation period increases as the shock radius increases, as seen clearly in the second column of Table 2.This is expected, for as the shock radius increases, the waves supported by the fluid must traverse a larger region, therefore each cycle will take longer.As a second general trend, we observe that, for a given PNS radius, the growth rate decreases as the shock radius increases.

Impact of GR
Next we discuss the impact of GR on the oscillation period and the growth rate.

Oscillation Period
First we discuss how the oscillation period varies with PNS compactness and shock radius.In Figure 6, we plot the amplitudes of the Fourier transform of Eq. ( 42), F r θ , versus T for four models, where T is defined as the inverse of the frequency determined by the FFT.No windowing was applied when computing the FFT of the signal.We see that the difference in the optimal period, T (defined as the T associated with the largest Fourier amplitude), between NR and GR increases with increasing compactness, with the GR period consistently longer than the NR period.This can be explained by differences in the structure of the post-shock solutions; in particular, the signal speeds, which are shown in Figure 3.Both the radial acoustic and advective signal speeds, as well as the angular acoustic signal speed, are consistently smaller in GR.Because of this, the period is longer for a given model when GR is used, regardless of whether the SASI is governed by an advective-acoustic or a purely acoustic cycle.
In Figure 7, we plot the SASI period, as provided by the Fourier analysis, versus initial shock radius for all of our models.The NR and GR models follow the same general trends.For ξ = 0.4, 0.7, 1.8, and 2.8, we find Note-Oscillation periods, growth rates, and their uncertainties for all fourteen models having the same accretion rate of 0.3 M⊙ s −1 .The uncertainties for the growth rates are defined as the square roots of the diagonal entries of the covariance matrix corresponding to the growth rate.The uncertainies for the oscillation period are defined as the full-width half-maximum values of the Fourier amplitudes (see Figure 6).The fourth column shows the product of the best-fit growth rate multiplied by the best-fit oscillation period.The fifth column shows the estimate for the period assuming an advective-acoustic origin of the SASI (Eq.( 27)), and the sixth column shows the estimate for the period assuming a purely acoustic origin of the SASI (Eq.( 30)), where we use Rac = 0.85 R sh , the midpoint of the shell in which we compute the power (see Section 5).

Growth Rate
Next we discuss how the growth rate varies with PNS compactness and shock radius.In Figure 8, we plot the power in the first Legendre mode, H 1 , during the linear regime, as a function of time in milliseconds for the same models as in Figure 6.For all models displayed here, the shock deviates from spherical symmetry by less than 10%, which we consider sufficient for characterizing the evolution as being in the linear regime.
We see a clear trend in the growth rate: the GR models display a slower SASI growth rate when compared to their NR counterparts.In Figure 8, it can be seen that the growth rate for the ξ = 0.4 model is practically unaffected by GR, while all the larger compactness models show non-negligible effects of GR.This effect is a function of shock radius, with a smaller shock radius leading to a larger difference in the NR and GR growth rates.This effect is drastically enhanced when looking at the ξ = 2.8 models.In these models, the power in the ℓ = 1 mode is some five orders of magnitude lower in the GR case by the end of the simulation (see Figure 8), which is the result of a factor of two difference in the growth rate along with exponential growth for approximately ten SASI cycles.The effects of GR can also be seen in Figure 9, which plots the growth rate for all models as a function of R sh .Again, in all cases, the growth rate is slower for the GR models and the difference between GR and NR growth rates increases with decreasing shock radius and increasing PNS mass -i.e., under more relativistic conditions -the ratio becoming as small as 0.47.
In Figure 10, we plot H 1 as a function of t for all models (left panels) and as a function of t/T for all models (right panels).Note that for an exponentially growing ℓ = 1 mode, H 1 ∝ exp(ωt); the gain in one cycle is |Q| = exp(ωT ) (e.g., see Janka 2017).Here, we simply refer to ωT = ln |Q| as the SASI efficiency.(When interpreting the right panels, note that T is different for each model.)These figures demonstrate that the SASI cycle efficiency is approximately constant for a given physical model (NR or GR) and compactness.This can be seen by, for example, first looking at the top-left panel of the left figure in Figure 10 and noting that the curves all grow at different rates and then comparing this to the top-left panel of the right figure, where the curves are brought much closer together when plotted against t/T .Further, all the NR models and the ξ ∈ {0.4,0.7} GR models reach about the same total power of ∼10 25 -10 26 after ten oscillations, while the ξ ∈ {1.8, 2.8} GR models reach a maximum of ∼10 23 , demonstrating that GR modifies both the oscillation period and the growth rate, but modifies them in such a way as to keep the SASI efficiency roughly constant for each of these groupings of models.
The SASI efficiency is further illustrated in Figure 11, which plots the efficiency for all of our models.It can be seen that the SASI efficiency of the NR models takes on values between 1.5 and 1.8, regardless of compactness.The GR models follow the same trend for lower com-pactnesses, but as the compactness increases the SASI efficiency drops.

SUMMARY, CONCLUSION, AND FUTURE WORK
We examined the effect of GR on 7 idealized, axisymmetric models of the standing accretion shock instability (SASI) by performing a parameter study in which we systematically varied the initial shock radius, along with the mass and radius of the proto-neutron star (PNS); i.e., the compactness of the PNS.We compared these runs by measuring the growth rate and oscillation period of the SASI for all models, which were run once with NRHD and Newtonian gravity and once with GRHD and GR gravity.We set up our simulations to excite a clean ℓ = 1 Legendre mode, computed the power of the SASI in that mode within a thin shell just below the shock, and, following Blondin & Mezzacappa (2006) the resulting signal to Eq. ( 41).From the fits, we computed the growth rates and their uncertainties, and using the FFT, we computed the oscillation periods and their uncertainties.
For all models, we found that the period of the SASI in the GR case is longer than that of the SASI in the NR case.For the ξ = 0.4, 0.7, 1.8, and 2.8 models, we find ratios, T gr /T nr , of 1.19, ≤ 1.15, 1.22, and ≤ 1.29, respectively.We explained these differences as resulting from differences in the post-shock flow structure for the GR and NR setups.For all models, we found that the growth rate is slower when GR is used, and significantly slower for the ξ = 1.8 and ξ = 2.8 models, with the ratio between the GR and NR cases, ω gr /ω nr , as low as 0.47.We found that both the growth rate and the oscillation period are practically independent of the accretion rate for the range of parameter space we considered.
The connection between our results and the results from realistic CCSN simulations can be made by considering the trends across all of the models considered here.First, our results suggest that CCSN simulations based on Newtonian hydrodynamics and Newtonian gravity may fail to predict correctly the growth rate of the SASI and its period of oscillation as conditions below the shock become increasingly relativistic; i.e., increased PNS compactness and decreased shock radius.Under such conditions, the growth rate in the Newtonian case may be overestimated by a factor of two, or more, relative to the GR case.Additionally, the period may be underestimated by about 20%.Second, as the conditions we considered became increasingly relativistic, the SASI growth rate continued to increase (Figure 9) and its period continued to decrease (Figure 7).Thus, our studies provide theoretical support for the conclusions reached in the study of the SASI development in higher-mass progenitors (see, e.g., Hanke et al. 2013;Matsumoto et al. 2022), where the time scales for the development of convection may be long and where the SASI, able to develop on shorter time scales, is able to provide support to the stalled shock, potentially sustaining neutrino heating, and, in turn, the development of neutrino-driven convection, with all of their anticipated benefits for generating an explosion.
Given the substantial impact of GR we find in our simulations, we stress the importance of GR treatments  in any future studies aiming to understand the SASI in a CCSN context.
Our results also suggest that CCSN simulations based on Newtonian hydrodynamics and Newtonian gravity may fail to predict correctly the emission of gravitational waves by the SASI (specifically, its frequency and amplitude), a primary target of gravitational wave astronomy given its anticipated existence in that part of frequency space where current-generation gravitational wave detectors such as LIGO and VIRGO are most sensitive.In addition, efforts to discern between contributions from different SASI modes -e.g., sloshing and spiral -may be affected, as well.
In light of the results presented here, further analysis is well motivated.Our assumption of axisymmetry will need to be lifted, since the SASI is known to have non-axisymmetric modes (Blondin & Mezzacappa 2007;Blondin & Shaw 2007;Fernández 2010Fernández , 2015)).A full three-dimensional comparison of the SASI in NR and GR, such as the one performed here for axisymmetry, should be conducted.Further analyses may also include adding a third type of model that uses NRHD and a GR monopole potential, similar to what is done in several CCSN simulation codes (e.g., Rampp & Janka 2002;Kotake et al. 2018;Skinner et al. 2019;Bruenn et al. 2020), in order to discern whether the use of an effective potential to capture the stronger gravitational fields present in GR is able to better capture the SASI growth rate and subsequent evolution.Of course, even if it does, nothing can replace a true GR implementation, as we have done here, and as must be done in future CCSN models.

Figure 1 .
Figure1.Post-shock, steady-state solutions for GR (solid) and NR (dashed) equations as functions of coordinate distance r ∈ [3 km, R sh ] for three models having the same accretion rate of 0.3 M⊙/s, the same PNS mass of 1.4 M⊙, and shock radii 120 km (blue), 150 km (orange), and 175 km (green).For this PNS mass, the Schwarzchild radius is ∼1 km.Quantities defined just ahead of the shock are denoted with a subscript "1".The top-left panel shows the comoving baryon mass density normalized to its value just ahead of the shock, the top-right panel shows the radial component of the fluid three-velocity normalized to its value just ahead of the shock, the bottom-left panel shows the comoving pressure normalized to the Newtonian ram pressure just ahead of the shock, ρ1 v 2 1 , and the bottom right panel shows the lapse function (solid) and the Newtonian approximation to the lapse function (dashed), 1 + Φ, with Φ the Newtonian gravitational potential, normalized to their values at the shock.

Figure 2 .
Figure2.Post-shock, steady-state solutions for GR (solid) and NR (dashed) equations as functions of coordinate distance r ∈ [3 km, R sh ] for three models having the same accretion rate of 0.3 M⊙/s, one having a PNS mass of 1.8 M⊙ and a shock radius of 70 km (blue), and the other two having PNS mass 2.8 M⊙ and shock radii 60 km (orange) and 70 km (green).For the 2.8 M⊙ PNS models, the Schwarzchild radius is ∼2.1 km.Quantities defined just ahead of the shock are denoted with a subscript "1".The top-left panel shows the comoving baryon mass density normalized to its value just ahead of the shock, the top-right panel shows the radial component of the fluid three-velocity normalized to its value just ahead of the shock, the bottom-left panel shows the comoving pressure normalized to the Newtonian ram pressure just ahead of the shock, ρ1 v 2 1 , and the bottom right panel shows the lapse function (solid) and the Newtonian approximation to the lapse function (dashed), 1 + Φ, with Φ the Newtonian gravitational potential, normalized to their values at the shock.

Figure 3 .
Figure3.Ratio of radial GR acoustic signal speed to radial NR acoustic signal speed (solid), ratio of angular GR acoustic signal speed to angular NR acoustic signal speed (dashed), and ratio of radial GR advective signal speed to radial NR advective signal speed (dotted), plotted as functions of η for all models.Lower-compactness models are shown in the top panel and higher-compactness models are shown in the bottom panel.All signal speeds have been multiplied with the appropriate scale factors.

Figure 6 .
Figure 6.Amplitudes of Fourier transforms of Eq. (42), normalized by their largest values, as functions of T in ms, where T is the inverse of the frequency returned by the FFT.The compactness of each model is shown in the bottom of the respective panels.NR results are shown in blue and GR results are shown in orange.

Figure 7 .
Figure 7. SASI oscillation period in ms as a function of PNS compactness (dimensionless) and shock radius in km, as determined by the FFT (see text for details).NR results are shown in open squares and GR results are shown in solid squares, with different colors corresponding to different compactnesses.

Figure 8 .
Figure8.H1 in cgs units versus t in ms.The layout of the panels is the same as in Figure6.NR results are shown in blue and GR results are shown in orange.Note that the horizontal axis limits are different in each panel.

Figure 9 .
Figure 9. SASI growth rate in ms −1 as a function of shock radius in km, as determined by the fit to Eq. (41) (see text for details).NR results are shown in open squares and GR results are shown in solid squares.

Figure 10 .Figure 11 .
Figure 10.H1 in cgs units versus t/ms (left) and versus t/T (right).In each four-panel figure, the top-left panel shows the ξ ∈ {0.4,0.7}, NR models; the top-right panel shows the ξ ∈ {0.4,0.7}, GR models; the bottom-left panel shows the ξ ∈ {1.8, 2.8}, NR models; and the bottom-right panel shows the ξ ∈ {1.8, 2.8}, GR models.A.M., and J.B. acknowledge support from the National Science Foundation's Gravitational Physics Program under grant 2110177.This research was supported in part by the Exascale Computing Project (17-SC-20-SC), a collaborative effort of the U.S. Department of

Table 1 .
Model ParametersNote-Model parameters chosen for the 7 models.All models were run with both GR and NR.