© ESO 2019

1 Introduction

Explosions of supernovae (SNe) are one of the most prominent events in the universe. They may provide excellent observational and theoretical material for studying a vast range of physics: from cosmological principles to particle physics. The interaction of the SN blast wave may serve as a probe of the progenitor mass loss and its circumstellar matter (CSM).

Most of the core-collapse (type II) SN progenitors are considered to be red supergiants (Smartt 2009, 2015; Pejcha & Prieto 2015; Müller et al. 2017), however, there is strong evidence that blue supergiants (BSG) can also be hydrogen-rich massive progenitors of type II SN events (e.g., Vanbeveren et al. 2013). Most famous example is SN 1987A, which in the time of explosion was a BSG star evolved from the previous stage of red supergiant; see, for example, Meynet et al. (2015). Furthermore, there are indications of a strong connection between BSG progenitors and superluminous SNe that are likely associated with very massive stars that retained their thick hydrogen envelopes and that explode within the dense CSM environment. There is also a strong evidence of rugged and anisotropic nature of the CSM, since some of these stars may have been similar to luminous blue variables that lost a large amount of mass, up to several M_⊙, prior to explosion (Smith et al. 2007, 2008; Ofek et al. 2013; Smith 2017).

The interaction of SNe with CSM may significantly power and modify the observed luminosity. After a radiation mediated shock goes through the stellar body of a SN progenitor, then it continues propagating through CSM. Depending on the optical thickness of CSM in the direction to the observer, the radiation is in the first phase more or less absorbed while an increasing amount of energy is accumulated in the shocked region (Svirski et al. 2012). After the optical depth of CSM drops, this energy is released as a breakout pulse that is much more energetic than a CSM-less breakout, however, the energy release in the presence of CSM happens over much longer timescales. Yet the energy behind the shock grows, if the CSM density does not fall abruptly, after this breakout, which is accompanied by the characteristic profile of the light curve (e.g., Chevalier & Irwin 2011; Svirski et al. 2012; Smith 2017).

It follows from the nature of (some) BSG stars that during their pre-explosion epochs a CSM may be formed that is anisotropic. In order to make a simulation of a hydrodynamic interaction of exploding SN progenitor with such aspherical CSM consisting of isotropic stellar wind component and a circumstellar equatorial disk, we consider some typical representatives of relevant types of SN progenitors that may be associated with BSG stars and with various types of circumstellar disks, such as Pop III stars, sgB[e] stars and luminous blue variables (Lee et al. 1991; Krtička et al. 2011; Kurfürst et al. 2014). We regard the structure of circumstellar disks for the purpose of modeling as that corresponding to viscous outflowing disks, which are typical for classical Be stars. However, we may assume to find this disk type in sgB[e] stars as well (for a review see, e.g., Zickgraf 1998; Hillier 2006), where a disk or rings of high density material have been detected (Kraus et al. 2013). Other types of dense equatorial disks or disk-like density enhancements may also be formed owing to, for example, magnetically compressed winds or binarity and accretion around certain classes of, for instance, B[e] stars (e.g., Hillier 2006), luminous blue variables (e.g., Schulte-Ladbeck et al. 1994; Davies et al. 2005), and post-AGB stars (e.g., Heger & Langer 1998).

Although there have been a number of studies regarding interactions of SNe expansion with various forms of spherically symmetric CSM, there is a lack of multidimensional hydrodynamic models that describe SNe interaction with dense aspherical CSM containing equatorial disks. The importance of examining such a CSM configuration is further enhanced by the fact that the disk densities can be significantly higher in the near-stellar region than the spherically symmetric wind densities in the case of the same stellar mass-loss rate Ṁ. This is because of their concentration only in the stellar equatorial area along with their much lower velocity of radial outflow. We calculate time dependent evolution of such SN–CSM interactions for three different values of pre-explosion mass-loss rate Ṁ, which led to different masses and densities of CSM and different ratios of isotropic wind and circumstellar disk densities. We use our models to explore and compare the profiles of basic hydrodynamic quantities under the significantly aspherical conditions. We ignore the initial density profile of the progenitor star and the effects of radiative cooling and other nonadiabatic processes in the models because our simulations cover only a relatively short post shock-emergence time.

2 Basic physics and parameterization of the circumstellar medium

We studied the interaction of a SNe blast wave with the CSM consisting in the equatorial disk and spherically symmetric CSM (stellar wind).

2.1 Spherically symmetric CSM

We assume a spherically symmetric CSM with a density profile given as a power law $${\rho }_{\text{sw}}={\rho }_{\text{0,sw}}{(R_{\star }/r)}^{w_{\text{sw}}},$$(1)

where r is radius in spherical coordinates and R_⋆ is the radius of the central object. We examined various values of spherically symmetric CSM base density ρ_0,sw as well as various values of spherically symmetric CSM density slope factor w_sw within the range from w_sw = 0 to w_sw = 5. However, the naturally expected structures of the spherically symmetric CSM that correspond to a steady mass loss Ṁ_sw = const. fulfill ρ_sw = Ṁ_sw∕(4πr²v_sw) ∝ r⁻² (Chevalier & Soker 1989; Chevalier & Fransson 1994; Chevalier & Irwin 2011; Moriya et al. 2013a, among others) where we input the average initial outflow velocity v_sw = 100 km s⁻¹ (Moriya et al. 2014). The adopted spherically symmetric CSM base density ranges from ρ_0,sw ≈ 10⁻¹⁰ kg m⁻³ (corresponding to stellar winds of OB stars; see, e.g., Kudritzki & Puls 2000) to ρ_0,sw ≈ 10⁻⁶ kg m⁻³, created, for example, by enhanced mass loss from violent convective motion caused by the unstable nuclear burning at latest stages before core collapse or by pre-explosion pulses of a progenitor prior to a SN event (cf. Moriya et al. 2014; Smith & Arnett 2014). Since the interaction of the SN ejecta with a spherically symmetric CSM is not a key problem in this study, we do not investigate finer details in this point.

2.2 Circumstellar equatorial disks structure

The disks around massive stars result either from the rotation of contracting matter and from the evacuation of its angular momentum (accretion disks, see, e.g., Pringle 1981; Frank et al. 2002; Maeder 2009) or they stem from the angular momentum loss from fast rotating central object (stellar decretion disks, e.g., Lee et al. 1991; Okazaki 2001; Krtička et al. 2011; Kurfürst et al. 2014).

In the following we briefly describe the equations that are essential for modeling the structure of a rotating Keplerian disk in this study: stationary, vertically integrated, and axisymmetric equation of continuity takes the cylindrical form $$\dot{M}=2\pi R\Sigma V_{\text{R}}=\text{\hspace{0.17em}const.,}$$(2)

where Ṁ is the (constant) disk mass-loss rate, R is radius in cylindrical coordinates, $\Sigma ={\displaystyle {\int }_{-\infty }^{\infty }\rho }\text{\hspace{0.17em}d}z$ is the disk surface density in 1D thin disk approximation, and V_R is the velocity of a radial flow of the disk matter. We hereafter distinguish between cylindrical (R, ϕ, z) and spherical (r, θ, ϕ) coordinates using uppercase and lowercase letters, i.e., radius and velocities in cylindrical coordinates are denoted using uppercase letters, while radius and velocities in spherical (and Cartesian) coordinates are denoted using lowercase letters.

Disk density, pressure, and temperature are determined in general by equations of radial momentum and angular momentum, which we present in their stationary form (e.g., Okazaki 2001; Krtička et al. 2011; Kurfürst et al. 2014), $$\begin{array}{l}{V}_{\text{R}}\frac{\partial V_{\text{R}}}{\partial R}-R{\Omega }^2+\frac{1}{\Sigma }\frac{\partial (a^2\Sigma )}{\partial R}-\mathcal{F}=0,\\ \frac{\dot{M}}{2\pi }\frac{\partial (R^2\Omega )}{\partial R}-\frac{\partial }{\partial R}\left({\alpha }_{\text{vis}}{a}^2{R}^2\Sigma \frac{\partial \ln \Omega }{\partial \ln R}\right)=0,\end{array}$$

where Ω is the angular velocity, a is the speed of sound, and α_vis is the parameter of viscosity (Shakura & Sunyaev 1973). The term $\mathcal{F}=-{\Sigma }^{-1}{\displaystyle {\int }_{-\infty }^{\infty }\rho }{g}_R\text{\hspace{0.17em}d}z$ is the density-weighted vertically integrated radial component of gravitational force per unit volume, ρg_R ≡ ρGM_⋆R∕r³, in the thin disk approximation, z ≪ R, where G is the gravitational constant and M_⋆ is the mass of a central object (star).

We assume a disk initial radial temperature profile for the purpose of the study in a form of the simple power law $$T=T_0{(R_{\star }/R)}^p,$$(5)

where T₀ is the temperature of the disk near the stellar surface (disk base temperature) and p is a free parameter (p ≥ 0) with theoretically estimated values in range between 0 and 0.4 (cf. Carciofi & Bjorkman 2008). For simplicity we omit the vertical (relatively shallow) variations of the temperature profile (e.g., Sigut & Jones 2007; Kurfürst et al. 2018). However, the temperature profile serves in our study only as a basis for determining the initial gas pressure of CSM (where ρ may be either ρ_sw or ρ_csd), $$P=\rho a^2,\text{\hspace{1em} where \hspace{1em}}{a}^2=kT/(\mu m_{\text{u}}),$$(6)

denoting k the Boltzmann constant, μ the mean molecular weight, and m_u the atomic mass unit. After doing various simulations using different indices 0 ≤ p≤ 0.4 with a minor impact on the effect of the counterpressure, we use in the simulations performed in the study the power-law index p = 0 (cf. also Sect. 5). For the same reason in the calculations we neglect the disk velocities V_R and Ω (which are negligible comparing to SN-CSM interaction dynamics) as well as the effects of gravitational force and viscosity.

For theanalytical approximation of the disk midplane density profile ρ_eq (R) we use the equation which relates ρ_eq to the vertically integrated density Σ(R) (e.g., Krtička et al. 2011), $$\Sigma =\sqrt{2\pi }{\rho }_{\text{eq}}H.$$(7)

Following the analytical approach of, for example, Okazaki (2001), Σ(R)~ R⁻² and by denoting H the disk vertical scale height, H = a∕Ω, in the Keplerian disk, we (analytically) obtain $${\rho }_{\text{eq}}(R)\approx {\rho }_{0,\text{csd}}{(R_{\star }/R)}^{3.5},$$(8)

where ρ_0,csd is the disk base density and H ~ R^3∕2 in the disk Keplerian region (Kurfürst et al. 2014). Numerical models usually find the power d of disk midplane radial density slope ρ_eq ~ R^−d between 3 and 4, and the disk base density ρ_0,csd is, for example, in case of classical Be stars, estimated between 10⁻¹⁰ kg m⁻³ to almost 10⁻⁶ kg m⁻³ (see, e.g., Granada et al. 2013). In case of much denser disks around, for instance, sgB[e]s, we assume disk base densities up to the values of on the order of 10⁻⁴ kg m⁻³ (cf. Kurfürst et al. 2018). The vertical profile of the thin circumstellar disk for z ≪ R, which is the case for disks up to the distance of several stellar radii from the central object, takes the Gaussian form (e.g., Pringle 1981; Okazaki 2001); this reflects the disk vertical hydrostatic balance (see also Fig. 1), $${\rho }_{\text{csd}}(R,z)={\rho }_{\text{eq}}(R){\text{e}}^{-\frac{G{M}_{\star }}{a^2{R}^3}\frac{z^2}{2}}.$$(9)

The initial structure of the disk is in this study basically parameterized by the disk mass-loss rate Ṁ. The (density independent) value of the disk base radial velocity V_R(R_⋆) is pre-calculated from a 1D disk model (see Kurfürst et al. 2014; Kurfürst 2015, for more details). The disk base surface density Σ_0,csd is fully determined from Ṁ and V_R(R_⋆) using the equation of mass conservation Eq. (2). The base surface density Σ_0,csd is subsequently converted to the base volume density ρ_0,csd using the Gaussian vertical hydrostatic equilibrium solution (see Eqs. (7) and (9)). The complete initial disk profile in pressure and density is obtained from Eqs. (5), (6), (8), and (9).

Fig. 1 Color map of the initial density of the SN progenitor with M⋆ = 45 M⊙, R⋆ = 80 R⊙ (with mean density ${\langle \rho \rangle }_{\star }\approx 0.125{\text{kg m}}^{-3}$) and CSM with assumed Ṁsw = 10−1 M⊙ yr−1 and Ṁcsd = 10−4 M⊙ yr−1 (corresponding to Model A in Sect. 5). The base density of the spherically symmetric part of CSM ρ0,sw ≈ 7 × 10−6 kg m−3, while the equatorial disk base density ρ0,csd ≈ 4 × 10−4 kg m−3 (see also Table 1 in Sect. 5).

3 Hydrodynamics and similarity solution of expanding supernova envelope

It is conventional to describe the evolution of a SN remnant in several consecutive phases (e.g., Matzner & McKee 1999). In the early phase the ejecta expands more or less ballistically into the surrounding medium. Ahead of the ejecta a strong shock is driven into the ambient medium while the high pressure behind the forward shock drives a reverse shock into the ejecta (Chevalier 1982; Chevalier & Soker 1989; Truelove & McKee 1999). This “ballistic” or free expansion phase ends when the total amount of material swept up by the forward shock has a mass comparable to that of the ejecta, or equivalently, the mean density of the ejecta drops to that of the surrounding medium.

The basic analytical solution of the SN explosive event follows the idea of homologous (u ∝ r∕t) energy-conserving hydrodynamic expansion, using the Sedov-Taylor blast wave solution (Sedov 1959) and the Sakurai solution of a shock wave in a nonuniform medium (Sakurai 1960) that has different density profiles of an internal gas and surrounding medium (Sedov 1977; Zel’dovich & Raizer 1967). Gravitational force is unimportant except for the very interior layers, and the radiative losses are almost negligible in the early phase of expansion when comparing to a kinetic energy of the gas. We thus may regard the process as adiabatic. The scaling parameters are the energy of explosion E_SN, mass M_ej of the expanding ejecta, and initial stellar radius R_⋆.

The main aim of our study is the modeling of interactions between the spherically expanding SN remnant and the surrounding CSM, which may be spherically symmetric (stellar winds and/or pre-explosion gaseous shells) or asymmetric (circumstellar disks, disk-like density enhancements, and bipolar lobes). Therefore, we neglect a pre-explosion SN interior density distribution, which is completely transformed during the explosion, and on the other hand is of little influence on the structure of the studied interaction. We also do not introduce the complete analytical formalism in this section, and we only briefly review the main relations in the way in which we use these for the comparative semi-analytical solution given in Appendix A.2. For the complete description see the fundamental solutions of, for example, Chevalier (1982) and Nadezhin (1985).

The basic hydrodynamic equations in an axially symmetric slice of spherical coordinates in the plane ϕ = 0, where the angle θ ranges therefore in full angle to include the plane ϕ = π, and where x = R sin θ, z = R cos θ, θ ∈⟨0, 2π⟩ (involving only the terms that are fundamental for nonviscous, non-self-gravitating, and axially symmetric (∂∕∂ϕ = 0) gas expansion), i.e., the continuity equation, radial and polar components of momentum equation, and energy equation, take the form $$\begin{array}{lll}\hfill & \hfill & \frac{\partial \rho }{\partial t}+\frac{1}{r}\frac{\partial }{\partial r}\left(r\rho u_r\right)+\frac{1}{r}\frac{\partial \left(\rho u_{\theta }\right)}{\partial \theta }=0,\hfill \\ \hfill & \hfill & \frac{\partial u_r}{\partial t}+u_r\frac{\partial u_r}{\partial r}+\frac{u_{\theta }}{r}\frac{\partial u_r}{\partial \theta }-\frac{u_{\theta }^2}{r}+\frac{1}{\rho }\frac{\partial P}{\partial r}-F_0=0,\hfill \\ \hfill & \hfill & \frac{\partial u_{\theta }}{\partial t}+u_r\frac{\partial u_{\theta }}{\partial r}+\frac{u_{\theta }}{r}\frac{\partial u_{\theta }}{\partial \theta }+\frac{u_r{u}_{\theta }}{r}+\frac{1}{\rho r}\frac{\partial P}{\partial \theta }=0,\hfill \\ \hfill & \hfill & \frac{\partial \left(P{\rho }^{-\gamma }\right)}{\partial t}+u_r\frac{\left(P{\rho }^{-\gamma }\right)}{\partial r}+\frac{u_{\theta }}{r}\frac{\left(P{\rho }^{-\gamma }\right)}{\partial \theta }=0,\hfill \end{array}$$

where ρ is the density, u_r is the radial component of the flow velocity, u_θ is the polar component of the flow velocity, P is the scalar pressure, F₀ is the body force per unit mass (typically gravitation), and γ is the ratio of specific heats. Equation (13) is obtained from the zero total derivative of the adiabatic transformation of a perfect gas Pρ^−γ = const. However, we further neglect the gravity term F₀ in the momentum Eq. (11), since gravitational force is of very little importance for the studied process.

In the reference frame that is co-moving with the forward shock (for simplicity the planar propagation of the very thin shocked layer) omitting viscosity and assuming constant adiabatic index γ, which is astrophysically relevant (Nadezhin 1985), we can rewrite the adiabatic hydrodynamic Eqs. (10)–(13) on both sides of the shock front into the simple conservative form (Zel’dovich & Raizer 1967), $$\begin{array}{l}{\rho }_1{u}_1={\rho }_0{u}_0,\\ {\rho }_1{u}_1^2+P_1={\rho }_0{u}_0^2+P_0,\\ \frac{\gamma }{\gamma -1}\frac{P_1}{{\rho }_1}+\frac{u_1^2}{2}=\frac{\gamma }{\gamma -1}\frac{P_0}{{\rho }_0}+\frac{u_0^2}{2},\end{array}$$

where the upstream and downstream quantities are labeled with subscripts 0 and 1, respectively. The same velocities in the laboratory frame are D − u₀ and D − u₁, where we denote D the propagation speed of the forward shock. To get the analytical solution, in the early phase after explosion we approximate the process as a free expansion with the internal velocity u of the gas proportional to radius. Assuming the density profile ρ_in in the envelope as a (time-dependent) power law and the density of the surrounding (initially static) medium ρ_out as an alternate power law, we obtain $$u=\frac{r}{t},{\rho }_{\text{in}}(t)=A{r}^{-n}{t}^{n-3},{\rho }_{\text{out}}=B{r}^{-w}.$$(17)

Assuming naturally the strong shock condition, P_out ≪ P_in, the proportionality for the inner pressure, P_in ~ ρ_inu² (in the shock co-moving frame, or alternatively P_in ~ ρ_outD² in the laboratory frame; see Zel’dovich & Raizer 1967), gives the inner and outer pressure in the form $$P_{\text{in}}(t)=\tilde{A}{r}^{2-n}{t}^{n-5},P_{\text{out}}=\tilde{B}{\rho }_{\text{out}}^{\gamma }.$$(18)

We describe the formalism of the analytical similarity solution based on Rankine–Hugoniot relations (see the principles given in Chevalier 1982; Nadezhin 1985) in the Appendix A.1, while we describe the principles of semi-analytical model based on analytical solution in Appendix A.2. We also describe in Appendix A.3 the analytical principles of the fundamental limitation of the adiabatic self-similar solution for only a certain time interval. The specific values of this limitation are given within the sections corresponding to the particular models.

4 Numerical method

For the modeling of the interaction with extremely sharp discontinuities characterized by extremely high pressure gradients we have developed and used our own version of single-step (unsplit, ATHENA-like, e.g., Stone et al. 2008) finite volume Eulerian hydrodynamic code based on Roe’s method (Roe 1981; Toro 1999; Kurfürst et al. 2018, see also the description of our geometrical modifications of the algorithm used for different astrophysical situations in Kurfürst & Krtička 2017).

For the time-dependent calculations we write the adiabatic hydrodynamic Eqs. (10)–(13) in two different manners. The primitive form generally is (see, e.g., Stone et al. 2008) $$\frac{\partial W}{\partial t}+\nabla \cdot F(W)=0,$$(19)

where the primitive variables are W = ρ, u, P, and the fluxes F(W) = ρu, ρu | u + P 1, ρuH, where | denotes the dyadic product of two vectors, H = (E + P)∕ρ is the enthalpy (+ specific kinetic energy), 1 is the unit matrix, and E = ρɛ + ρu²∕2 is total energy (while ɛ is specific internal energy and $u^2=u_r^2+u_{\theta }^2$). The primitive (as well as conservative) form of continuity equation is given by Eq. (10), while Eqs. (11) and (12) represent the primitive form of radial and polar components of momentum equation. We transform Eq. (13) into explicit primitive form of energy equation by inserting Eq. (10), $$\frac{\partial P}{\partial t}+u_r\frac{\partial P}{\partial r}+\frac{u_{\theta }}{r}\frac{\partial P}{\partial \theta }+\gamma P\left(\frac{\partial u_r}{\partial r}+\frac{u_r}{r}+\frac{1}{r}\frac{\partial u_{\theta }}{\partial \theta }\right)=0,$$(20)

where γ = 4∕3 for the radiation dominated gas. Equations (10)–(12) together with Eq. (20) are used in the primitive numerical scheme (Eq. (19)).

The analogous conservative form is (see, e.g., Norman & Winkler 1986; Hirsch 1989; Stone & Norman 1992; Feldmeier 1995; LeVeque et al. 1998) $$\frac{\partial U}{\partial t}+\nabla \cdot F(U)=0,$$(21)

where the conservative variables are U = ρ, M, E and F(U) = ρu, M | u + P 1, MH (where M = ρu is the momentum density) for the mass, momentum, and energy equations, respectively. Abbreviating γ′ = γ∕(γ − 1), the explicit adiabatic form of the enthalpy H used in the flux terms is H = γ′P∕ρ + u²∕2. The explicit conservative form of radial (omitting the gravitational force; see Sect. 3) and polar components of momentum equation and of the energy equation, respectively, is $$\begin{array}{l}\frac{\partial M_r}{\partial t}+\frac{1}{r}\frac{\partial }{\partial r}\left(r{M}_r{u}_r\right)+\frac{1}{r}\frac{\partial (M_{\theta }{u}_r)}{\partial \theta }-\rho \frac{u_{\theta }^2}{r}+\frac{\partial P}{\partial r}=0,\\ \frac{\partial M_{\theta }}{\partial t}+\frac{1}{r}\frac{\partial }{\partial r}\left(r{M}_r{u}_{\theta }\right)+\frac{1}{r}\frac{\partial (M_{\theta }{u}_{\theta })}{\partial \theta }+\rho \frac{u_r{u}_{\theta }}{r}+\frac{1}{r}\frac{\partial P}{\partial \theta }=0,\\ \frac{\partial E}{\partial t}+\frac{1}{r}\frac{\partial }{\partial r}\left(r{M}_{r}H\right)+\frac{1}{r}\frac{\partial }{\partial \theta }\left(M_{\theta }H\right)=0.\end{array}$$

Equation (10) together with Eqs. (22)–(24) are used in the conservative numerical scheme (Eq. (21)).

The polar form of Roe matrices A_W = ∂F(W)∕∂W and A_U = ∂F(U)∕∂U, resulting from the linearization of Eqs. (19) and (21), is in this case identical and is supplemented by formalism introduced in Stone et al. (2008) and Skinner & Ostriker (2010). We used the (fully capable) second-order piecewise linear reconstruction algorithm, without any necessity to alternate the Roe’s method using the HLLE (Harten, Lax, van Leer, and Einfeldt; cf. Einfeldt et al. 1991; Stone et al. 2008) solver or similar positive-definite solvers in case of unexpected density or pressure drop to negative values in the intermediate states of computation (see Stone et al. 2008, for the details), even in extremely thin zones of very strong shocks.

To derive the initial conditions for hydrodynamic quantities at the explosion time t₀ we assumed that the interior shock wave basically rearranges the original density structure of the progenitor. The density gradient between the inner and outer interstellar region is reduced by the rapid expansion of the inner regions (Arnett 1996, see also, e.g., Fig. 1 in Arnould & Prantzos 1986 showing the evolution of the density structure from core collapse to mantle ejection). We therefore neglect the details in a SN interior density structure that can barely affect the main features of examined effects in a complex interaction between SNe and CSM. We thus assume for simplicity that the stellar density is equal to the “averaged” homogeneous initial density ${\langle \rho \rangle }_{\text{}\star }=3M_{\star }/(4\pi R_{\star }^3)$ of an exploding star (according to the stellar parameters introduced in Sects. 3, 5, and Appendix A.2) in the range 0.1 − 1 kg m⁻³ (${\langle \rho \rangle }_{\text{}\star }\approx 0.125\text{\hspace{0.17em}kg\hspace{0.17em}}{\text{m}}^{-3}$; cf. Fig. 1).

The injected SN explosion energy is E_SN = 10⁴⁴ J. A SN explosion drives a radiation dominated shock, i.e., the internal energy inside the sphere of shock wave that propagates through the star is dominated by radiation (Nakar & Sari 2010). For that reason we assume the averaged homogeneous initial internal pressure of a photon gas; we consider the whole stellar volume as a “thermal bomb” and neglect the details of the processes connected with the shock wave propagation through the star which, in fact, have a secondary effect on the studied interaction of SN with the aspherical CSM. The averaged homogeneous initial internal pressure of a photon of gas takes the form, ${\langle P\rangle }_{\text{}\star }=E_{\text{SN}}/(3V_{\text{}\star })=E_{\text{SN}}/(4\pi R_{\star }^3)$, i.e., in the range between 10¹⁰ and 10¹¹ Pa (${\langle P\rangle }_{\text{}\star }\approx 4.63\times {10}^{10}\text{\hspace{0.17em}Pa}$).

We insert the initial state for the CSM in the following explicit form: the density of a spherically symmetric CSM is expressed simply by Eq. (1) while its outflow velocity v_sw = 100 km s⁻¹ (see Sect. 2.1 and also Moriya et al. 2014) is used only for initial density estimate and is neglected for the calculation itself. We assume an isothermal initial temperature structure of the CSM (of both the components) within the computational domain, given as T_sw = T_csd = 0.7T_eff (cf. Carciofi & Bjorkman 2008; Kurfürst et al. 2018). This also indicates the initial counterpressure of the CSM as P = ρa², where the square sound speed a² is given by Eq. (6). The initial circumstellar equatorial disk density structure (Eqs. (8) and (9)) is in the polar coordinates specified in the following explicit form as $${\rho }_{\text{csd}}(r,\theta )={\rho }_{0,\text{csd}}{\left(\frac{R_{\star }}{r\left|\text{sin\hspace{0.17em}}\theta \right|}\right)}^{3.5}\exp \left(-\frac{G{M}_{\star }}{2a^{2}r}\frac{{\text{cos}}^2\theta }{\left|{\text{sin}}^3\theta \right|}\right).$$(25)

We neglect the velocity field of the circumstellar equatorial disk (see Sect. 2.2), since we regard its effects during the studied process as negligible. The total CSM density profile is given as ρ_sw (Eq. (1)) + ρ_csd (Eq. (25)).

Within our considerations we may set R₀ = R_⋆, while t₀ is certainly less than one day, for example, for SN 1987A the time t₀ is estimated to be 53 min after core collapse (Arnett 1996). We set the free boundary conditions at the inner (originally interstellar region) as well as at the outer boundary of the computational domain during the expansion phase. The values for the grid boundary and ghost zones are extrapolated from mesh interior values as a 0th order extrapolation.

We performed all the calculations on a center-symmetric area domain (2D axially symmetric slice of the spherically symmetric 3D domain) with a radius of 20 R_⋆, and with a numerical polar grid uniformly scaled in both directions r and θ. The star is located in the center of the computational domain. We also tested the Cartesian uniform grid, interlaced by an equatorial plane of the spherical coordinate system, for comparison with the calculated model; however, we achieved smoother profiles and the optimal computational cost using the polar grid. The number of spatial grid cells was 400 in radial and 300 in azimuthal directions for the models performed up to radius 20 R_⋆ and in 2π polar domain. The complete physical time of the simulations corresponds in this case to approximately 165 h. We also performed the calculation on 2D grid containing one polar quadrant (π∕2 polar domain) with 2400 grid cells in radial and 480 grid cells in azimuthal direction with the complete physical time of the simulations corresponding to approximately 180 h. The selected grid aspect ratio has a stabilizing effect because the cells that are stretched along the tangent to the shock have a damping effect on the development of the carbuncle numerical instability (Pandolfi & D’Ambrosio 2001). All the models are described in Sect. 5 (see also Kurfürst et al. 2017, 2018; Kurfürst & Krtička 2017).

5 Numerical models

We selected a BSG star with M_⋆ = 45 M_⊙, R_⋆ = 80 R_⊙, T_eff = 25 000 K as the progenitor for all studied models; all of the parameters of this star are described in Sect. 4 (see also Fig. 1). The assumed properties of CSM of studied models (denoted as Model A, B, and C), i.e., the adopted mass-loss rates Ṁ and the approximate base densities ρ_0,sw and ρ_0,csd derived from the mass conservation equation (see Sects. 2.1 and 2.2), are given in Table 1. In all models we assume v_sw = 100 km s⁻¹ (Moriya et al. 2014, see also Sect. 4), which is an order of magnitude latitudinally averaged estimate considering the effects of possible high rotation (Maeder & Meynet 2000; Kraus et al. 2008, 2010) and high radiation pressure.

5.1 Model A: high density of the surrounding media

We first study the SN–CSM interaction for the case of high density of surrounding media. The CSM consists in the spherically symmetric component and the aspherical circumstellar disk described in Sect. 2. The parameters of the CSM and progenitor star are given in Table 1 and in the first paragraph of Sect. 5.

Integrating the density profiles, we estimate the total mass of the CSM components: for spherically symmetric component of CSM we obtain the total mass $M_{\text{sw}}=4\pi R_{\star }^2{\rho }_{0,\text{sw}}{\displaystyle {\int }_{R_{\star }}^{r_1}\text{d}}r\approx 0.15M_{\odot }$ within the radius r₁ = 20 R_⋆ (within the radius of the computational domain), while we obtain the total mass for circumstellar disk component within the same domain as $M_{\text{csd}}=2\pi R_{\star }^2{\Sigma }_{0,\text{csd}}{\displaystyle {\int }_{R_{\star }}^{r_1}{R}^{-1}}\text{\hspace{0.17em}d}R\approx 0.65M_{\odot }$ (cf. Krtička et al. 2011; Kurfürst et al. 2014, 2018).

We show in Figs. 2 and 3 the 2D snapshots of density, velocity, and temperature, in equatorial (x)–polar (z) plane, at selected times t = 16 h (only density) and t = 66 h. We also add illustratively separated 1D slices of the graphs of selected quantities in equatorial plane and in polar direction in Fig. 4 at the time t ≈ 66 h, supplemented by a graph of specific radiation entropy s = 16σT³∕(3cρ), where σ is the Stefan-Boltzmann constant and c is the speed of light (Zel’dovich & Raizer 1967). The figures show the regions of decelerated expansion velocity toward the equatorial disk while the expansion is unlimited (relatively, according to the density of the spherically symmetric CSM component) in the polar direction. All the characteristics show enhanced peaks in the SN-disk contact region and indicate shoulders and strips of overdensity and enhanced temperature that propagate around the equatorial disk.

We show in the upper panel of Fig. 5 the density slopes n and w (see Eqs. (17) and (A.1)) and in the lower panel of Fig. 5 the temperature slopes p (cf. Eq. (5)), calculated in the numerical model at the same time t ≈ 66 h since shock emergence. However, the outer shock, which roughly separates the outer region (denoted in Sect. 3 by subscript“out”) from all the other regions between the two shock waves, moves to the radius approximately 11 R_⋆ in the equatorial plane and 16.5 R_⋆ in the polardirection. The inner envelope density slope ranges in the polar direction between n = 0 and n ≈ 9, while in the equatorial plane the density slope is highly discontinuous and becomes much steeper. However, we cannot simply apply the analytical similarity relations in the equatorial plane with respect to the constraint w <3 given by Eq. (A.21). The temperature slope in the smooth region of the inner envelope (excluding the areas of sharp discontinuities) ranges between p = 0 and p ≈ 2.7 in the equatorial and in the polar direction. In order to avoid inadequate computational difficulties, we calculate the temperature in all the domain as radiative (T ~ P^0.25; the same applies for the other models). We may regard the displayed distances reached by the forward shock wave in the equatorial and in the polar direction as a certain scale of the level of asphericity (except the slopes of the selected quantities) of the SN explosion; the same also applies for another models (see Figs. 11 and 18).

In Figs. 6 and 7 we demonstrate the semi-analytically found (see Sect. A.2) similarity solution of density, velocity, pressure, and temperature for two inner density slopes n = 7 and n = 9 that correspond to the average and maximum numerically calculated inner density slopes in the polar direction where the similarity solution holds and where the parameter w < 3 (see Eq. (A.1) and its consequences). In the disk-residing equatorial plane the analytical similarity solution is however violated because of the disk equatorial plane density profile ρ_eq (R) ~ R^−3.5 (w > 3).

Comparing slices of 2D numerical calculations in Fig. 4 with the 1D semi-analytical solution (cf. Sect. A.2) in polar direction in Figs. 6 and 7, we find that the solutions well agree in global parameters. That is, density ρ decreases from the value ρ ≈ 5 × 10⁻⁶ kg m⁻³ at the radius r = 11.5 R_⋆ to ρ ≈ 3 × 10⁻⁷ kg m⁻³ behind the outer shock, while it drops to almost 10⁻⁸ kg m⁻³ outside the outer shock; the radial (expansion) velocity u grows in the same domain from u ≈ 3000 km s⁻¹ to approximately 3500 km s⁻¹ behind the shock wave zone. We hereafter denote the radial velocity of expansion u for simplicity; its value in fact does not much differ from the total magnitude of velocity since the nonradial velocity components are more than an order of magnitude smaller. The pressure approaches P ≈ 10⁶ Pa, while the corresponding temperature T is slightly higher than 10⁵ K (compare Figs. 6 and 7 with Fig. 4). The semi-analytical solution is directly compared with 1D slices of our 2D numerical solution for the same time of calculation. We also add the 1D synchronous profile calculated using the SNEC-1.01 code (Morozova et al. 2015), which takes into account more realistic shape of the density enhancement zone between reverse shock and contact discontinuity (between 13.5 and 13.8 r∕R_⋆) that is formed during propagation of the shock wave through the stellar body (Sedov 1959; Sakurai 1960). Since our current models use the whole stellar body as an initial thermal bomb, the shock wave directly propagates (like a Riemann-Sod shock wave) to the region of dropped density. However, this detail in the shock wave structure has only a minor impact on the global solution within the study. For the SNEC calculation we input the same initial profile of the progenitor star and the CSM as for all the other calculations with radiation included, the initial thermal bomb zone (Morozova et al. 2015) is however limited only to a relatively small region near the stellar center.

We also constructed similar comparative diagrams for other values of parameter n lower than n = 7. However, in this case the radial distance between the inner and outer shock waves becomes too large and the values of analytically calculated quantities differ significantly from SNEC and from our 2D numerical result. In any case, we found that the analytical versus numerical solutions are closest (relatively) at this stage for the internal density slope parameter n ≈ 7 to n ≈ 9.

We also checked the constraints on the applicability of the adiabatic self-similar solution within the semi-analytical calculation (see Appendix A.3). For the model with high density (Sect. 5.1) and the density slope n = 7 we found the upper limit t_max ≈ 20 yr, while the lower limit t_min is on the order of several minutes at the start of the calculation. The lower limit continuously grows during the time evolution, however, the values of t_min remain several orders of magnitude lower than the actual time of simulation.

Figure 3 indicates that the zone of shear between the slower expansion into the region of circumstellar equatorial disk, where the vertical motion occurs owing to the violation of disk hydrostatic equilibrium caused by the interaction of the two masses. Figure 3 also shows that the spherically symmetric fast SN expansion may provoke the development of Kelvin-Helmholtz instabilities. This is further demonstrated in Fig. 8, where we show the snapshot of 2D density and expansion velocity field in later time t ≈ 165 h (up to the radius 50 R_⋆) and in higher numerical resolution. We thus get a comparison with the model in Fig. 3, where we can distinguish the details of evolution of the selected quantities in yet more evolved and more structured stage. The shoulders of overdensity and significantly decelerated expansion velocity at the zone of shear between the disk and spherically symmetric CSM are illustrated in a more detailed picture. There is also a clear indication of developing Kelvin-Helmholtz and Rayleigh-Taylor instabilities along the zone of shear interaction and SN-disk contact region.

Fig. 2 Model A: color map of the density structure of the interaction of SN (with progenitor and CSM parameters introduced in Sect. 3, Table 1, and in Fig. 1) with asymmetric CSM that forms a dense equatorial disk, at time t ≈ 16 h since shock emergence. Contours denote the densities ρ = 10−7, 10−6, 10−5, and 10−4 kg m−3.

Fig. 3 Model A: interaction of SN with asymmetric CSM at time t ≈ 66 h since shock emergence. Upper panel: color map of the density structure. Contours denote the densities ρ = 10−6, 10−5, and 10−4 kg m−3. Middle panel: color map of the velocity structure. Contours denote the radial (expansion) velocities u = 1000, 2000, and 3000 km s−1. Lower panel: color map of the temperature structure. Contours denote the temperatures T = 105, 5 × 105, and 106 K. Characteristic 1D sections of the model in the equatorial and polar plane are shown in Fig. 4.

Fig. 4 Model A: one-dimensional slices of variations of hydrodynamical variables in the equatorial plane (x-coordinate, solid line) and in the polar direction (z-coordinate, dashed line), at the same time as in Fig. 3. Panel A: density ρ. Panel B: expansion velocity u. Panel C: pressure P and temperature T. Panel D: specific entropy s.

Fig. 5 Model A: upper panel: slope parameters n and w of the density (see Eq. (17)) in the equatorial (solid line) and polar (dashed line) direction corresponding to the panel A in Fig. 4, at the same time. The position of the outer shock wave is at r ≈ 11 R⋆ in the equatorial direction, while it is at r ≈ 16.5 R⋆ in the polar direction. There are preserved the initial density slope parameters of CSM, wsw = 2 in the polar direction and wcsd = 3.5 in the equatorial direction, above these radii. The inner envelope density slope parameter increases from n = 0 to n ≈ 9 in the polar direction. The numbers nmax and nmin above and below the graph denote the maximal and minimal slopes within the peaks of the density gradient discontinuities. Lower panel: slope parameters p of the temperature (cf. Eq. (5)) in the equatorial (solid line) and polar (dashed line) direction corresponding to the temperature graph (panel C) in Fig. 4 at the same time. The inner envelope (smooth) temperature slope parameter increases to p ≈ 2.7 at r ≈ 7.8 R⋆ in the equatorial direction while it increases to approximately the same value at r ≈ 11 R⋆ in the polar direction. The slope parameter p outside the outer shock is p = 0.875 in the equatorial and p = 0.5 in the polar direction. The numbers pmax above the graph denote the maximal temperature slopes within the peaks of the T gradient discontinuities.

Fig. 6 Model A: results of semi-analytical solution of the density ρ, velocity u, pressure P, and temperature T (black line), with the slope parameters n = 7 and wsw = 2, in the polar direction. The method of calculation is described in Sect. A.2. The time is the same as in Figs. 3–5. The velocity at the outer shock wave region (at approx. 16 R⋆) correspondsto the polar expansion velocity in the panel B of Fig. 4 while the other quantities and semi-analytical model input variables are described in Appendix A.2. The blue line depicts the corresponding calculation from our 2D model. The corresponding 1D profile calculated using the SNEC-1.01 code is depicted with a red line (see the description in Sect. 5.1).

Fig. 7 Model A: as in Fig. 6, with the slope parameters n = 9 and wsw = 2, in the polar direction (the same method of calculation as in Fig. 6). The semi-analytical profiles are compared to the same polar expansion velocity as in Figs. 4 and 6.

5.2 Model B: intermediate density

We further studied interaction of an SN blast wave with the CSM for the case of intermediate density (Model B). The progenitor and CSM parameters are introduced in the first paragraph of Sect. 5 and in Table 1. Integrating the CSM density profiles, we obtain the following total mass of the CSM components: the total mass M_sw ≈ 0.015 M_⊙ within the radius of the computational domain r₁ = 20 R_⋆ and M_csd ≈ 0.65 M_⊙ (cf. Krtička et al. 2011; Kurfürst et al. 2014, 2018) within the same domain.

We show (analogous to Sect. 5.1) in Fig. 9 the 2D snapshots of density profiles in equatorial (x)–polar (z) plane at selected times t = 33 h and t = 66 h, respectively.The lower panel in Fig. 9 shows significant layers of overdensity that are formed near the zone where the expanding matter is forced to move along the interface between the spherically symmetric component of CSM and the denser circumstellar equatorial disk region of CSM. The development of Kelvin-Helmholtz instabilities is also remarkably indicated near the shear zones. We also add in Fig. 10 the separated 1D slices of graphs of density, velocity, pressure, temperature, and specific radiation entropy in equatorial and polar plane at the time t = 33 h. The model shows higher difference in expansion velocity between the equatorial and polar direction than in the previous model because of the lower density (lower mass injection) of spherically symmetric CSM. The front of the polar SN expansion reaches the velocity approximately 4000 km s⁻¹ while it approaches 2000 km s⁻¹ in the equatorial direction. Because the expansion shock front in the polar direction is too close to the limit of computational area in the time for which we plotted the previous models (t ≈ 66 h), we provide for this model the profiles of most of the characteristics in earlier time (t ≈ 33 h) while, on the other hand, the graphs of the previous model with high density would be in this time poorly illustrative.

We also show (analogously to Sect. 5.1) in the upper panel of Fig. 11 the density slopes n and w, while in the lower panel we show the temperature slopes p, calculated at the same time t ≈ 33 h since shock emergence. The outer shock moves at the radius of approximately 5.5 R_⋆ in the equatorial and 11 R_⋆ in the polar direction. The inner envelope density slope ranges in the polar direction between n = 0 and n ≈ 16, the temperature slope in the smooth region of the inner envelope ranges between p = 0 and p ≈ 1.9 in the equatorial direction and between p = 0 and p ≈ 3.6 in the polar direction.

The calculation of the constraints of applicability of the adiabatic self-similar solution (see Appendix A.3) gives for the model with intermediate density the upper limit t_max ≈ 21 yr for the density slope n = 7, while for n = 12 it increases to t_max ≈ 35 yr. For the lower limit t_min applies the same as in previous model with high density in Sect. 5.1.

Analogous to the previous model, we show in Fig. 12 the snapshots of 2D density, radial expansion velocity, and temperature in time t ≈ 180 h since shock emergence, up to the radius 50 R_⋆, in high numerical resolution, performed on 2D grid with 2400 grid cells in radial and 480 grid cells in azimuthal direction. Comparing this with Fig. 9, we distinguish the details of evolution of the selected quantities, including the even more significant shoulders of overdensity and over-heated gas near disk – spherically symmetric CSM shear zone. Although the forward shock front is in the snapshot already out of the computational domain, it particularly well illustrates the development of Kelvin-Helmholtz and Rayleigh-Taylor instabilities which is clearly indicated.

We also add in this intermediate density model the 1D equatorial, polar slices of density and expansion velocity (Fig. 13), and the temperature and specific entropy (Fig. 14), which correspond to the instant time shown in Fig. 12. A qualitative comparison with the profiles in Fig. 10 (in less advanced time) shows the decrease of density, however in a (roughly) self-similar mode. This is accompanied by an increase of the equatorial velocity peak. The peak of the polar velocity is currently already outside the computational domain, nevertheless the results of the calculation before it reaches the outer boundary show that the polar velocity increases very little (cf. the increase of polar velocity in the model with high density (Sect. 5.1) by comparing the values of the polar velocity maximum in Figs. 3 and 4 with Fig. 8).

Fig. 8 Model A: detailed color maps of the density (upper panel) and the radial velocity (lower panel) structure, up to the distance 50 R⋆ at time t ≈ 165 h since shock emergence. Contours denote the densities ρ = 10−7, 10−6, and 10−5 kg m−3 and the velocities 2500 and 3000 km s−1. The resolution of the simulation is 2400∕480 grid cells in the radial/azimuthal direction within the quadrant.

Fig. 9 Model B: upper panel: color map of the density structure at time t ≈ 33 h since shock emergence. Contours denote the densities ρ = 10−8, 10−7, 10−6, 10−5, and 10−4 kg m−3. Characteristic 1D sections of the model in the equatorial and polar plane are shown in Fig. 10. Lower panel: color map of the density structure at time t ≈ 66 h since shock emergence. Contours denote the densities ρ = 10−6, 10−5, and 10−4 kg m−3.

Fig. 10 Model B: 1D slices of hydrodynamical variables in the equatorial plane (x-coordinate, solid line) and polar plane (z-coordinate, dashed line) at time t ≈ 33 h since shock emergence (as in Fig. 9). Panel A: density ρ. Panel B: expansion velocity u. Panel C: pressure P and temperature T. Panel D: specific entropy s.

Fig. 11 Model B: upper panel: slope parameters n and w of the density (analogous to Fig. 5) in the equatorial (solid line) and polar (dashed line) direction, corresponding to panel A in Fig. 10 at the same time t ≈ 33 h. The positionof the outer shock wave is at r ≈ 5.5 R⋆ in the equatorial direction, while it is at r ≈ 11 R⋆ in the polar direction. The inner envelope density slope parameter increases from n = 0 to n ≈ 16 in the polar direction. Lower panel: slope parameters p of thetemperature in the equatorial (solid line) and polar (dashed line) direction corresponding to the temperature graph (panel C) in Fig. 10 at the same time. The inner envelope (smooth) temperature slope parameter increases to p ≈ 1.9 at r ≈ 4.1, R⋆ in the equatorial direction while it increases to p ≈ 3.6 at r ≈ 6.6 R⋆ in the polar direction.

Fig. 12 Model B: detailed color maps of the density (upper panel), radial velocity (middle panel), and temperature (lower panel) structure up to the distance 50 R⋆ at time t ≈ 180 h since shock emergence, with the same numerical resolution as in Fig. 8. Contours denote the densities ρ = 10−7, 10−6, 10−5, and 10−4 kg m−3, the velocities 2500, 3000, 3500, and 4000 km s−1, and the temperatures 7 × 104, 8 × 104, 105, and 5 × 105 K.

Fig. 13 Model B: 1D slices of density (left panel) and velocity (right panel)in the equatorial plane (x-coordinate, solid line) and in polar direction (z-coordinate, dashed line) at time t ≈ 180 h since shock emergence (corresponding to Fig. 12). The maximum of the polar velocity, which is already outside the computational domain, is approximately 4400 km s−1.

Fig. 14 Model B: 1D slices of temperature (left panel) and specific entropy (right panel) in the equatorial (x-coordinate, solid line) and in polar direction (z-coordinate, dashed line)at the same time as in Fig. 13.

5.3 Model C: low density

The asymetry of SN–CSM interaction zone is significant even in the case of low density CSM (Model C). The progenitor and CSM parameters are introduced in the first paragraph of Sect. 5 and in Table 1. The total masses of the spherical wind and circumstellar equatorial disk components within the radius of the computational domain r₁ = 20 R_⋆ are M_sw ≈ 1.5 × 10⁻⁴ M_⊙ and M_csd ≈ 6.5 × 10⁻³ M_⊙.

We show (analogous to Sects. 5.1 and 5.2) in Figs. 15 and 16 the 2D snapshots of density, velocity, and temperature, in equatorial (x)–polar (z) plane, at selected times t = 33 h and t = 66 h (only density). We again also add in Fig. 17 the separated 1D slices of graphs of density, velocity, pressure, temperature, and specific radiation entropy in equatorial plane and in polar direction at the time t = 33 h. The figures in this case show the regions of significantly decelerated SN expansion velocity toward the equatorial disk in greater distance from the star than in the previous cases; in the corresponding time t ≈ 66 h it is approximately 13 R_⋆ (see Fig. 16), while it is approximately 11 R_⋆ in the previous models (cf. Figs. 4 and 10). This is because the higher absolute value of expansion velocity even in this direction of maximum deceleration compared to previous models (cf. also Figs. 3, 9, and 17). The front of the polar SN expansion propagates with much higher velocity, i.e., > 7000 km s⁻¹ while it is ≈ 4000 km s⁻¹ in the model with intermediatedensity and ≈3000 km s⁻¹ in the model with high density in the same direction at the same time; the front should be within the same time t ≈ 66 h obviously in much larger distance than in the previous models. This is why, outside the computational area, we plot in this model the profiles of most of the characteristics in time t ≈ 33 h.

We show (analogous to Sects. 5.1 and 5.2) in the upper panel of Fig. 18 the density slopes n and w and in the lower panel the temperature slopes p, calculated at the same time t ≈ 33 h since shock emergence. The outer shock moves at the radius of approximately 9 R_⋆ in the equatorial plane and 18 R_⋆ in the polar direction. The inner envelope density slope ranges in the polar direction between n = 0 and n ≈ 17; the temperature slope in the smooth region of the inner envelope ranges between p = 0 and p ≈ 3.6 in the equatorial and also in the polar direction.

The calculation of the applicability limits of the adiabatic self-similar solution (see Appendix A.3) gives for the model with low density the upper limit t_max ≈ 190 yr for the density slope n = 7, while for n = 12 it gives t_max ≈ 3 × 10⁴ yr. The lower limit t_min is still practically unimportant (cf. the model with high density in Sect. 5.1).

As in both the previous models, we show in Fig. 19 the snapshots of 2D density and expansion velocity in a more evolved stage (in time t ≈ 105 h since shock emergence), up to the radius 50 R_⋆, performed on 2D grid with 2400 grid cells in radial and 480 grid cells in azimuthal direction. Comparing this with Fig. 15, we distinguish again the details of evolution of the selected quantities within the model.

Fig. 15 Model C: color maps of the hydrodynamical variables at time t ≈ 33 h since shock emergence. Upper panel: color map of the density structure. Contours denote the densities ρ =10−8, 10−7, 10−6, 10−5, and 10−4 kg m−3. Middle panel: color map of the velocity structure. Contours denote the velocities u = 5000 and 7000 km s−1. Lower panel: color map of the temperature structure. Contours denote the temperatures T = 105 and 106 K. Characteristic 1D sections of the model in the equatorial plane and in the polar direction are shown in Fig. 17.

Fig. 16 Model C: color map of the density structure of the interaction at time t ≈ 66 h since shock emergence. Contours denote the densities ρ = 10−6, 10−5, and 10−4 kg m−3.

Fig. 17 Model C: 1D slices of hydrodynamical variables in the equatorial (x-coordinate, solid line) and in polar direction (z-coordinate, dashed line) at time t ≈ 33 h since shock emergence (as in Fig. 15). Panel A: density ρ. Panel B: expansion velocity u. Panel C: pressure P and temperature T. Panel D: specific entropy s.

Fig. 18 Model C: upper panel: slope parameters n and w of the density (analogous to Figs. 5 and 11) in the equatorial (solid line) and polar (dashed line) direction, corresponding to panel A in Fig. 17 at the same time t ≈ 33 h. The positionof the outer shock wave is at r ≈ 9 R⋆ in the equatorial direction while it is at r ≈ 18 R⋆ in the polar direction. The inner envelope density slope parameter increases from n = 0 to n ≈ 17 in the polar direction. Lower panel: slope parameters p of thetemperature in the equatorial (solid line) and polar (dashed line) direction corresponding to the temperature graph (panel C) in Fig. 17 at the same time. The inner envelope (smooth) temperature slope parameter increases to p ≈ 3.6 at r ≈ 6.3–6.5 R⋆ in the equatorial and also in the polar direction.

6 Conclusions and future work

We studied the hydrodynamic behavior of interaction between the expanding envelope of SNe and aspherical CSM containing an spherically symmetric stelar wind component and an aspherial circumstellar equatorial disk (or disk-like) component that is axially symmetric and resides in the equatorial plane of the star. In the present model we used simplified initial conditions, assuming homogeneous initial distribution of the stellar density and pressure without gravity, which has only a minor impact on the studied process of external interaction.

We discussed three particular 2D models, supplemented by 1D illustrative slices of basic hydrodynamic quantities and of gradients (slopes) of density and temperature, for three different values of stellar pre-explosion mass-loss rate that led to formation of CSM. All the models significantly demonstrate the higher rate of deceleration of SN expansion in case of higher density of CSM and the aspherical evolution of the density, velocity, and temperature structure of the SNe ejecta, where the mass preferably expands to the area outside the dense equatorial disk (cf. McKee & Cowie 1975; McKee et al. 1978; Cowie et al. 1981). We conclude that in case the disk initial base density values ρ_0,csd are of an order of magnitude (or more) higher than the base densities ρ_0,sw of spherically symmetric CSM, such circumstellar equatorial disks effectively slow down the equatorial SN expansion (compared to the SN expansion into regions outside the disk) with significant peaks of density, pressure, and temperature in the contact region of SN-disk interaction. This is according to the presented models and other models with even lower initial CSM densities down to Ṁ_csd ≈ 10⁻⁹ M_⊙ yr⁻¹ with ρ_0,csd of an order of 10⁻⁹ kg m⁻³, where the effects of asphericity vanish, and down to Ṁ_sw ≈ 10⁻⁷ M_⊙ yr⁻¹ with ρ_0,sw of an order of 10⁻¹²–10⁻¹¹ kg m⁻³ below whichthe effects of CSM become negligible. The comparison of power-law gradients of density and temperature (basically described in Figs. 5, 11, and 18), at least for the smooth parts of the expanding area (excluding the regions of steep shocks and discontinuities), indicates their higher values for lower pre-explosion mass-loss rates. The slopes also decrease with time, however, this can be expected from the similarity nature of expansion. We also checked the time limits of similarity solution applicability, where all the models match the time frame with considerable reserve.

The models also indicate the development of Kelvin-Helmholtz instabilities in the zone of shear where the expanding matter flow is distorted along the interface between the spherically symmetric component and circumstellar equatorial disk component of CSM, forming the nearby layers of overdensity.

In our models, we did not consider the effects of radiative cooling, the influence of magnetic fields, or the energy input from the compact remnant of the explosion. These effects are typically mitigated and overlapped by the enormous SNe energy and their impact on the dynamics of the problem is very small (Truelove & McKee 1999), or are connected only with a specific situation or a type of the remnant (like magnetars; see, e.g., Metzger et al. 2017). The reverse shock formed in uniform ejecta is assumed to be radiative for a certain time (Chevalier 1977), but its velocity becomes eventually so high and the ejecta are so tenuous that the cooling timescale for the gas behind the shock exceeds the dynamical time of the SN expansion (Truelove & McKee 1999). It is therefore relevant to use the nonradiative approximation. In this point, Truelove & McKee (1999) showed that the fraction of ejecta mass that cools is only a few percent in the case of n < 5 ejecta, in case of typical values of parameters. Moreover, Chevalier & Fransson (1994) studied the early radiative period for n > 5 ejecta and found it had very little effect upon the dynamics of the SN expansion.

Because of the inclusion of radiative cooling, and therefore a violation of the energy conservation, the shocked shell becomes thinner and denser than in the case of the used self-similar analytical description (Chevalier 1982). Such a shell is however expected to be unstable because of developing (the Rayleigh-Taylor) instabilities (e.g., Chevalier & Blondin 1995), whose effect reduces the rate of conversion of the kinetic energy to radiation (Blinnikov et al. 1998; Moriya et al. 2013a). These instabilities are usually included in 1D models via an approximative “smearing term” (similar in principle to the numerical viscosity). We currently omit such details owing to the complexity of the SN–CSM interaction morphology, but they may be included in the future calculations (cf. the study of Steinberg & Metzger 2018).

Although the detailed description of observational signatures connected with the conclusions of the models is not the subject of this paper, we briefly comment on that issue. From an analytic solution suggested, for example, by McDowell et al. (2018; see also Moriya et al. 2013b or Chatzopoulos et al. 2013), it follows that in cases in which the mass of the disk M_csd is comparable with the mass of the ejecta M_ej, the ratio of the luminosity of the light curve peak is approximately indirectly proportional to the SN diffusion timescale τ_sn; the time when the opacity drops so that the radiation can freely escape. In cases in which the mass of the disk is much smaller (two orders of magnitude or more) than the ejecta mass, then the peak luminosity ratio of the light curve is significantly reduced; the peak luminosities for the same ratio of the SN diffusion times are much closer. Comparing the first two models A and B (where M_csd is comparable with M_ej ≡ M_⋆) with model C, we obtain a ratio of the diffusion timescales of about 3∕2. We may thus expect the light curve peak ratio to be approximately 1.2–1.5 in favor of model C. We may also assume that in case of yet lower disk masses the difference in luminosities would be even less distinct. However, the situation may be complicated according to the line of sight of an observer. We expect that in the case of pole-on observational direction the difference in light curve luminosities would be much smaller (if any) while in case of equator-on (disk-on) observational direction the effects of the disk may correspond to the above estimate. We also expect that because of the asymetry of the outflow the resulting line profiles should depend on the direction of observation showing higher outflow velocities in the polar direction than in the equator-on direction.

As an immediate future step we plan to modify the initial stellar profiles of the density and temperature into a more realistic form, including the internal structure of the progenitor pre-calculated from stellar evolution code. We also compare the current 2D models calculated without gravity with the models with the stellar gravity included. However, in this point we checked the influence of the gravitational force on the expansion profiles using the 1D SNEC calculations, and there is a remarkable impact only in the central region of the original progenitor, which expands self-similarly together with the expansion of the whole envelope. For future 2D models we will implement the equilibrium-diffusion radiation transport solver into the code, taking into account recombination effects in the envelope and the effects of radioactive heating. This will enable us to perform the calculations in radiation-hydrodynamic mode with radiation-matter coupling.

We expect that subsequent study of the SN light curves powered by SN thermal energy excess will lead to a better understanding of the mass and density distribution of the CSM. It may also provide more precise estimates of the disk mass-loss rates and a deeper knowledge of the geometry and mechanism of the mass loss of massive stars.

Fig. 19 Model C: detailed color maps of the density (upper panel) and radial velocity (lower panel) up to the distance 50 R⋆ at time t ≈ 105 h since shock emergence, with the same resolution as in Fig. 8. Contours denote the densities ρ = 10−7, 10−6, 10−5, and 10−4 kg m−3 and the velocities 5000 and 7000 km s−1.

Acknowledgements

We thank Dr. Ondřej Pejcha for helpful comments on details of the studied (or similar) process. We are also gratefulto Milan Prvák for technical support in preparation of the draft of this paper. The access to computing and storage facilities owned by parties and projects contributing to the National Grid Infrastructure MetaCentrum, provided under the program “Projects of Large Infrastructure for Research, Development, and Innovations” (LM2010005) is appreciated. This work was supported by the grant GA ČR 18-05665S and by the grant Primus/SCI/17.

Appendix A Analytical solution of SN expansion

A.1 Basics ofadiabatic similarity solution

Contact discontinuity separates the SN remnant and CSM and propagates among the forward and reverse shock wave. An analytical approach requires a similarity solution (see, e.g., the formalism in Chevalier 1982 or Nadezhin 1985) with the coefficient of proportionality λ(r, t) defined as Kλ(r, t) = rt^−α (where the parameter α is different from α_vis in Sect. 2and K is the scaling constant). Denoting λ_I = 1 the value of λ at the inner (reverse) shock wave, λ_C at the contact discontinuity, and λ_II at the outer forward shock wave (where obviously λ_in(r < r_I) < λ_I < λ_C < λ_II < λ_out(r > r_II)), coupling of the densities in Eq. (17) gives $$\alpha =\frac{n-3}{n-w}.$$(A.1)

Substituting the equations for λ and α into Eq. (17), we introduce the similarity variables U(λ), V (λ), W(λ), and C(λ), which correspond to velocity, density, pressure, and sound speed. The equations for u, ρ, and P are $$\begin{array}{l}u(r\le r_{\text{II}})=K{t}^{\alpha -1}\lambda U(\lambda ),\\ u(r>r_{\text{II}})=0,\\ \rho (r\le r_{\text{C}})=A{K}^{-n}{t}^{-\alpha w}V(\lambda ),\\ \rho (r>r_{\text{C}})=B{K}^{-w}{t}^{-\alpha w}V(\lambda ),\\ P(r\le r_{\text{C}})=A{K}^{2-n}{t}^{\alpha (2-w)-2}{\lambda }^{2}W(\lambda ),\\ P(r>r_{\text{C}})=B{K}^{2-w}{t}^{\alpha (2-w)-2}{\lambda }^{2}W(\lambda ).\end{array}$$

The similarity variable C(λ) is defined as C² = γW(λ)∕V (λ). Following the relation for the coefficient λ, we may express the velocity u_C of the contact discontinuity, $$u_{\text{C}}=\frac{\text{d}{r}_{\text{C}}}{\text{d}t}=\alpha \frac{r_{\text{C}}}{t},U({\lambda }_{\text{C}})=\alpha .$$(A.5)

Substituting Eqs. (A.2)–(A.4) for the domain r ≤ r_II into Eqs. (10)–(13), assuming spherical symmetry, and partially differentiating with respect to t and λ, we obtain the following similarity relations (cf. Nadezhin 1985) for the continuity, momentum, and energy equations, respectively, $$\lambda \left(U-\alpha \right)\frac{V^{\prime }}{V}+3U+\lambda U^{\prime }-\alpha w=0,$$(A.6) $$\lambda \frac{V^{\prime }}{V}+\gamma \lambda U^{\prime }\frac{U-\alpha }{C^2}+2\lambda \frac{C^{\prime }}{C}+\gamma U\frac{U-1}{C^2}+2=0,$$(A.7) $$\lambda \left(U-\alpha \right)\left[(1-\gamma )\frac{V^{\prime }}{V}+2\frac{C^{\prime }}{C}\right]+2(U-1)+(\gamma -1)\alpha w=0.$$(A.8)

The explicit form of Eqs. (A.6)–(A.8) for each particular similarity variable derivative is $$\lambda U^{\prime }=-\frac{\mathcal{S}{C}^2+\left(U-\alpha \right)[3C^2-U(U-1)]}{C^2-{\left(U-\alpha \right)}^2},$$(A.9) $$\lambda \frac{V^{\prime }}{V}=\frac{\eta \mathcal{S}{C}^2+\left(U-\alpha \right)[\left(U-\alpha \right)(3U-\alpha w)-U(U-1)]}{\left(U-\alpha \right)[C^2-{\left(U-\alpha \right)}^2]},$$(A.10) $$\begin{array}{lll}\hfill & \hfill & \lambda \frac{C^{\prime }}{C}=-\frac{\eta \mathcal{S}{C}^2}{2\left(U-\alpha \right)[C^2-{\left(U-\alpha \right)}^2]}\hfill \\ \hfill & \hfill & +\frac{2C^2+(\gamma -1)U(U-1)+\left(U-\alpha \right)[2-(3\gamma -1)U]}{2[C^2-{\left(U-\alpha \right)}^2]},\hfill \end{array}$$(A.11)

where we substitute the constant expressions η and $\mathcal{S}$, $$\eta =\frac{(\gamma -1)\alpha w-2(1-\alpha )}{\alpha (3\gamma -w)-2(1-\alpha )}<1,$$(A.12) $$\mathcal{S}=\frac{\alpha (3\gamma -w)-2(1-\alpha )}{\gamma }>0.$$(A.13)

Multiplying Eq. (A.6) by an auxiliary constant β and subtracting it from Eq. (A.8), the integration gives the implicit adiabaticity integral (see, e.g., Sedov 1959; Nadezhin 1985) $${\lambda }^{2+3\beta }|U-\alpha {|}^{\beta }{V}^{1+\beta -\gamma }{C}^2=Q=\text{\hspace{0.17em}const.},$$(A.14)

where the constant β is $$\beta =\frac{2(\alpha -1)+(\gamma -1)\alpha w}{\alpha (w-3)}.$$(A.15)

The Rankine–Hugoniot relations (e.g., Zel’dovich & Raizer 1967) give the λ-independent similarity variables for the inner strong shock, $$U_{\text{I}}=\frac{2\alpha +\gamma -1}{\gamma +1},V_{\text{I}}=\frac{\gamma +1}{\gamma -1},W_{\text{I}}=\frac{2{(\alpha -1)}^2}{\gamma +1},C_{\text{I}}^2=\gamma \frac{W_{\text{I}}}{V_{\text{I}}}.$$(A.16)

The similarity variables for the outer strong shock are $$U_{\text{II}}=\frac{2\alpha }{\gamma +1},W_{\text{II}}=\frac{2{\alpha }^2}{\gamma +1},C_{\text{II}}^2=\gamma \frac{W_{\text{II}}}{V_{\text{II}}},$$(A.17)

where only the variable V_II has to be evaluated numerically using the principles given in Appendix A.2.

We also express the similarity variables in the zones inside the inner shock wave (λ_in < 1) and outside the outer shock wave (λ_out > λ_II). Combining Eq. (17) with Eqs. (A.2)–(A.4) and assuming u_out, P_in, P_out = 0, we obtain $$\begin{array}{l}{U}_{\text{in}}=1,V_{\text{in}}={\lambda }_{\text{in}}^{-n},C_{\text{in}}=0,\\ U_{\text{out}}=0,V_{\text{out}}=\frac{\gamma -1}{\gamma +1}{\left(\frac{{\lambda }_{\text{II}}}{{\lambda }_{\text{out}}}\right)}^w{V}_{\text{II}},C_{\text{out}}=0.\end{array}$$

From Eq. (A.9) it follows that in very proximity of the contact discontinuity (where U − α → 0) $\lambda U^{\prime }=-\mathcal{S}$. Integrating the left-hand side of this relation from U(λ_C) to a nearby U and the right-hand side from λ_C to a nearby λ, as well as integrating analogously Eqs. (A.10) and (A.11), gives the similarity relations close to r_C in the form $$U-\alpha =-\mathcal{S}\ln \frac{\lambda }{{\lambda }_{\text{C}}},V=\xi {\left|U-\alpha \right|}^{-\eta },C^2=\zeta {\left|U-\alpha \right|}^{\eta },$$(A.20)

where ξ and ζ are arbitrary constants that differ in general on both sides of r_C. However, since u and P are continuous (Zel’dovich & Raizer 1967) through r_C (implying a constant product V C² in Eq. (A.20)), we may write ${(\xi \zeta )}_{\text{i}}={(\xi \zeta )}_{\text{o}}$, where the subscript “i” denotes the zone between the inner shock wave and the contact discontinuity while the subscript “o” denotes the zone between the contact discontinuity and the outer shock wave. To determine the constants ξ and ζ, Eq. (A.20) have to be calculated numerically within the semi-analytical solution (see Appendix A.2).

We obtain the relations for total mass M_I,II and total energy E_I,II confined between the two shocks (Nadezhin 1985) by integration of Eq. (17), with use of the formalism in this section, $$\begin{array}{l}{M}_{\text{I,C}}=\frac{4\pi A}{(n-3)K^{n-3}}{t}^{\alpha (3-w)},M_{\text{C,II}}=\frac{\gamma \text{}-\text{}1}{\gamma \text{}+\text{}1}\frac{4\pi B{K}^{3-w}}{3\text{}-\text{}w}{V}_{\text{II}}{\lambda }_{\text{II}}^3{t}^{\alpha (3-w)},\\ E_{\text{I,II}}=\frac{2\pi A}{(n-5)K^{n-5}}{t}^{(1-\alpha )(n-5)}.\end{array}$$(A.21)

This leads to the important constraint w < 3 and n > 5 for the similarity solution, while Eq. (A.1) implies the condition α < 1.

To evaluate the constants A, B, and K, we identify the mass M_ej of the SN envelope as M_I,C (which tends to be relevant in an advanced time) and the energy E_ej of the envelope as E_I,II. Integrating the first Eq. (A.3) from r_I to infinity, with use of Eq. (A.21), gives the constant A. The constant B results fromthe mass conservation equation for CSM (cf. Eq. (17)) as $B={\dot{M}}_{\text{sw}}/(4\pi R_{\star }^{2-w}{v}_{\text{sw}})$, assuming the stationary pre-explosion CSM velocity v_sw ≪ u(r < r_II). The constant K results from Eq. (A.4) following the continuous pressure at the contact discontinuity r_C. This also implies BK^j−w ≡ A∕K^n−j, where j is an arbitraryinteger.

A.2 Principles of semi-analytical solution of the similarity equations

We basically adopt the principles of the numerical integration of the analytically expressed solution given in Sect. 3 from Nadezhin (1985). We first numerically integrate Eqs. (A.9)–(A.11) with use of constant adiabaticity integral, Eq. (A.14), in the region “i” (where r_I < r < r_C), until the condition 0 < U − α ≪ α is satisfied. The input values of the similarity variables for evaluation of the constant Q = Q_i in Eq. (A.14) are given by Eq. (A.16). We put the similarity variables U, V, C, and λ, calculated in the previous step, to Eqs. (A.20). We thus determine ξ_i, ζ_i, and λ_C to complete the integration between r_I and r_C. Then we integrate the region “o” (where r_C < r < r_II) in a similar way by inserting an input value of λ slightly larger than λ_C, until the condition 0 < U − α ≪ α is satisfied (Eq. (A.20)) again.

We select a trial value of the constant ξ_o to determine the constant ζ_o. The analogous initial conditions for the constant Q = Q_o in Eq. (A.14), noting that Q_o≠Q_i in general, are given by Eq. (A.17). Integration continues until we reach the equality given by the second Eq. (A.4) for λ = λ_II. Because the value of C_II differs in this case from the value given by Eq. (A.17), we use an algorithm that finds such a constant ζ_o, for which Eq. (A.17) is satisfied with the required accuracy while the second equation (A.4) holds as well. Once this ζ_o value is found, the problem is completely resolved because in this case the values of the similarity variables V_II and C_II are also uniquely determined. We study the cases with γ_o = γ_i (see explanation in Sect. 3) and μ_o = μ_i; the possible difference between the mean molecular weight μ of the matter expelled from the star and that of the ambient medium has no effect on the structure of the similarity solution (Nadezhin 1985).

A.3 Constraints of the adiabatic self-similar solution

Regarding the early stages of SN expansion, we may neglect the energy losses due to the volume radiation (Nadezhin 1985). However, the similarity solution is yet limited by the applicability of the adiabatic gas dynamics. The first constraint results from the condition that the length of the mean free path of the particles must be small compared to the characteristic scale Δr of the expanding envelope. In other words, the width ℓ of the shock front must be smaller than the characteristic scale, ℓ < Δr. Following the study of a strong shock wave front in hydrogen plasma by Imshennik (1975; cf. also Nadezhin 1985), we adopt the value (in SI units) $$\mathcal{l}\simeq \frac{8.1\times {10}^{22}}{n_{\text{in}}}{\left(\frac{{\tilde{u}}_{\text{in}}}{{10}^3\text{\hspace{0.17em}km\hspace{0.17em}}{\text{s}}^{-1}}\right)}^4\text{\hspace{0.17em}m},$$(A.22)

where n_in = ρ_in∕m_H is the number density of hydrogen atoms before the shock wave front (m_H is the mass of a hydrogen atom) and ${\tilde{u}}_{\text{in}}$ is the relative flow velocity of the expanding gas before the corresponding shock wave (in the rest frame of the shock wave). Taking the distance r_I,C = |r_I − r_C| between the inner shock wave and the contact discontinuity as the characteristic scale Δr, where we do not consider the outer shock wave because it propagates outside the SN gas, r_II > r_C, following the definition of the coefficient λ gives $$\mathcal{l}<({\lambda }_{\text{C}}-1)K{t}^{\alpha }.$$(A.23)

We evaluate the velocity ${\tilde{u}}_{\text{in}}$ as u_I − u_in following the strong shock Rankine–Hugoniot equations (e.g., Zel’dovich & Raizer 1967). Employing Eqs. (A.2) and (A.16) gives $${\tilde{u}}_{\text{in}}=(1-\alpha )\frac{r_{\text{I}}}{t}.$$(A.24)

The second constraint obviously results from the finite amount of the ejected matter M_ej. The total mass in the SN envelope must obey the inequality $$M_{\text{I,C}}<M_{\text{ej}},$$(A.25)

where M_I,C is given by Eq. (A.21). The adiabatic similarity solution is thus applicable within the time interval t_min < t < t_max, where we obtain t_min from Eq. (A.23) and t_max results from Eq. (A.25).

References

All Tables

All Figures

Fig. 1 Color map of the initial density of the SN progenitor with M⋆ = 45 M⊙, R⋆ = 80 R⊙ (with mean density ${\langle \rho \rangle }_{\star }\approx 0.125{\text{kg m}}^{-3}$) and CSM with assumed Ṁsw = 10−1 M⊙ yr−1 and Ṁcsd = 10−4 M⊙ yr−1 (corresponding to Model A in Sect. 5). The base density of the spherically symmetric part of CSM ρ0,sw ≈ 7 × 10−6 kg m−3, while the equatorial disk base density ρ0,csd ≈ 4 × 10−4 kg m−3 (see also Table 1 in Sect. 5).

In the text

	Fig. 2 Model A: color map of the density structure of the interaction of SN (with progenitor and CSM parameters introduced in Sect. 3, Table 1, and in Fig. 1) with asymmetric CSM that forms a dense equatorial disk, at time t ≈ 16 h since shock emergence. Contours denote the densities ρ = 10−7, 10−6, 10−5, and 10−4 kg m−3.
In the text

Fig. 3 Model A: interaction of SN with asymmetric CSM at time t ≈ 66 h since shock emergence. Upper panel: color map of the density structure. Contours denote the densities ρ = 10−6, 10−5, and 10−4 kg m−3. Middle panel: color map of the velocity structure. Contours denote the radial (expansion) velocities u = 1000, 2000, and 3000 km s−1. Lower panel: color map of the temperature structure. Contours denote the temperatures T = 105, 5 × 105, and 106 K. Characteristic 1D sections of the model in the equatorial and polar plane are shown in Fig. 4.

In the text

	Fig. 4 Model A: one-dimensional slices of variations of hydrodynamical variables in the equatorial plane (x-coordinate, solid line) and in the polar direction (z-coordinate, dashed line), at the same time as in Fig. 3. Panel A: density ρ. Panel B: expansion velocity u. Panel C: pressure P and temperature T. Panel D: specific entropy s.
In the text

Fig. 5 Model A: upper panel: slope parameters n and w of the density (see Eq. (17)) in the equatorial (solid line) and polar (dashed line) direction corresponding to the panel A in Fig. 4, at the same time. The position of the outer shock wave is at r ≈ 11 R⋆ in the equatorial direction, while it is at r ≈ 16.5 R⋆ in the polar direction. There are preserved the initial density slope parameters of CSM, wsw = 2 in the polar direction and wcsd = 3.5 in the equatorial direction, above these radii. The inner envelope density slope parameter increases from n = 0 to n ≈ 9 in the polar direction. The numbers nmax and nmin above and below the graph denote the maximal and minimal slopes within the peaks of the density gradient discontinuities. Lower panel: slope parameters p of the temperature (cf. Eq. (5)) in the equatorial (solid line) and polar (dashed line) direction corresponding to the temperature graph (panel C) in Fig. 4 at the same time. The inner envelope (smooth) temperature slope parameter increases to p ≈ 2.7 at r ≈ 7.8 R⋆ in the equatorial direction while it increases to approximately the same value at r ≈ 11 R⋆ in the polar direction. The slope parameter p outside the outer shock is p = 0.875 in the equatorial and p = 0.5 in the polar direction. The numbers pmax above the graph denote the maximal temperature slopes within the peaks of the T gradient discontinuities.

In the text

Fig. 6 Model A: results of semi-analytical solution of the density ρ, velocity u, pressure P, and temperature T (black line), with the slope parameters n = 7 and wsw = 2, in the polar direction. The method of calculation is described in Sect. A.2. The time is the same as in Figs. 3–5. The velocity at the outer shock wave region (at approx. 16 R⋆) correspondsto the polar expansion velocity in the panel B of Fig. 4 while the other quantities and semi-analytical model input variables are described in Appendix A.2. The blue line depicts the corresponding calculation from our 2D model. The corresponding 1D profile calculated using the SNEC-1.01 code is depicted with a red line (see the description in Sect. 5.1).

In the text

	Fig. 7 Model A: as in Fig. 6, with the slope parameters n = 9 and wsw = 2, in the polar direction (the same method of calculation as in Fig. 6). The semi-analytical profiles are compared to the same polar expansion velocity as in Figs. 4 and 6.
In the text

	Fig. 8 Model A: detailed color maps of the density (upper panel) and the radial velocity (lower panel) structure, up to the distance 50 R⋆ at time t ≈ 165 h since shock emergence. Contours denote the densities ρ = 10−7, 10−6, and 10−5 kg m−3 and the velocities 2500 and 3000 km s−1. The resolution of the simulation is 2400∕480 grid cells in the radial/azimuthal direction within the quadrant.
In the text

Fig. 9 Model B: upper panel: color map of the density structure at time t ≈ 33 h since shock emergence. Contours denote the densities ρ = 10−8, 10−7, 10−6, 10−5, and 10−4 kg m−3. Characteristic 1D sections of the model in the equatorial and polar plane are shown in Fig. 10. Lower panel: color map of the density structure at time t ≈ 66 h since shock emergence. Contours denote the densities ρ = 10−6, 10−5, and 10−4 kg m−3.

In the text

	Fig. 10 Model B: 1D slices of hydrodynamical variables in the equatorial plane (x-coordinate, solid line) and polar plane (z-coordinate, dashed line) at time t ≈ 33 h since shock emergence (as in Fig. 9). Panel A: density ρ. Panel B: expansion velocity u. Panel C: pressure P and temperature T. Panel D: specific entropy s.
In the text

Fig. 11 Model B: upper panel: slope parameters n and w of the density (analogous to Fig. 5) in the equatorial (solid line) and polar (dashed line) direction, corresponding to panel A in Fig. 10 at the same time t ≈ 33 h. The positionof the outer shock wave is at r ≈ 5.5 R⋆ in the equatorial direction, while it is at r ≈ 11 R⋆ in the polar direction. The inner envelope density slope parameter increases from n = 0 to n ≈ 16 in the polar direction. Lower panel: slope parameters p of thetemperature in the equatorial (solid line) and polar (dashed line) direction corresponding to the temperature graph (panel C) in Fig. 10 at the same time. The inner envelope (smooth) temperature slope parameter increases to p ≈ 1.9 at r ≈ 4.1, R⋆ in the equatorial direction while it increases to p ≈ 3.6 at r ≈ 6.6 R⋆ in the polar direction.

In the text

Fig. 12 Model B: detailed color maps of the density (upper panel), radial velocity (middle panel), and temperature (lower panel) structure up to the distance 50 R⋆ at time t ≈ 180 h since shock emergence, with the same numerical resolution as in Fig. 8. Contours denote the densities ρ = 10−7, 10−6, 10−5, and 10−4 kg m−3, the velocities 2500, 3000, 3500, and 4000 km s−1, and the temperatures 7 × 104, 8 × 104, 105, and 5 × 105 K.

In the text

	Fig. 13 Model B: 1D slices of density (left panel) and velocity (right panel)in the equatorial plane (x-coordinate, solid line) and in polar direction (z-coordinate, dashed line) at time t ≈ 180 h since shock emergence (corresponding to Fig. 12). The maximum of the polar velocity, which is already outside the computational domain, is approximately 4400 km s−1.
In the text

	Fig. 14 Model B: 1D slices of temperature (left panel) and specific entropy (right panel) in the equatorial (x-coordinate, solid line) and in polar direction (z-coordinate, dashed line)at the same time as in Fig. 13.
In the text

Fig. 15 Model C: color maps of the hydrodynamical variables at time t ≈ 33 h since shock emergence. Upper panel: color map of the density structure. Contours denote the densities ρ =10−8, 10−7, 10−6, 10−5, and 10−4 kg m−3. Middle panel: color map of the velocity structure. Contours denote the velocities u = 5000 and 7000 km s−1. Lower panel: color map of the temperature structure. Contours denote the temperatures T = 105 and 106 K. Characteristic 1D sections of the model in the equatorial plane and in the polar direction are shown in Fig. 17.

In the text

	Fig. 16 Model C: color map of the density structure of the interaction at time t ≈ 66 h since shock emergence. Contours denote the densities ρ = 10−6, 10−5, and 10−4 kg m−3.
In the text

	Fig. 17 Model C: 1D slices of hydrodynamical variables in the equatorial (x-coordinate, solid line) and in polar direction (z-coordinate, dashed line) at time t ≈ 33 h since shock emergence (as in Fig. 15). Panel A: density ρ. Panel B: expansion velocity u. Panel C: pressure P and temperature T. Panel D: specific entropy s.
In the text

Fig. 18 Model C: upper panel: slope parameters n and w of the density (analogous to Figs. 5 and 11) in the equatorial (solid line) and polar (dashed line) direction, corresponding to panel A in Fig. 17 at the same time t ≈ 33 h. The positionof the outer shock wave is at r ≈ 9 R⋆ in the equatorial direction while it is at r ≈ 18 R⋆ in the polar direction. The inner envelope density slope parameter increases from n = 0 to n ≈ 17 in the polar direction. Lower panel: slope parameters p of thetemperature in the equatorial (solid line) and polar (dashed line) direction corresponding to the temperature graph (panel C) in Fig. 17 at the same time. The inner envelope (smooth) temperature slope parameter increases to p ≈ 3.6 at r ≈ 6.3–6.5 R⋆ in the equatorial and also in the polar direction.

In the text

	Fig. 19 Model C: detailed color maps of the density (upper panel) and radial velocity (lower panel) up to the distance 50 R⋆ at time t ≈ 105 h since shock emergence, with the same resolution as in Fig. 8. Contours denote the densities ρ = 10−7, 10−6, 10−5, and 10−4 kg m−3 and the velocities 5000 and 7000 km s−1.
In the text