# Effects of rotation and magnetic field on the revival of a stalled shock in supernova explosions

We investigate axisymmetric steady solutions of (magneto)hydrodynamics equations that describe approximately accretion flows through a standing shock wave and discuss the effects of rotation and magnetic field on the revival of the stalled shock wave in supernova explosions. We develop a new powerful numerical method to calculate the 2-dimensional (2D) steady accretion flows self-consistently. We first confirm the results of preceding papers that there is a critical luminosity of irradiating neutrinos, above which there exists no steady solution in spherical models. If a collapsing star has rotation and/or magnetic field, the accretion flows are no longer spherical owing to the centrifugal force and/or Lorentz force and the critical luminosity is modified.In fact we find that the critical luminosity is reduced by about 50 70 fields, depending on the mass accretion rate. These results may be also interpreted as an existence of the critical specific angular momentum or critical magnetic field, above which there exists no steady solution and the standing shock wave will revive for a given combination of mass accretion rate and neutrino luminosity.

## Authors

• 4 publications
• 4 publications
• 4 publications
• 5 publications
04/17/2020

### A Multistate Low-dissipation Advection Upstream Splitting Method for Ideal Magnetohydrodynamics

We develop a new numerical scheme for ideal magnetohydrodynamic (MHD) si...
03/29/2021

### Higher dimensional generalization of the Benjamin-Ono equation: 2D case

We consider a higher-dimensional version of the Benjamin-Ono (HBO) equat...
01/06/2022

### A Divergence-Conforming Hybridized Discontinuous Galerkin Method for the Incompressible Magnetohydrodynamics Equations

We introduce a new hybridized discontinuous Galerkin method for the inco...
04/05/2022

### A gyroscopic polynomial basis in the sphere

Standard spectral codes for full sphere dynamics utilize a combination o...
04/28/2020

### A new 3D model for magnetic particle imaging using realistic magnetic field topologies for algebraic reconstruction

We derive a new 3D model for magnetic particle imaging (MPI) that is abl...
02/02/2021

### Global Earth Magnetic Field Modeling and Forecasting with Spherical Harmonics Decomposition

Modeling and forecasting the solar wind-driven global magnetic field per...
02/24/2021

### Survival of the Fittest: Numerical Modeling of Supernova 2014C

Initially classified as a supernova (SN) type Ib, ∼ 100 days after the e...
##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Introduction

Core-collapse supernovae (CCSNe) play important roles in various fields of astrophysics such as the star formation, galactic evolution and the acceleration of cosmic-ray through their nucleosynthesis, energetic shock waves and high luminosities. The large gravitational energy liberated in the collapse makes CCSNe promising sites for the emissions of neutrinos and gravitational waves as well as heavy elements. The understanding of the physical processes in and the mechanism of CCSNe is important for the development of multi-messenger astronomy.

CCSNe commence with the gravitational collapse of the core of massive stars at the ends of their lives. When the central density reaches the nuclear saturation density, core bounce occurs and a shock wave is generated. As it propagates outward, the shock loses energy via photo-dissociations of nuclei and neutrino cooling and is eventually stagnated inside the core by the ram pressure of accreting matter in addition to these energy losses. The stalled shock wave should be revived somehow to produce a successful explosion. It is widely expected that it will be re-energized by the irradiation of matter behind it by neutrinos diffusing out of a proto-neutron star (PNS) (Wilson, 1985). This neutrino heating scenario is currently the most favored mechanism of CCSNe (Janka et al., 2012; Kotake et al., 2012; Burrows, 2013; Foglizzo et al., 2015; Müller, 2016; Burrows et al., 2018 for recent reviews).

Intensive and extensive efforts in numerical simulations have revealed so far that no successful explosion obtains in spherical symmetry (e.g. Liebendörfer et al., 2001, 2005; Rampp & Janka, 2002; Sumiyoshi et al., 2005) and multi-dimensional effects are crucially important (Burrows et al. 2006; Marek & Janka 2009; Bruenn et al. 2009; Suwa et al. 2010; Müller et al. 2012; Takiwaki et al. 2012; Hanke et al. 2013; Murphy et al. 2013; Couch 2013; Couch & Ott 2013; Lentz et al. 2015; Nakamura et al. 2015; Melson et al. 2015; Roberts et al. 2016; Bruenn et al. 2016; O’Connor & Couch 2018). Important effects among them are rotation (Fryer & Heger 2000; Kotake et al. 2003; Thompson et al. 2005; Marek & Janka 2009; Iwakami et al. 2014a; Nakamura et al. 2014; Takiwaki et al. 2016; Summa et al. 2018), magnetic field (Akiyama et al. 2003; Yamada & Sawai 2004; Kotake et al. 2004; Sawai et al. 2005; Obergaulinger et al. 2006; Burrows et al. 2007; Takiwaki et al. 2009; Sawai & Yamada 2014; Obergaulinger et al. 2014; Mösta et al. 2015; Sawai & Yamada 2016; Obergaulinger et al. 2018), non-spherical structures of the progenitor (Couch & Ott 2013; Takahashi & Yamada 2014; Couch et al. 2015; Takahashi et al. 2016), turbulence (Murphy & Burrows 2008; Murphy & Meakin 2011; Murphy et al. 2013; Couch & Ott 2015; Mabanta & Murphy 2018), (magneto)hyrdodynamical instabilities (Blondin et al. 2003; Scheck et al. 2006; Blondin & Mezzacappa 2007; Iwakami et al. 2008; Guilet et al. 2010; Wongwathanarat et al. 2010; Takiwaki et al. 2014; Fernández et al. 2014; Fernández 2015), general relativistic gravity (Dimmelmeier et al. 2002; Shibata & Sekiguchi 2004, 2005, Müller et al. 2012; Ott et al. 2012; Kuroda et al. 2012, 2016) and neutrino transport (Nagakura et al. 2014; Dolence et al. 2015; Pan et al. 2016; Nagakura et al. 2017, 2018). It is true that large-scale dynamical simulations have played a crucial role in the recent progresses in our understanding of these ingredients, but we believe that a more phenomenological approach that employs toy models still plays an indispensable and complementary role to understand each effect more deeply.

Burrows & Goshy (1993) took such an approach and introduce the concept of the critical neutrino luminosity. They approximated the accretion flows after the shock wave is stagnated with spherical steady solutions with constant mass accretion rates of the hydrodynamics (HD) equations that incorporate neutrino heating of matter; they found then that for a given mass accretion rate, there is a critical neutrino luminosity, above which no steady solution exists; they surmised that the revival of the stalled shock wave should occur there. The existence of such a critical neutrino luminosity was confirmed both in similar phenomenological methods (Yamasaki & Yamada, 2005, 2007; Keshet & Balberg 2012; Murphy & Dolence 2017) and in numerical simulations (Janka & Mueller, 1996; Ohnishi et al. 2006; Iwakami et al. 2008; Murphy & Burrows 2008; Nordhaus et al. 2010; Iwakami et al. 2014b, a; Suwa et al. 2016).

For examples in the phenomenological models, Pejcha & Thompson (2012) contended that the ante-sonic condition for the ratio of the sound velocity to the escape velocity is an important diagnostic for explosion. (see also Raives et al. 2018). Nagakura et al. (2013) argued, on the other hand, that there is a critical fluctuation for a given accretion flow, above which the shock wave can sustain its outward propagation once launched. Gabay et al. (2015) applied a phase-space analysis to the shock acceleration in the shock radius and shock velocity plane. More recently Murphy & Dolence (2017) contended that five parameters (mass accretion rate, neutrino luminosity, neutrino temperature, PNS mass and radius) instead of two (mass accretion rate and neutrino luminosity alone) are more appropriate in considering the condition for explosion. They also derived a diagnostic dimensionless quantity from the integrability condition, which roughly measures the balance between pressure and gravity. It should be stressed here that all these conditions are based on the analysis of spherical accretion flows and should be modified if the background flow is non-spherical.

In fact, it was shown that the critical luminosity is lowered by rotation both in the phenomenological method (Yamasaki & Yamada, 2005) and dynamical simulations (Murphy & Burrows, 2008; Nordhaus et al. 2010; Hanke et al. 2013; Couch 2013; Iwakami et al., 2014b, a). Yamasaki & Yamada (2005) extended their spherically-symmetric, steady accretion-flow to axisymmetric ones, taking rotation into account consistently. They showed that the accretion flows become non-radial in the presence of rotation and the critical luminosity is lowered. Iwakami et al. (2014a), on the other hand, performed 3-dimensional (3D) dynamical simulations of rotational accretions and obtained similar results: the critical luminosity is lowered by rapid rotation and there is a critical value of specific angular momentum, above which the shock wave is revived, for a given combination of mass accretion rate and neutrino luminosity. This implies that an explosion may be triggered by rapid rotation even if the neutrino luminosity is sub-critical without rotation.

It is emphasized here that Yamasaki & Yamada (2005) is the only one that has so far taken the phenomenological approach to the non-spherical accretion flows. This is in sharp contrast to the fact that there have been many numerical simulations that are meant to address the issue (Janka & Mueller, 1996; Ohnishi et al. 2006; Iwakami et al. 2008; Murphy & Burrows 2008; Nordhaus et al. 2010; Hanke et al. 2013; Iwakami et al. 2014b, a; Suwa et al. 2016

). The reason for this is probably the difficulty to obtain numerically non-spherical steady solutions of the HD equations. Even in

Yamasaki & Yamada (2005), the number of solutions was not large. Not to mentions, there has been no attempt to incorporate magnetic field in such approaches.

In this paper, we develop a new numerical scheme to numerically obtain axisymmetric steady solutions of HD and magnetohydrodynamics (MHD) equations that describe self-consistently non-spherical post-shock accretion flows in the presence of rotation and/or magnetic field. Based on these solutions we study the effects of rotation and magnetic field on the critical neutrino luminosity systematically. This paper is organized as follows. In Sec. 2 we formulate the problem and give the basic equations to be solved. We then explain our new numerical method. In Sec. 3 we present numerical solutions and the critical luminosities obtained from them. Finally we give some discussions and summarize this paper in Sec. 4.

## 2 Formulations

### 2.1 Assumptions and basic equations

After the shock stalls, the mass accretion rate and neutrino luminosity change slowly and the post-shock accretion flows may be approximated as steady states with constant mass accretion rates and neutrino luminosities (Burrows & Goshy 1993; Yamasaki & Yamada 2005). We make the following assumptions to construct idealized models in this paper.

1. The equatorial symmetry is assumed in addition to the stationary () and axisymmetry ().

2. We ignore both viscosity and magnetic resistivity and consider only ideal HD or MHD equations.

3. We consider the accretion flows only in the domain between the stalled shock wave and the neutrino sphere and impose the outer and inner boundary conditions there.

4. The PNS is treated as a point mass with . Only its Newtonian gravitational forces are considered and the self-gravity of accreting matter is ignored for simplicity.

5. We neglect the neutrino transfer again for simplicity and the luminosity and energy spectrum of neutrinos are assumed to be constant in radius (see Equations 4 and 5).

6. A simplified equation of state (EOS) that neglects the photo-dissociations of nuclei is adopted and the convection is not taken into account.

Under these assumptions, the basic equations are given as

 1r2 ∂∂r(r2ρur)+1rsinθ∂∂θ(sinθρuθ)=0, (1a) ur ∂ur∂r+uθr∂ur∂θ−u2θ+u2φr=−1ρ∂p∂r−GMr2 −14πρ[Bθ∂Bθ∂r+Bφ∂Bφ∂r+B2θ+B2φr], (1b) ur ∂uθ∂r+uθr∂uθ∂θ+uruθr−u2φcotθr=−1ρr∂p∂θ +14πρ[Br∂Bθ∂r−Brr∂Br∂θ −Bφr∂Bφ∂θ+BrBθ−B2φcotθr], (1c) ur ∂uφ∂r+uθr∂uφ∂θ+uφurr+uθuφcotθr =14πρ[Br∂Bφ∂r+Bθr∂Bφ∂θ+BrBφ−BθBφcotθr], (1d) ur (∂ε∂r−pρ2∂ρ∂r)+uθr(∂ε∂θ−pρ2∂ρ∂θ)=˙q, (1e)

where and denote, respectively, the density, velocity, magnetic field, pressure, specific internal energy and neutrino-heating rate; and are the gravitational constant and the mass of PNS. In the absence of the magnetic field, is simply set to zero. In its presence, on the other hand, we also solve the Maxwell equations under the conditions of stationary and ideal MHD (i.e. ) as follows:

 1r2 ∂∂r(r2Br)+1rsinθ∂∂θ(sinθBθ)=0, (2a) − uθ∂Br∂θ+ur∂Bθ∂θ−Br∂uθ∂θ +Bθ∂ur∂θ+urBθcotθ−uθBrcotθ=0, (2b) r uθ∂Br∂r−rur∂Bθ∂r−rBθ∂ur∂r +rBr∂uθ∂r+Bruθ−Bθur=0, (2c) r Br∂uφ∂r−rBφ∂ur∂r+ruφ∂Br∂r−rur∂Bφ∂r +Bθ∂uφ∂θ−Bφ∂uθ∂θ+uφ∂Bθ∂θ −uθ∂Bφ∂θ+Bruφ−Bφur=0, (2d)

where is the electric field. Following Yamasaki & Yamada (2005), we employ a simplified EOS given as follows:

 p=11π2180k4c3ℏ3T4+ρkTmN, (3a) ε=11π260k4c3ℏ3T4ρ+32kTmN, (3b)

where , and are the Boltzmann constant, speed of light, reduced Planck constant, nucleon mass and matter temperature, respectively. Note that this EOS neglects photo-dissociations of nuclei and, as a consequence, the critical luminosity tends to be overestimated. Since we are interested in the qualitative effects of rotation and/or magnetic field on the critical luminosity in this paper, the simplification is justified. The net neutrino-heating rate is also simplified as in Herant et al. (1992)

 ˙q = 4.8×1032⎡⎣1−√1−r2νr2⎤⎦Lν2πr2νT2ν (4) − 2.0×1018T6 [ergs s−1 g−1],

where the neutrino luminosity and temperature are model parameters, which are constant both in time and space, and common to and ; we assume further that they are related with the radius of neutrino sphere , which is coincident with the inner boundary in our models, as

 Lν=74πr2νσST4ν, (5)

where denotes the Stefan-Boltzmann constant. The neutrino temperature, which is supposed to characterize the neutrino energy spectrum, is set to in the following (Yamasaki & Yamada 2005).

### 2.2 Boundary conditions

The region of our interest is the one between the neutrino sphere and the stalled shock wave, and the inner and outer boundary conditions are imposed there, respectively. Following Yamasaki & Yamada (2005), we impose on the inner boundary, which approximately corresponds to the condition that the optical depth to the neutrino sphere from infinity is equal to . The latter condition was adopted by Burrows & Goshy (1993) and Murphy & Burrows (2008) (see Keshet & Balberg 2012 for a comparison of these two conditions). The radius of the neutrino sphere is obtained from the neutrino luminosity and temperature in Equation (5).

The main difficulty in obtaining steady accretion flows through the stalled shock wave in multi-dimensions stems from the simple fact that the shock surface is non-spherical. In this paper, we assume that the shape of shock surface expressed as in the polar coordinates is given with the Legendre function as

 rs(θ)= rsph{1−ϵP2(cosθ)}, (6a) P2(cosθ)= 12(3cos2θ−1), (6b)

where and

are the equatorial radius and deformation amplitude of the shock surface, respectively. Instead of attempting to obtain the shock surface directly, we employ these two degrees of freedom alone, which turn out to be sufficient. (see Sec.

2.4). Once the shock shape is determined, the post-shock quantities are obtained from the pre-shock quantities via the Rankine-Hugoniot relations, which are given in Appendix A. We assume that the pre-shock flows are spherically symmetric free falls from infinity with different mass accretion rates. We do not attempt in this paper to elaborate on the outer flow so that rotation and/or magnetic field would be taken into account consistently.

Other boundary conditions are administered at and . As in Yamasaki & Yamada (2005), we impose Neumann conditions for at , and for at while we set Dirichlet conditions for , , , at and for and also at as follows:

 uθ=Bθ=uφ=Bφ=0 (7a) ∂ρ∂θ=∂ur∂θ=∂T∂θ=∂Br∂θ=0 (7b)

at and

 uθ=Bθ=0 (8a) ∂ρ∂θ=∂ur∂θ=∂T∂θ=∂Br∂θ=∂uφ∂θ=∂Bφ∂θ=0 (8b)

at respectively.

### 2.3 Rotation law and magnetic field profiles

In the presence of magnetic field, the rotation law and the magnetic field profile should be given consistently with each other, since there are some conserved quantities along the streamline, which constrain the variations of rotation and magnetic field, for the axisymmetric and steady ideal MHD flows (see Lovelace et al. 1986; Fujisawa et al. 2013). The details of their derivations are described in Appendix B.

If a magnetic field is purely toroidal (), the specific angular momentum and a quantity related with magnetic torque are conserved along the streamline:

 ℓ(ψ)=rsinθuφ, (9a) σ(ψ)=Bφrsinθρ, (9b)

where is a stream function defined as

 ur≡14πρr2sinθ∂ψ∂θ, uθ≡−14πρrsinθ∂ψ∂r. (10)

Instead of fixing the functional forms of and explicitly, we give the rotational velocity and magnetic field strength at the outer boundary to satisfy Equations (7) and (8). Following Yamasaki & Yamada (2005), we assume that a spherical shell is rotating rigidly at a radius of as

 uφ(r1000km,θ)=2πr1000kmsinθf1000km. (11)

Since the specific angular momentum is conserved along an individual streamline, the rotational velocity just above the shock surface is given as

 uφ(r,θ)=2πr21000kmrsinθf1000km, (12)

if the streamlines are assumed to be radial. Yamasaki & Yamada (2005) obtained steady solutions only for (slow rotation) and (fast rotation). We consider a larger variety than Yamasaki & Yamada (2005) in this paper.

As for the toroidal magnetic field, we assume the following profile at the same radius of as

 Bφ(r1000km,θ) =B1000kmsinθ, (13)

where givens the strength of the toroidal magnetic field there. Since the in Equation (9) is a conserved quantity, the toroidal magnetic field just outside the shock surface is determined as

 Bφ(r,θ) =B1000km√r1000kmrsinθ, (14)

if we again assume a radial infall upstream the shock wave (see Appendix B). Note that these functional forms in Equations (11) and (13) satisfy the boundary conditions in Equations (7) and (8).

If the magnetic field has both poloidal and toroidal components, both of them should be fixed at the outer boundary. If we again assume a radial free fall outside the shock, the poloidal magnetic field lines should be also radial. Then the continuity equation of the magnetic field gives the poloidal magnetic field just outside the shock surface as

 Br(r,θ)=r21000kmr2B0,Bθ(r,θ)=0, (15)

where is a constant giving the strength of the poloidal magnetic field at the radius of . When the poloidal magnetic field is non-zero, the specific angular momentum and the toroidal magnetic field are no longer independent of each other (Fujisawa et al. 2013) and

Equations (12) and (14) are not valid. The values of and just ahead of the shock surface are then obtained by solving Equations (B13a) and (B15a) numerically in Appendix B.

### 2.4 Numerical scheme and parameter settings

We summarize our new numerical scheme to obtain steady solutions briefly here. More details are given in Appendix C. Each steady solution is characterized by five parameters, , , , and . We discretize Equations (1) and (2) and obtain non-linear systems of equations. Setting the above five parameters and assuming the values of and in Equation (6), we solve these equations numerically from the shock surface to down the neutrino sphere. If the density at the neutrino sphere so obtained does not satisfy the requirement , we modify the values of and and repeat the procedure. and

are hence regarded as the eigenvalues in the boundary-value problem for this system. In order to obtain the value of the critical neutrino luminosity, we calculate a sequence of solutions for fixed values of four parameters (

, , and ), changing the value of until a steady state solution is no longer obtained. We develop a new numerical method dubbed the W4 method (Okawa et al. 2018), which is meant to solve the non-linear systems of equations. The details of the W4 method are described in Appendix D (see also our recent paper Okawa et al. 2018).

Information on the rotation and magnetic field deep inside the progenitor is scarce. In general, massive stars, such as B and O type stars, tend to have rapid rotation (e.g. Hunter et al. 2008). Approximately 10 per cent of them have surface rotational velocities larger than 300 (Ramírez-Agudelo et al. 2013). On the other hand, recent surveys for massive stars indicate that some Galactic O and B type stars have magnetic fields of 100 – 1000 at their surfaces (Wade & MiMeS Collaboration, 2015). Their total magnetic fluxes roughly coincidence with those of a typical magnetar whose magnetic field is about on their surfaces. According to the stellar evolution models, the toroidal magnetic field is likely to be larger by orders of magnitude than the poloidal one due to the differential winding inside the massive stars (Heger et al. 2005). Some progenitors may hence have a rapidly rotating and/or strongly magnetized core.

On the other hand, the stellar evolution models also indicate that the transport of angular momentum during the quasi-static phase of the progenitor reduces the angular momentum of the core, particularly if the magnetic torque is taken into account (e.g. Heger et al., 2005; Potter et al., 2012). The rotation might not be so rapid after all then. If this is really the case, rotation does not play an important role in the dynamics of core collapse. One should bear in mind, however, that almost all stellar evolution calculations are based on spherically averaged one-dimensional models and have uncertainties in the mechanism and formulation of angular momentum transport and magnetic field.

The aim of this paper is to systematically study the effects of rotation and magnetic field on the critical neutrino luminosity in CCSNe using the new numerical scheme. We consider rapid rotation and/or strong toroidal magnetic field, setting the parameters as , and . These values roughly correspond to and on the PNS surface, which are the typical rotation frequency of a milisecond pulsar and the canonical strength of the surface magnetic field of a magnetar respectively. The PNS neutrino temperature and the mass are fixed to and , respectively in this paper.

## 3 Numerical results

### 3.1 Streamlines of steady accretion flow

Figure 1 displays streamlines in a meridian plane for a model either with rotation (left panel) or with toroidal magnetic field (right panel). In these models the neutrino luminosity and accretion rate are set as and . The innermost black-dotted curve indicates the neutrino sphere whereas the outermost black-dashed one corresponds to the stalled shock surface. We also draw the contour lines for the azimuthal components of velocity and magnetic field in the left and right panels, respectively. The values of and (see Equation 6) are  km and respectively, in both models. Interestingly, although the shapes of the shock surfaces are almost the same, being oblate, the flow pattern are different from each other. In fact in the left panel of Figure 1, the flow is pushed toward the equatorial plane as it approaches the PNS. This is due to centrifugal force and was previously seen in Figure 6 in Yamasaki & Yamada (2005). In the right panel, on the other hand, the flow is bent toward the symmetry axis. This is due to Lorentz force exerted by the toroidal magnetic field, which indeed behaves as a negative centrifugal force (e.g. Fujisawa 2015). Noted that the flow patterns of both models are similar near the rotation axis because of the boundary condition for at the axis (Equation 7).

The left panel in Figure 2 shows the result for a model with both rotation and toroidal magnetic field. The mass accretion rate, neutrino luminosity and strength of toroidal magnetic field are set to the same as values in Figure 1, but the rotation frequency is somewhat smaller. We find , slightly larger than in the previous cases in Figure 1. The centrifugal force and the Lorentz force almost cancel each other in this solution. The streamlines are nearly radial in this model except for the inner region where the Lorentz force is dominant over the centrifugal force and the streamlines look similar to those for purely toroidal magnetic field. This is understood from the conserved quantities and along the streamline. As a matter of fact, does not depend on the density while does as is clear in Equation (9). Since the density is low near the shock surface and gets higher toward the neutrino sphere, where . The Lorentz force tends to be dominant in the inner region and vice versa.

The right panel in Figure 2 shows streamlines for a model that incorporates poloidal magnetic fields in addition to toroidal ones. The neutrino luminosity and accretion rate are set as and , respectively. The rotational frequency and the strength of magnetic fields are , and , respectively. Since the value of is smaller by an order than that of , the poloidal magnetic field in this model is much weaker than the toroidal magnetic field. We find that the value of is again. Note that the poloidal magnetic field lines are aligned with the streamlines in ideal MHD. They are slightly curved near the neutrino sphere similarly to the left panel in Figure 1 because of Lorentz force mainly exerted by the toroidal magnetic field.

Figure 3 displays the radial and the azimuthal components of velocity (left) and magnetic field (right) at and . The radial components of the flow velocities and of the poloidal magnetic field are almost identical at these zenith angles. This is mainly because of our assumption that the flow and the poloidal magnetic field outside the shock surface are radial and independent of . If we had assumed a -dependent functional form in Equation (15), for example, we would have found -dependent radial profiles of . In contrast, both and depend on , being larger on equator () than at . This dependence again (), however, simply reflects of the functional forms for and in Equations (11) and (13).

The rotational velocity and toroidal magnetic field both increase inwards from the shock surface to the neutrino sphere. It is apparent, however, that the toroidal magnetic field rises more steeply than the rotational velocity. This trend is almost independent of the functional forms for and . In fact, it is dictated by the conservation of and , the latter of which depends on the density profile as is explicit in Equation (B11) in Appendix B. One may hence roughly say that the dependence of the flow velocities and magnetic fields is mainly determined by their dependence just ahead of the shock wave, while their radial profiles are largely constrained by the conservation laws (Equation B11), being almost independent of the outer boundary conditions.

### 3.2 Effects of rotation and magnetic field on the critical luminosity

Now we look into the changes in the critical luminosity that rotation and/or magnetic field produce. In this subsection we neglect poloidal magnetic fields () because the toroidal magnetic field is supposed to be dominant. Figure 4 displays the critical luminosities either for models with rotation alone (left panel) or for models with toroidal magnetic field alone with no rotation (right panel), which are plotted as a function of the mass accretion rate. As a reference, we also plot the critical luminosity of the spherical symmetric models without rotation and toroidal magnetic field as black circles.

It is evident among other things that the critical luminosity for all these models with either rotation or toroidal magnetic field are lower than those for the spherical models. Non-spherical forces, i.e., centrifugal forces and hoop stresses, tend to reduce the critical luminosity although the reduction rate depends on the mass accretion rate. As a matter of fact, for a given rotation velocity or magnetic field strength at the outer boundary, the critical luminosity gets lowered more strongly as the mass accretion rate becomes smaller: the reduction rate is as high as about for the rapid rotation with and for the strong toroidal magnetic field with at the mass accretion rates lower than .

Interestingly, the situation is other way around in the right panel for the purely toroidal magnetic field. This is because such fields exert hoop stress, which behaves like a negative centrifugal force. As a result, the radial inflow velocity is smaller on the equator than on the axis near the neutrino sphere. According to the same argument, this may imply that only the equatorial flow becomes critical state locally and shock revival will be initiated on the equator. It is intriguing that the existence/non-existence of steady solutions that satisfy a certain boundary condition is itself a global issue while the critical state is realized locally as a result. Although some recent studies strongly deduce the presence of global conditions under critical state in spherical symmetry (e.g. Murphy & Dolence 2017), one may hence claim that the critical neutrino luminosity is determined both globally and locally.

Figure 6 displays the streamlines (left panel) and the radial profiles of the net heating rates along some radial rays (right panel) in the meridian section for one of the purely rotational models at its critical luminosity. The thick orange curve in the left panel indicates the boundary between the cooling and heating regions, or the 2D analog of the gain radius. One finds in the right panel that the cooling region near the pole (see the green line) is wider and deeper than that on the equator although the difference is not very large. Noted that for simplicity we assume in this paper that the neutrino sphere is spherical (see also Equation 4), which is certainly not true and tends to reduce the difference. We emphasize, however, that this slightly deformed cooling region results in the reduction of the critical neutrino luminosity as observed earlier.

Figure 7 shows distribution of the plasma defined as

 β=8πp|Bφ|2. (16)

in the meridian section for one of the non-rotating but purely toroidally-magnetized models close to its critical luminosity. The value of is very high, , almost everywhere in the computational domain, meaning that the magnetic pressure is much smaller than the matter pressure except near the equatorial area on the neutrino sphere where the value of is much smaller. This locally strong magnetic pressure (or locally low ) modifies the accretion flow there and gives rise to the critical state locally, leading to the reduction of the critical luminosity, as shown in the right panel of Figure 5.

Finally we show in Figure 8 the critical luminosity as a function of the rotational rate (left) and the strength of the toroidal magnetic field (right). We fix the mass accretion rate as . The critical luminosities of solutions with or are less than . The critical luminosity is reduced by for the most rapid rotation and by for the strongest toroidal magnetic field, respectively. It is interesting that we cannot find a steady state solution for the rotation parameter or the toroidal magnetic field strength because the shock surface comes too close to the neutrino sphere. These numerical results may imply that the existence of the critical specific angular momentum and the critical strength of the toroidal magnetic field, above which there actually exist no steady solutions.

## 4 Discussion and summary

We have derived numerically steady, non-spherical accretion flows through the standing shock wave onto the PNS solutions of stalled shock wave in the CCSNe core and effects of rotation and magnetic field on the shock revival in this paper. In order to obtain these steady solutions, we have developed a new numerical scheme to solve a system of nonlinear equations numerically, which we named the W4 method for non-linear equations, and we have indeed succeeded in generating various accretion flows with both rotation and magnetic field incorporated self-consistently in axisymmetric 2D. It should be noted that our new method can handle both poloidal and toroidal magnetic fields simultaneously.

Our main findings are summarized as follows.

1. The shock surface and the flow pattern become non-spherical by rotation and/or magnetic field in general. The shock surface is always oblate () whereas the stream lines are bent either toward the equatorial plane by rotation and poloidal magnetic field or toward the pole plane by toroidal magnetic field (Figures 1 and 2).

2. The toroidal magnetic field is dominant in the inner region (Figure 7) while the rotation is dominant in the outer region (Figures 2 and 3) because of the conserved quantities and in Equations (9) and (B11).

3. The critical luminosity is lowered by rotation and/or magnetic field in general although the degree of reduction depends on the mass accretion rate (Figure 4).

4. In the presence of rotation and/or magnetic field the critical state is realized locally: either on the symmetry axis for rotation or on the equatorial plane for toroidal magnetic field (Figure 5) despite the solution itself is determined globally according to the boundary conditions.

5. The gain region in the accretion flow is also deformed. The cooling region is widened near the pole by the rotation (Figure 6) whereas it is affected mostly near the equator by the toroidal magnetic field (Figure 7). Although the deviation from the spherical symmetry is small, it results in the substantial reduction of the critical neutrino luminosity.

6. It is suggested the existence of the non-vanishing critical specific angular momentum, above which no steady solution exists irrespective of the neutrino luminosity (left panel in Figure 8). Our results are roughly consistent with the dynamical simulations in Iwakami et al. (2014a). We found that there also exists the critical strength of the toroidal magnetic field, above which no steady solution exists (right panel in Figure 8).

We ignore other non-spherical physical processes such as the turbulence in progenitors in this paper. Turbulence in the accretion flow enables explosion by changing the flow property. Murphy & Burrows (2008) performed 1D and 2D dynamical simulations and suggested that the reduction of critical luminosity is caused by turbulence. Murphy & Meakin (2011) examined many turbulent models using Reynolds decomposition and proposed a global turbulent model that reproduces the profiles and evolution of the simulations (see also 2D and 3D dynamical simulations by Murphy et al., 2013). Recently, Mabanta & Murphy (2018) considered 1D spherical steady models with the turbulence and contended that the turbulent dissipation rather than the ram pressure from the turbulence reduces the critical luminosity. The turbulent dissipation is a significant contributor to successful supernova explosion. Note that these works ignored the effect of magnetic fields. Masada et al. (2015), on the other hand, performed high resolution local simulations with magnetic field in a 3D thin layer to calculate the magneto-rotational-instability near the neutrino sphere accurately and found that the convectively stable layer around the neutrino sphere becomes fully turbulent due to the MRI. Magnetic field plays an important role on the turbulence.

The steady state problems concerning these physical processes, however, are considered in only 1D spherical models (e.g. Mabanta & Murphy 2018) and nobody has succeeded in obtaining the 2D steady accretion models with them. We will study the steady flows with them in 2D self-consistently using our numerical procedure and analyze the effect from them on the critical neutrino luminosity in the future.

## Acknowledgements

The authors thank Dr. Kazuya Takahashi and Dr. Wakana Iwakami for helpful discussions. This work was supported by JSPS KAKENHI Grant Numbers 16K17708, 16H03986, 17K18792. KF was supported by JSPS Postdoctoral Fellowship for Research Fellowship (16J10223).

## Appendix A Rankine-Hugoniot relation on the deformed shock wave

We need to impose the junction condition on the shock surface for each variable. We first define two unit vectors,

and , each normal and tangential, respectively, to the shock surface given as . Then is defined to be the angle between the radial unit vector and as shown in the left panel of Figure 9. The value of is obtained as follows:

 →u∥⋅→er= |→u∥|cosχ, (A1a) cosχ= drsdθ{r2s+(drsdθ)2}−12. (A1b)

It is also useful to adopt the so-called surface-fitted coordinates defined from the spherical coordinates () as

 q≡r−rνrs(θ)−rν, (A2)
 θ′=θ. (A3)

This transformation maps the region between the neutrino sphere and the deformed shock surface into a domain given simply as

 q∈[0,1], θ′∈[0,π]. (A4)

See the right panel of Figure 9. Noted that this is the domain of our main concern in this paper. The derivatives are also transformed as follows:

 ∂Q∂r= 1rs−rν∂Q∂q, (A5a) ∂Q∂θ= ∂Q∂θ′−qrs−rνdrsdθ∂Q∂q, (A5b)

where is a physical variable such as density or velocity. Henceforth, we will employ these surface-fitted coordinates alone and use the same notation instead of to denote the new angle coordinate.

The MHD Rankine-Hugoniot relation (e.g. Winterhalter et al. 1984; Takahashi & Yamada 2013) is written as follows:

 [ρu⋅n]=0, (A6a) (A6b) [u⋅n(ρε+p+12ρu2+B24π) −14π(B⋅n)(B⋅u)]=0, (A6c) [B⋅n]=0, (A6d) [n×(u×B)]=0, (A6e)

where the background denotes as usual the difference between the upstream and downstream values of the quantity given inside the bracket and we use the ordinary notation for the normal vector to the shock surface. One can also express these relations as

 ρu⊥=^ρ^u⊥, (A7a) ρu⊥u+(p+B28π)n−14πB⊥B =^ρ^u⊥^u+(^p+^B28π)n−14π^B⊥^B, (A7b) u⊥(ρϵ+p+12ρu2+B24π)−14πB⊥(B⋅u) =^u⊥(^ρ^ϵ+^p+12^ρ^u2+^B24π)−14π^B⊥(^B⋅^u), (A7c) B⊥=^B⊥, (A7d) B⊥u−u⊥B=^B⊥^u−^u⊥^B, (A7e)

where the hat means the upstream value and the subscript denotes the normal component to the shock surface. We use these relations to obtain the downstream values of physical quantities. The rotation and magnetic field upstream the shock are given in Section 2.3 (see also Appendix B).

Since we have assumed a free-falling cold, spherical flow outside the shock surface, the radial velocity , density and pressure just ahead of the shock surface are given as follows:

 uf(r)= −√2GMr, (A8a) ρf(r)= −˙M4πr2uf=˙M√32π2GMr−32, (A8b) pf(r)= 0. (A8c)

Note that the mass accretion rate is positive () in this paper.

## Appendix B The rotation law and magnetic-field profile outside the shock wave

Here we explain in detail the ideas behind the setting of the outer boundary conditions for rotation and magnetic field. If the magnetic field is purely toroidal (), there are two conserved quantities: specific angular momentum and a quantity related with the magnetic torque . The former comes from the component of the Euler equation (Equation 1) and the latter results from the component of the ideal MHD condition. They are constant along each streamline. It should be noted that the specific entropy is not conserved, since heating and cooling are taken into account in our formulation unlike the previous works that assumed the barotropic (adiabatic) condition (e.g. Lovelace et al. 1986; Fujisawa et al. 2013). The conservations of and are expressed as

 ℓ(ψ)=rsinθuφ, (B1a) σ(ψ)=Bφrsinθρ, (B1b)

where and are arbitrary functions of , which is a so-called stream function defined in the following relations:

 ur≡14πρr2sinθ∂ψ∂θ, uθ≡−14πρrsinθ∂ψ∂r. (B2)

Since we have assumed that the flow is spherical outside the shock surface, the stream function there is given explicitly as

 ψ=˙Mcosθ. (B3)

We need to specify the functional forms of and to fix the rotation law and the profile of toroidal magnetic field. In this paper following Yamasaki & Yamada (2005) we specify the angular distributions of and at . Since the stream function depends only on alone (see Equation B3), and are also the functions of alone and the and at are expressed as follows:

 uφ(r1000km,θ) =ℓ(θ)r1000kmsinθ, (B4a) Bφ(r1000km,θ) =r1000kmsinθρ1000kmσ(θ) =˙M√32π2GMr1000kmsinθσ(θ), (B4b)

where is the density at . Just as in Yamasaki & Yamada (2005) we assume uniform rotation at with a rotational frequency . The is given as

 ℓ(θ)=2πr21000kmf1000kmsin2θ. (B5)

The azimuthal component of velocity is also obtained as follows:

 uφ(r,θ)=2πr21000kmrsinθf1000km. (B6)

As for , we simply assume that it is independent of at :

 σ(θ)=−√32π2GM˙MB1000km√r1000km. (B7)

Then the profile of the toroidal magnetic field is given as follows:

 Bφ(r,θ)=B1000km√r1000kmrsinθ, (B8)

where the is the strength of toroidal magnetic field at . Note that the profiles of and thus obtained apparently satisfy the boundary conditions on the symmetry axis and the equator in Equation (7).

Next we consider the case, in which the magnetic field has both poloidal and toroidal components. Then the situation is more complicated with and being no longer independent of each other (Fujisawa et al. 2013). From the continuity equation of magnetic field (), the magnetic flux function is defined to give the poloidal magnetic-field components as

 Br=1r2sinθ∂Ψ∂θ,Bθ=−1rsinθ∂Ψ∂r. (B9)

Since the poloidal magnetic-field lines are coincide with the streamlines in ideal MHD, the magnetic flux function is also constant along each streamline,

 Ψ≡Ψ(ψ). (B10)

Then and are expressed in terms of the magnetic flux function, which is in turn a function of the stream function alone (see Fujisawa et al. 2013) as

 ℓ(ψ) =rsinθuφ−rsinθdΨ(ψ)dψBφ, (B11a) σ(ψ) =Bφrsinθρ−4πrsinθuφdΨ(ψ)dψ. (B11b)

Note that the specific angular momentum () is no longer conserved along a streamline but the is sill and too remains a conserved quantity although it also contains a term that depends on .

Just as in the previous case, in which the purely toroidal magnetic field is considered, we set the functional form of outside the shock surface from the condition that the streamlines and hence the poloidal magnetic-field line are spherical as follows:

 Ψ =−r21000kmB0cosθ, (B12a) Br(r,θ) =r21000kmr2B0, (B12b) dΨdψ =−r21000kmB0˙M=constant, (B12c)

where is a constant corresponding to the strength of poloidal magnetic field at . It is apparent from Equation (B12b) that the radial magnetic field does not depend on outside the shock surface. Since the component of magnetic field vanishes outside the shock surface, the component of the ideal MHD condition given in Equation (2) becomes

 ∂∂r(rBruφ−rBφur)=0 (B13a) ⇒(rBruφ−rBφur)=C1(θ), (B13b) ⇒uφrsinθdΨdψ∂ψ∂θ−Bφ4πρrsinθ∂ψ∂θ=C1(θ). (B13c)

In the above equations, is a function of alone and is related with as follows (see Equation B11b):

 σ(θ)=−4π(∂ψ∂θ)−1C1(θ). (B14)

The component of the equations of motion (Equation 1) can be recast in a similar way as

 ur∂uφ∂r+uruφr−Br4πρ(∂Bφ∂r+Bφr)=0 (B15a) ⇒∂∂r(ruφ)−Br4πρur∂∂r(rBφ)=0 (B15b) ⇒ruφ−Br4πρur(rBφ)=C2(θ), (B15c) ⇒ruφ−rdΨdψBφ=C2(θ), (B15d)

where we use the fact that is constant in the current case. is again related with (see Equation B11a) as

 ℓ(θ)=sinθC2(θ). (B16)

Using and instead of and , we express and as

 Bφ =(rurBr−r3Br˙M)−1(C2−C1Br), (B17a) u