© The Authors 2026

1 Introduction

Interactions of supernovae (SNe) with sufficiently dense circumstellar matter (CSM) can significantly alter the evolution of light curves (LCs) and other observable characteristics, such as spectra and polarisation properties; readers can refer to, for example, Chevalier (1984), Chevalier & Fransson (1994, 2006), and Chevalier & Fransson (2003) as well as Fraser (2020) for reviews. In most notable cases, this process gives rise to the super-luminous SNe (SLSNe, e.g. Quimby et al. 2007; Gal-Yam 2012). These can be 1-2 magnitudes brighter at maximum than the brightest SNe powered only by the progenitor’s thermal energy and the radioactive decay of newly synthesised elements; although this may not be the only scenario for producing such extraordinary luminosity (e.g. Woosley 2010; Kasen & Bildsten 2010; Nicholl et al. 2013; Metzger et al. 2015). Another observed feature that reveals SNe interacting with their surroundings is the occurrence of narrow spectral lines, at least in certain phases; these SNe are classified as type IIn (Schlegel 1990).

We performed numerically consistent 2D simulations of axisymmetric radiation hydrodynamics (RHD) of spherical SNe interacting with CSM with two aspherical morphologies (cf. McDowell et al. 2018; Kurfürst et al. 2020, hereafter PK20): a circumstellar disc and a system of bipolar lobes, scaled in individual models with different initial densities. Based on numerically calculated observables, our goal is to identify qualitative differences between the various CSM morphologies, in combination with their density. Once such a correlation is convincingly found, it can contribute to a more precise determination of the mass-loss rate of massive stars before the SN (e.g. Smith 2014). It may also determine the spatial orientation of ejected material during violent eruptions in the pre-SN period or the rate of mass transfer between components in binary systems (cf., e.g. Zasche et al. 2020).

The analysis of the SN ejecta interaction with the surrounding medium is also relevant in the context of denser galactic environments, such as molecular clouds and galactic nuclei. Previous studies of the SN interaction with molecular clouds (Chevalier 1999; Slane et al. 2015) pointed out their X-ray, radio, and γ-ray signatures and the relevance in terms of cosmic ray contribution via diffuse shock acceleration. One of the key signatures of the interaction of the propagating SN shock with the molecular cloud is the compact OH (1720 MHz) maser emission (Frail et al. 1996). The SN interaction with central galactic accretion discs is expected to lead to luminous X-ray/UV/optical flares that could have been mistaken for some tidal disruption events (Grishin et al. 2021; Li et al. 2023; Mondek et al. 2026). In the context of star-jet collisions (Zajaček et al. 2020; Kurfürst et al. 2025), the SN-jet interaction (Vieyro et al. 2019; Bosch-Ramon 2023; Longo et al. 2025) can be the source of high-energy γ-rays and can lead to the production of ultra-high-energy cosmic rays.

Since this work primarily seeks to explore the relationships between the observables during SN expansion and the various morphologies and densities of pre-SN CSMs, we have not created specific models of individual SN events for the purposes of this study. However, we compare the theoretically found observables, especially LC shapes, with some actual observed events and formulate conclusions based on these comparisons.

Not many detailed multidimensional simulations describing the aspherical interactions of SNe with their surroundings, consistently including RHD and non-local thermodynamic equilibrium (NLTE) radiation transfer (RT), have been performed to date (see e.g. Blondin et al. 1996; Vlasis et al. 2016; Orlando et al. 2024, and others). Our papers Kurfürst & Krticka (2019, hereafter PK19), PK20, and Pejcha et al. (2022) provide certain insights into this problem using detailed 2D gas dynamic numerical simulations, while the latter performed calculations of observables semi-analytically, based on the determination of the shock-heating power per unit volume q from the adiabatic gas dynamic simulations, and the subsequent conversion of the volume-integrated $\dot{Q}=\int \dot{q}\text{d}V$ to bolometric luminosity L_SN(t) using Arnett’s law (Arnett 1980, 1982, see Eqs. (8)-(11) in PK20). For this study, we fully employed 2D grey flux diffusion RHD, including radiative cooling using the CASTRO code (e.g. Almgren et al. 2010b, 2011), as well as post-processed 2D Monte Carlo radiation transfer (MC-RT) calculations of observables using the codes SEDONA (Kasen et al. 2006; Almgren et al. 2010a) and SIROCCO (Matthews et al. 2025).

The article is organised as follows: In Sect. 2, we describe the principles used to solve the RHD equations and the numerical solution using the MC-RT codes SEDONA and SIROCCO. Section 3 presents the results of RHD models with stratified density levels. Section 4 demonstrates the calculated LCs, either pseudo-bolometric or in selected colour bands from three viewing angles. Section 5 shows the calculated spectra as a function of the viewing angle, including the temporal evolution of selected lines. Section 6 demonstrates the semi-analytically calculated time evolution of relative polarisation for all models. In Sect. 7, we discuss the results and their context, and propose explanations for possible discrepancies, while Sect. 8 summarises the article and outlines directions for further research in this area. In the Appendix, we describe the technical details and provide figures not included in the main article.

2 Numerical setup

We describe, to a necessary extent, the equations of the RHD involved by CASTRO, as well as the equations of RT entering SEDONA and SIROCCO, including the initial and boundary conditions.

2.1 Description of RHD calculations

We solve the 2D RHD equations using the CASTRO grey radiation solver (Zhang et al. 2011), under the assumption of axial symmetry and the zero azimuthal velocity component of CSM (a negligible effect compared to the SN expansion). The fundamental equations in conservative form employed to solve the evolution of the RHD structure are $$\frac{\mathrm{\partial }\rho }{\mathrm{\partial }t}+\mathrm{\nabla }\cdot \left(\rho \mathbf{u}\right)\mathbf{=}\mathbf{0}\mathbf{,}$$(1) $$\frac{\mathrm{\partial }\left(\rho \mathbf{u}\right)}{\mathrm{\partial }t}+\mathrm{\nabla }\cdot \left(\rho \mathbf{u}\mathbf{u}\right)\mathbf{+}\mathrm{\nabla }\mathbf{p}\mathbf{+}\lambda \mathrm{\nabla }{\mathbf{E}}_{\mathbf{r}}\mathbf{=}\mathbf{0}\mathbf{,}$$(2) $$\frac{\mathrm{\partial }E}{\mathrm{\partial }t}+\mathrm{\nabla }\cdot \left[(\mathbf{E}\mathbf{+}\mathbf{p})\mathbf{u}\right]\mathbf{+}\lambda \mathbf{u}\cdot \mathrm{\nabla }{\mathbf{E}}_{\mathbf{r}}\mathbf{=}\mathbf{-}\mathbf{c}{\kappa }_{\mathbf{P}}(\mathbf{a}{\mathbf{T}}^{\mathbf{4}}\mathbf{-}{\mathbf{E}}_{\mathbf{r}}^{\mathbf{(}\mathbf{0}\mathbf{)}})\mathbf{,}$$(3) $$\frac{\mathrm{\partial }{E}_r}{\mathrm{\partial }t}+\mathrm{\nabla }\cdot ({\mathbf{f}}^{\mathrm{\prime }}\mathbf{u}{\mathbf{E}}_{\mathbf{r}}\mathbf{-}{\chi }^{\mathrm{\prime }}\mathrm{\nabla }{\mathbf{E}}_{\mathbf{r}})\mathbf{-}\lambda \mathbf{u}\cdot \mathrm{\nabla }{\mathbf{E}}_{\mathbf{r}}\mathbf{=}\mathbf{c}{\kappa }_{\mathbf{P}}(\mathbf{a}{\mathbf{T}}^{\mathbf{4}}\mathbf{-}{\mathbf{E}}_{\mathbf{r}}^{\mathbf{(}\mathbf{0}\mathbf{)}})\mathbf{.}$$(4)

Hence, the continuity (1), momentum density (2), and total energy density (3) equations, respectively, with the grey radiation heating and cooling source term on the right-hand side of Eq. (3), are supplemented by the radiation energy density Equation (4). In these equations, ρ, u, p, T, and E are the mass density, velocity, gas pressure, temperature, and total energy density. In radiative terms, E_r is the radiation energy density, λ is the flux limiter (Levermore & Pomraning 1981), c is the speed of light, a = 4σ/c is the radiation constant (σ is the Stefan-Boltzmann constant), and κ_P is the Planck mean opacity. The quantities f′ and χ′ denote the terms f′ = (3 - f)/2 and χ′ = cλ/χ_R, respectively, where f is the Eddington factor and χ_R is the summed absorption and scattering mean Roseland opacity coefficient. The values of these coefficients and the relevance of their implementation are described further in Sect. 2.3. The terms with superscript (0) ($E_r^{(0)}$) denote the corresponding quantity in the co-moving frame, while the terms without a superscript denote the quantity in the lab frame (see Almgren et al. 2010b; Zhang et al. 2011, for details). The set of equations is closed with the ideal gas equation of state (EOS) e = p/(γ - 1) (gamma-law EOS in CASTRO), where e is the gas internal energy density, along with the monatomic gas ratio of specific heats γ = 5/3.

2.2 Initial state of SN, spherical, and aspherical CSM components

Similar to PK20, we first use the 1D SNEC code (Morozova et al. 2015) to propagate the SN shock wave through the progenitor, a red supergiant (RSG) star with a mass M_* ≈ 12.3 M_⊙ and a radius R_* ≈ 1000 R_⊙ at the moment of core collapse. The explosion is initiated as a thermonuclear bomb with energy 10⁵¹ erg, with M = 0.05 M_⊙ of radioactive nickel, and Paczynski EOS (radiation included). Regarding the higher masses of the CSM in the denser models, unlike PK20, we now include the gravitational force using the ’Poisson’s gravity’ implicit solver in CASTRO; all other parameters are set to their default values. We do not consider the star’s rotation due to its negligible effect on the investigated process. At the moment of shock breakout, we remap the following quantities from SNEC to CASTRO using its 1D spherically symmetric interpolation routine: density, velocity, temperature, pressure, and the abundances of H, He, C, O, Si, Fe, and Ni. However, in the hydrodynamic evolution using CASTRO, we no longer take radioactive heating into account.

We insert two components of the gaseous surroundings of an SN as an initial state: the spherically symmetric stellar wind and an aspherical component. Two types of aspherical CSM are implemented separately within five different models: two models of a circumstellar disc, labelled as DMIN and DMAX (hereafter DM), differing by two orders of magnitude in initial density values; more parameters and initial scaling of the DM models are described in Appendix B. Three models of a bipolar lobe system are stratified by different initial density levels, formally labelled as models LMIN, LMOD, and LMAX (hereafter LM). We list the basic parameters of all models in Table 1: The first column containing numbers, Ṁ, shows the progenitor s mass loss rate that produces the aspherical formation in the models. The second column M_w+d quantifies the total mass of the surrounding CSM prior to the SN explosion; the third column ρ_0,w lists the base densities of the spherically symmetric stellar CSM (wind), while the last column ρ_{0, aspherical CSM} shows the base densities in the circumstellar disc midplane in models DMIN and DMAX only because, in the case of bipolar lobes, such a base density cannot be identified.

Both types of CSM models (disc vs. bipolar lobes) thus represent relatively large transfers of matter or outbursts in the pre-SN epoch. We do not simulate an exact, observationally motivated model of the CSM, such as the Homunculus Nebula, in this work; we use its iconic structure only as a representation of a non-trivial morphology of lobes, forming transverse density barriers to the SN expansion and forcing SN-CSM interaction to be even more asymmetric (see also other similar models in Kurfürst & Krticka 2023; Kurfürst 2024; Kurfürst & Holoubek 2024). In the RHD model, as well as in RT calculations, we do not consider dust grains, which can be relevant for the evolution of complex plasma (Czerny et al. 2023b) and for RT calculations, especially at longer infrared wavelengths. We will include dust in the future models of the SN-CSM interaction.

2.3 Boundary conditions, RHD code grid, and its setup

The main RHD code parameters are as follows: we set the 2D cylindrical radial-vertical (ϖ-z) grid, where ϖ ∈ 〈0, R_max〉, z ∈ 〈-R_max, R_max〉, and R_max = 3.475 × 10¹⁶ cm ≈ 500R_* are defined with 4200 × 8400 grid cells in the two specified directions, respectively; the numerical reconstruction scheme is the standard piecewise parabolic method (PPM). The boundary conditions (BC) for the gas-dynamic solution are as follows: symmetry BC at the left boundary (at ϖ = 0, which represents the axis of the cylindrical system), while all the other three boundaries are outflow. The ideal-gas EOS, implementation of gravitation, and chemical elements are mentioned in Sects. 2.1 and 2.2; we do not include a nuclear reaction network.

Within the RHD solution, we employed the grey radiation solver stabilised by the flux limiter of the Levermore & Pomraning (1981) type, with the Dirichlet left BC, where E_r = aT⁴ at ϖ = 0, and Neumann BCs with the radiation flux F_r = - (cλ/κ_R) ∇E_r at the three other boundaries. We parameterise, for simplicity, the initial constant Rosseland mean opacity κ_R = 0.34 cm⁻¹ (fit for electron scattering with the typical hydrogen abundance X = 0.7), while the Planck mean opacity is not relevant here since the total opacity in the diffusion coefficient for the grey solver is set in CASTRO by the Rosseland mean opacity κR, and the Planck mean opacity is employed when the multigroup radiation solver is used (Almgren et al. 2011). Incidentally, the value κ_R = 0.334 from (web Los Alamos National Library) corresponds to our situation at the moment of the SN shock breakout, when the values at the shock front from the SNEC calculation are ρ ≈ 3.12 × 10⁻⁹ g cm⁻³, T ≈ 2.9 × 10⁵ K, and the relative abundances of key elements are 0.67 for H, 0.313 for He, 0.0056 for C, 0.0075 for O, etc. We list additional parameters in Appendix B.

2.4 Description of post-processed MC-RT calculations

We use the data calculated in the RHD code as input for the MC-RT code SEDONA (Kasen et al. 2006; Almgren et al. 2010a), which is primarily designed for calculations of either LTE or NLTE radiation transfer in SNe. We do not use the simplified spherically symmetric hydrodynamic procedure in SEDONA with homologous expansion v ∝ r; instead, we use individual output data files from the RHD code (a total of 250 time steps across all models) as a sequence of stationary MC-RT pre-calculations from which the SEDONA model file reads the geometrical parameters, density, temperature, radiative energy, velocity components, and the abundance of chemical elements involved (listed in Sect. 2.2). The final emergent spectra and light curves of the SN model are obtained by collecting all escaping photon packets; they are binned in time of arrival, observed wavelength, and escape direction (i.e. viewing angle). A large number of bins in each dimension must be used to achieve the requisite resolution, along with enough packets collected in each bin to provide adequate photon statistics (Kasen et al. 2006), outputting the L_ν(t_n,μ) quantity, i.e. the monochromatic, directionally ’observed’ ’luminosity’. We use the frequency grid in the range ν_min = 10¹⁴ Hz, ν_max = 2 × 10¹⁶ Hz for the LC calculations, corresponding to 150-30 000 Å, logarithmically scaled with an increment of q = 0.001 (where dν = qν), resulting in 5300 frequency grid points, with 5 × 10⁷ photon packets (‘particles’) to initialise the simulations. We set 25 angular directions μ = cos θ, where θ is the standard spherical polar angle, to cover a large range of observer’s directions, from (nearly) polar to equatorial. We employ the NLTE regime for collisions and for calculations of the following opacities: electron scattering, bound-bound, free-free, and bound-free; for the elements listed in Sect. 2.2, in variant configurations demonstrated in Fig. A.4. Since we include the gas temperature from the RHD CASTRO calculations, we do not update the gas temperature with opacity iterations; an additional test showed that whether we include this temperature update or not makes a negligible difference; moreover, a temperature update that includes only radiation would be inconsistent. We do not include the artificial SN core parameters. The bolometric (pseudo-bolometric, within the given frequency range) LCs (see Figs. 3 and 5) and AB magnitudes (Oke 1965, see Figs. 4 and 6) for each model time step n are calculated as $$L_{\text{bol}}(t_n,\mu )={\int }_{{\nu }_{\text{min}}}^{{\nu }_{\text{max}}}{L}_{\nu }(t_n,\mu )\text{d}\nu .$$(5)

The AB magnitude for a frequency band [b] and an angular direction μ is $$M_{\text{AB}}(t_n,\mu ,\text{b})=-2.5{\log }_{10}\left[\frac{\int T_{\text{b}}(\nu ){\nu }^{-1}{f}_{\nu }(t_n,\mu )\text{d}\nu }{\int T_{\text{b}}(\nu ){\nu }^{-1}\text{d}\nu }\right]-48.60,$$(6)

where T_b(ν) is the frequency-dependent transmission function for a given filter band [b]; the term ν⁻¹ (in fact (hν)⁻¹, where the Planck constant h cancels) indicates a photon energy-counting detector (CCD) response (e.g. Blanton & Roweis 2007), and the monochromatic spectral flux density f_ν = L_ν/(4πd²), where d ≡ 10 pc. We note that within the numerically processed flux-filter convolution in Eq. (6), the tabulated frequency values of T_b(ν) for all selected bands must first be interpolated to the frequency grid of f_ν using the spline interpolation function.

We calculate the broadband spectra (see Fig. A.4) in a shorter frequency range, covering the near-infrared, optical, and UV regions, roughly up to the Lyman jump, namely from 3.3 × 10¹⁴ Hz to 3.3 × 10¹⁵ Hz, with 2305 frequency grid points and an increment of q = 0.001 on a logarithmic scale. Unlike the LC calculation, we now involve the ’last-iteration pump’ function, where, in the last computational run, there are 10 times more photons initialised to achieve a higher final resolution of the spectral features. We also demonstrate the time evolution of different line-of-sight profiles of a hydrogen spectral line, represented by Ly-α, in greater detail (see Fig. A.5); we calculate the individual line profiles using the Sobolev-based MC-RT code SIROCCO for 2.5-dimensional axisymmetric astrophysical outflows (Matthews et al. 2025; Long et al. 2025) in the ‘classic’ mode (see Sect. 3.2 in Matthews et al. 2025), where, unlike in the ‘hybrid macro-atom’ mode, the wind itself may become a radiation source. We set the range 1150-1280 Å with 10 000 frequency-grid points and 10⁷ photons per computational cycle. We inserted the mandatory central object with a mass 1.3 M_⊙, radius 7 × 10¹³ cm, and a surface temperature 10⁴ K (roughly mimicking a neutron star remnant mass, surrounded by a relatively small optically thick environment), while the analytically imported parameters of a disc were excluded. The reason for choosing these parameters - particularly the temperature of 10⁴ K - in combination with the other parameters is that they approximate (albeit in a somewhat idealised manner) the actual conditions at a given distance from the centre of the simulation, and at the same time, best generate the Ly-α spectral line in the models at all the specified times. We note that the SIROCCO code requires the mandatory selection of a ’central’ object with a given mass, size, and temperature; however, its selected radius is so small (relative to the total extent of the computational domain) that it does not interfere significantly with the possible nebulosity of the entire configuration. For a similar calculation of the Hα spectral line, a larger radius of 1.1 × 10¹⁴ cm was chosen in order to achieve a somewhat greater contribution to the total radiative flux, which proved to be essential for the clear generation of the Hα line.

We entered a one-component radiative wind of the imported type (whose parameters of coordinates, velocity components, density, and temperature are loaded from a file) in a cylindrical system. The calculation of each spectral line was done in 25 ionisation cycles and 100 spectrum cycles (whose flowcharts are shown in Fig. 2 of Matthews et al. 2025). The atomic data were imported from the SIROCCO ‘kurucz91.ls’ library (Long et al. 2025), containing ≃10⁴ line transitions selected from the Kurucz CD-ROM 23 line-list with ~500 000 lines (Kurucz & Bell 1995).

3 RHD models of interaction

We have calculated the SN interactions with two morphologies of aspherical CSM: circumstellar disc and bipolar lobes (see PK20), now with various levels of their pre-explosion density. We calculated our new RHD models using the CASTRO code (Almgren et al. 2010b, 2011). In the main article, we show images of the RHD models with the lowest mass of aspherical CSM (DMIN and LMIN), as we consider them the most realistic; models with higher masses provide an interesting comparison and insight into the manifestations of SN-CSM interactions; however, given the observations, they are apparently mass-overestimated.

3.1 SN interaction with circumstellar disc - DM models

We have previously studied the interaction of SN ejecta with the circumstellar disc in PK19; PK20 (also e.g. Vlasis et al. 2016; McDowell et al. 2018, and a few others), but only in an adiabatic hydrodynamic approximation, without including radiation hydrodynamics, which accounts for radiative cooling and heating.

Analogous to previous studies, we see that when the SN ejecta interact with the disc, the local shock wave propagates to some extent through the disc itself, but mainly into more ‘polar’ regions outside the disc, generating an over-density shoulder near the vertical edge of the disc. In the DMIN model (see Fig. 1), the ratio of the disc mass to the ejecta mass is relatively small (comparable to model A in PK20), so the SN ejecta propagate relatively unimpeded into the disc region in the equatorial direction. The DMAX model (see Fig. A.1), with disc density parameters two orders of magnitude higher than in the DMIN model, shows (as expected) much greater deceleration of expansion in the disc direction, while the expansion of material in the polar direction is roughly similar to that in the DMIN model. However, there is a significantly greater over-density of a spiralarm-like clump of accumulated material within the disc opening angle, which is not redirected in the polar direction.

However, an important difference from the previous study appears to be the lower expansion velocity of the ejecta; for example, the shock region in the polar direction reaches an outer edge of the computational domain R_max = 450R_* (≈3.12 × 10¹⁶ cm) in approximately 400 days in the adiabatic solution in PK20, while it reaches a slightly greater distance ≈3.5 × 10¹⁶ cm in approximately 700 days within the RHD solution in the current study. Other hydrodynamic phenomena, such as the development of Kelvin-Helmholtz instabilities near the interface between freely expanding ejecta in the polar direction and the ’edge’ of the disc, remain essentially similar.

The left panel of Fig. A.8 shows the density map of helium and iron in two more evolved times in model DMIN. We subsequently input the abundances of specific chemical elements (listed in Sect. 2.2) into the code SEDONA for the calculation of broadband spectra, illustrating the intensity and time evolution of the corresponding spectral lines.

3.2 SN interaction with bipolar lobes - LM models

We also study the hydrodynamic evolution of SN ejecta colliding with nebulae consisting of bipolar lobes, resembling, for example, the structure of the well-known Homunculus Nebula surrounding the binary star η Car. However, its size is downscaled, and its total mass (and thus the density structure) is scaled in three levels, from lowest to highest: models LMIN, LMOD, and LMAX. We can see that in the LMIN model, the SN ejecta is only partially deformed according to the shape of the original CSM; for example, in the equatorial direction, the ejecta can expand almost without slowing down, and its shock wave front in this direction is outside the computational domain at the 730-day ultimate simulation time. On the other hand, the densest lobe model LMAX forces the ejecta to expand almost exclusively according to the shape of the original lobe structure, and the expansion velocity is significantly lower in this case. Furthermore, we can see that in the case of the thinnest LMIN model, the expansion velocity (in the initial and middle periods of the simulation) is highest in the equatorial direction; in the LMAX model, it is highest in the polar direction; and in the intermediate LMOD model, it is highest in the intermediate direction of approximately 45 °. Other hydrodynamic phenomena, except for the lower average expansion velocity (cf. Sect. 3.1), correspond to those in the PK20 study, including reflecting shock waves and the fact that, due to the complex geometry of the CSM, it is difficult to identify various hydrodynamic instabilities, as is the case in the simpler geometry of the circumstellar disc. Similar to the left panel of Fig. A.8, and alongside other contexts for calculating observables, the right panel of Fig. A.8 illustrates the density map of helium and iron at two later times in model LMIN.

Fig. 1 Model DMIN. Selected evolutionary stages of SN interaction with a ’thinner’ circumstellar disc calculated by the RHD code. The four columns show the snapshots at times t = 0, 240, 480, and 730 days. The individual rows shows four different quantities: from top to bottom it is density ρ, velocity magnitude v, temperature T, and radiative energy density Er. The resolution of the bitmaps had to be significantly reduced compared to resolution of the simulations in order to reduce the file size. The animated version of this image is available online.

4 Light curves and AB magnitudes

Unlike in PK20, where we calculated the LCs using the semi-analytical approach (Arnett’s law, Arnett 1980, 1982), taking into account the semi-analytically calculated shock heating rate Q̇, we used the 3D MC-RT code SEDONA (Kasen et al. 2006; Almgren et al. 2010a) in this study to calculate robust observables (LCs, AB magnitudes, broadband spectra), while also having the ability to compare their appearance from different viewing directions.

Fig. 2 Model LMIN. Selected evolutionary stages of the SN interaction with the ’thinnest’ bipolar lobes structure calculated by the RHD code. The four columns show the snapshots at times t = 0, 240, 480, and 730 days. The individual rows show four different quantities: from top to bottom it is density ρ, velocity magnitude v, temperature T, and radiative energy density Er. The animated version of this image is available online.

4.1 SN interaction with circumstellar disc - DM models

Bolometric (pseudo-bolometric) as well as visual LCs calculated for models DM (see Fig. 3) show a relatively smooth evolution over time without significant drops, re-brightening, or bumps related to irregularities and asymmetries in the surrounding CSM, as well as the resulting changes in the energy balance, especially during the passage of shock waves in the SN ejecta as they pass through denser or rarer environments. This corresponds to the smooth nature of the surrounding CSM configuration and, in principle, also to the results of PK20. In model DMIN, we see a significant difference in the polar and equatorial views only during the early phase, up to approximately 100 days. Equatorially ‘measured’ LCs are even brighter at this time, which is probably due to the larger viewing angle at which the intense shock region with a higher temperature is unobscured by the optically thick region (cf. e.g. the third row’s second panel in Fig. 1, even though this is a later time). This occurs regardless of the potentially higher shielding by the disc, as this shock area with high temperature is effectively shielded even in the polar direction to a greater distance by an expanded and, therefore, thicker layer of material ejected from the SN. Later, the two directional views show the same profiles due to the relatively rapid drop in the luminosity of shorter wavelengths, including UV, in this model (cf. Fig. 4), and the subsequent reprocessing of radiation into visible light, IR, and even longer wavelengths, where the re-emitted luminosity is practically the same in different directions. It is also evident that for up to approximately 50 days, the luminosity of this model is greater in shorter wavelengths, including UV, than in longer wavelengths.

For comparison, we have also included in Fig. 3 an analytically calculated LC (similar to that in PK20, where we used Arnett’s law, Eq. (10)) of an SN powered by the decay of 0.28 M_⊙ of radioactive nickel (we remind that radioactive heating is not included in the current numerical models of SN-CSM interactions). In this context, it should be noted that the interaction energy that subsequently powers the LC is generated only in a certain section of the full spatial angle; therefore, the LC consists of two components: one generated by the interaction with dense, aspherical CSM and the other produced by noninteracting SN ejecta in the region outside the dense CSM. If, for example, the denser CSM were located further away from the progenitor star, there would be a double maximum in the LC: the first generated by the thermal energy of the SN explosion itself, plus the chain of radioactive decay. Only after the ejecta collide with the denser CSM would a secondary increase in luminosity caused by the SN-CSM interaction appear, as illustrated, for example, by the B series models in PK20. In the theoretical models presented in this work, however, both of these components are combined into one because dense, aspherical CSM is located directly next to the progenitor star.

The LCs (bolometric as well as visual) are expected to be brighter in the DMAX model due to the greater interaction energy. The DMAX model also shows a larger difference between the LCs ‘observed’ in the polar and equatorial directions, even in much later times. The reason for the higher luminosity in the equatorial direction is likely the same as that in the DMIN model. The reason for the much longer-lasting difference between the poles and the equator in terms of apparent luminosity is the much slower decrease in brightness at shorter wavelengths, where the difference is greater, as well as the increased brightness during the middle period, around 400-500 days, in longer wavelengths, apparently due to a greater redistribution of radiative energy in denser CSM.

Fig. 3 Pseudo-bolometric, directionally ‘observed’ LCs of the two models of the SN interaction with circumstellar disc-like structures, denoted as model DMIN (left panel) and model DMAX (right panel), up to approximately 2 years (see Sects. 3.1 and 4.1). The upper pair of lines shows the LCs calculated within the ’largest’ frequency range 1014-2 × 1016 Hz (see Sect. 2.4), while the bottom pair of lines shows the same within the range (3.875-7.825) × 1014 Hz, which roughly corresponds to the frequency range of visible light (for LCs corresponding to AB magnitudes in main filters, see Fig. 4). We demonstrate here the LC profiles for two different viewing angles in the graphs, the pole-on (polar) view marked in blue and the equatorial view marked in dashed black line for the bolometric LC (labelled Lbol in the legend), and the polar view marked in red dash-dotted line and the equatorial view marked in dotted black line for the visual LC (labelled jvis in the legend). The black solid line labelled ’no CSM’ corresponds to an SN without CSM interaction, heated with 0.28 M⊙ of radioactive Ni which roughly corresponds to the brightest normal H-rich SNe (Pejcha & Prieto 2015; Müller et al. 2017) or the population mean of stripped SNe (Anderson 2019).

Fig. 4 Calculated AB magnitudes for the two models of the SN interaction with a circumstellar disk-like morphology. Each row in the figure corresponds to one of the models DMIN or DMAX labelled on the right. The three columns demonstrate the three various viewing angles marked at the top. The AB magnitudes are calculated for nine filters listed in the left top and bottom panels.

Fig. 5 Pseudo-bolometric, directionally ‘observed’ LCs of the three models of the SN interaction with a bipolar lobes structure, denoted as model LMIN (left panel), model LMOD (middle panel), and model LMAX (right panel), up to approximately 2 years (see Sects. 3.2 and 4.2). The meaning and structure of the image, the order and types of lines, and the legend are the same as in Fig. 3.

4.2 SN interaction with bipolar lobes - LM models

The bolometric and visual LCs in models LM show a significant two-peaked profile with an intermediate drop at around 100 days. This ‘bumpy’ evolution corresponds to the collisions of an SN shock expansion with relatively dense lobes positioned transversely to the direction of SN expansion propagation. The increase in luminosity as a result of the interaction is relatively steep from the outset due to the combination of the shock wave passing through the dense CSM in the polar direction and the inner lobes in the equatorial direction. The second luminosity peak is caused by the shock wave passing through the outer lobe in the equatorial direction. The later interaction of the SN with the lobes in the polar direction is manifested only by a slight undulation in the shape of the LC around 400-500 days.

All other aspects are, in principle, similar to those in the case of interaction with the disc in Sect. 4.1. The model LMIN shows a sharp early maximum, particularly in the equator-ward luminosity (left panel in Fig. 5), caused by a dominant luminosity in short (UV) wavelengths, which is clearly visible from the top row of Fig. 6. The model LMOD shows the most prominent yet relatively smooth drops in luminosity between SN-lobe interactions (see the middle panel in Fig. 5 and the middle row in Fig. 6). On the other hand, the model LMAX exhibits remarkable wave-like disturbances in shorter wavelengths and in the medium simulation time around 400 days, most likely caused by numerical instabilities in the MC calculation of radiative transfer in a relatively optically thick environment (see the right panel in Fig. 5 and the bottom row in Fig. 6).

The AB magnitude graphs in Fig. 6 clearly show much more pronounced peaks, rises, and falls in brightness at shorter wavelengths, especially in the UV region, which can thus perfectly trace any irregularities in the surrounding CSM and their distribution at these wavelengths. The detailed calculations of the AB and other types of magnitudes represent an important theoretical prerequisite for comparing luminosity, for example, in the UV range, particularly in connection with the planned UV missions, such as the wide-field survey mission ULTRASAT (near-UV photometry; Ben-Ami et al. 2022; Shvartzvald et al. 2024), the transient follow-up mission QUVIK (Quick Ultra-Violet Kilonovae surveyor; Werner et al. 2022) for near-UV and far-UV high-cadence photometry, which is scientifically coordinated by the Masaryk University stellar and high-energy astrophysics teams (Werner et al. 2024; Krtička et al. 2024; Zajaček et al. 2024, 2025), and the planned UV photometry and spectroscopy mission - Ultraviolet Explorer (UVEX; Kulkarni et al. 2021).

5 Spectral patterns

5.1 Broadband spectral evolution

We calculate the time evolution of the broadband spectral features for the five specific models using the SEDONA code (Kasen et al. 2006, see also the technical description of the code in Sect. 2.4 and the connection between the output from the CASTRO code and the SEDONA input therein). Figure A.4 demonstrates the broadband spectra in 3 or 4 selected instant snapshots labelled in the pictures. Each spectral profile in each selected time is also represented by a split graph: coloured and black, showing either polar (coloured) or equatorial (black) viewing direction.

For a better overview of the presence or absence of certain spectral lines, we show in Fig. A.4 these spectra divided into three columns, including a gradually increasing number of specific chemical elements: only H and He in the left column, medium-mass elements C, O, and Si in the middle column, and Fe and Ni in the right column. In the overall spectrum profile, especially in the case of heavier CSMs, jumps in the continuum are clearly visible, particularly the Balmer jump near 3650 Å (the continuum curve is partially affected by the fact that the graphs are scaled in log-log form).

In this sense, the display of broadband spectra does not aim to show the detailed shape of spectral lines, as the current version of the SEDONA code does not include all the physics necessary for accurate nebular spectrum calculations (D. Kasen, private communication, October 10, 2025). For this reason, the SEDONA code apparently displays only the absorption curve of spectral lines without the corresponding emission.

Fig. 6 Calculated AB magnitudes for the three models of the SN interaction with a bipolar lobes-like morphology. Each row in the figure corresponds to one of the models LM, labelled on the right. The other features of the figure are analogous to Fig. 4.

5.2 Individual spectral-line evolution

We also attach a brief study of the temporal evolution of a typical hydrogen spectral line, represented here by the Ly-α line. The technical details of the calculation process using the code SIROCCO (e.g. Matthews et al. 2025, also the private communication with S. Sim and Ch. Knigge, November 03 and December 10, 2025) and the justification for the choice of the Ly-α line are described in Sect. 2.4. Figure A.5 shows the evolution of this line in all five models, arranged in individual rows for three different viewing directions, represented by the three columns.

The Ly-α emission in the models is very characteristic, as it predominantly exhibits an asymmetric profile with multiple peaks. During the relatively early phase of up to approximately 100 days, Ly-α develops in the SN-disc (DM) and SN-lobes (LMIN and LMOD) models from the polar view as a somewhat broad emission plateau, slightly shifted to the blue and concentrated in the range of approximately -4000 to 4000 kms⁻¹. In the intermediate view of 45°, such a plateau is clearly narrower for the given period, in the range of approximately -2500 to 2500 km s⁻¹, and is present only for SN-disc interactions, while for SN-lobes interactions, a more spiky profile prevails here. In the case of the equatorial view, this flat profile is very narrow to indistinguishable; especially for the SN-disc interaction model, it has a strikingly narrow central absorption corresponding to a dense, static environment in the direction of view. In the following period, after approximately 100 days, this ‘flattish’ spectral line profile disappears and is replaced by a more or less doublepeaked (horn-like) profile, which is very pronounced in the case of the SN-disc interaction (cf. e.g. the theoretical spectral line profiles in Jerkstrand 2017). At the same time, it gradually transitions from slightly blue-shifted asymmetry to a red-shifted peak; we attribute this to the increasing rate of the resonant scattering of Ly-α photons in neutral hydrogen moving outward in the SN ejecta or in the surrounding CSM (further discussion of this effect can be found in Sect. 7).

Figure A.6 demonstrates the typical, expected time-evolution of the correspondingly calculated Hα spectral line. In the early (photospheric) phase, due to the relatively dense and extensive surrounding layers of the CSM lasting up to approximately 100 days, the Hα emission profile exhibits only a very narrow emission peak, typical of Type IIn SNe, caused by the ionisation of the outer layers by the SN’s short-wavelength radiation. Later, after the transition from the photospheric to the nebular phase, the line broadens and forms a double-peaked profile due to emission from exposed, deeper, faster-moving layers, while the narrow emission peak produced by the ionised, static, or much slower-moving surrounding medium remains (cf. e.g. Yan et al. 2015; Inserra et al. 2018, etc.). In addition, the Hα profile, like the profiles of the other Balmer lines, is also predicted by these models to be much weaker (in terms of its equivalent width) than the Ly-α. This is also indicated by Fig. A.4, even though the lines are shown in absorption there.

The manifested time evolution of calculated spectral line details corresponds well to the measured profiles during an SN expansion. For example, Smith et al. (2015) demonstrates the case of measured Hα spectral line evolution in the PTF11iqb SN, up to approximately 1100 days, with very similar evolutionary features in principle (see fig. 5 therein). The semi-analytical models A and C in PK20 also indicate the characteristic similarity in the basic properties of the evolution of a typical spectral line for the corresponding SN-CSM interaction morphologies.

6 Polarisation signatures

Similar to our PK20 study, we again estimate the degree of relative polarisation P_R of the axially symmetric SN ejecta using the analytical formulas of Brown & McLean (1977) and Brown et al. (1978), including Thomson scattering in an optically thin environment irradiated by an optically thick internal source. Due to the complexity of calculating and estimating the position of the spatially fragmented photosphere within the expanding material, we do not calculate the polarisation numerically at this stage. The semi-analytical approximations mentioned above provide some insight into the morphologies of axisymmetric CSM and allow for a comparison of polarisation signatures depending on the density and/or mass of CSM in individual models.

The relative polarisation of Thomson scattering P_R (%) is given by the relation (Brown & McLean 1977, Eq. (23)) $$P_R\simeq ⟨\tau ⟩(1-3\gamma ){\sin }^2\theta ,$$(7)

where the spherical polar angle θ = 0 indicates the polar viewing direction, and θ = π/2 indicates the equatorial direction. We therefore provide the highest values according to the direction of observation (see Table 2 and Fig. A.7), ’seen’ from the equatorial direction, while the values for other polar angles can be easily calculated using Eq. (7). The shape factor γ and the averaged Thomson scattering optical depth 〈τ〉 of the envelope (Brown et al. 1978) are $$\gamma =\frac{{\int }_{r_{\text{min}}}^{\mathrm{\infty }}{\int }_{-1}^{1}n{\mu }^2\text{d}r\text{d}\mu }{{\int }_{r_{\text{min}}}^{\mathrm{\infty }}{\int }_{-1}^{1}n\text{d}r\text{d}\mu },⟨\tau ⟩=\frac{3}{\text{16}}{\sigma }_T{\int }_{r_{\text{min}}}^{\mathrm{\infty }}{\int }_{-1}^{1}n\text{d}r\text{d}\mu ,$$(8)

where σ_T is the Thomson scattering cross-section, n(r,μ) is the electron number density, and μ = cos θ. The lower limit r_min in radial integrals means either a fixed radius of 10¹⁶ cm or a limit of τ ≤ 1 (for the two versions, see further in this paragraph).

Assuming a completely ionised pure hydrogen envelope for simplicity, the electron number density is n(r,μ) = ρ(r,μ)/m_H, where m_H is the mass of a hydrogen atom. Equation (8) implies γ = 1/3 for a spherically symmetric morphology, while γ < 1/3 (P_R > 0) is for an oblate and γ > 1/3 (P_R < 0) for a prolate mass distribution.

Figure A.7 shows the relative polarisation P_R calculation graph, including the optically thin region τ ≤ 1 for the polar angle θ = π/2 viewing direction, calculated using the short characteristics solver (Kurfürst et al. 2018, PK20). Since the formula (7) applies exclusively to optically thin media where τ ≪ 1, the estimate of the relative polarisation may be only marginally relevant; whereas including layers with even greater optical depth could yield non-physical results. However, since we consider this solution to be a starting point serving as a preliminary semi-analytical estimate for further, more advanced calculations, we include it here. The values of γ, 〈τ〉, and P_R for four selected times t ≈ 100, 300, 500, and 730 days are listed in Table 2. All models also show a prolate distribution of the ejecta with a relative polarisation P_R < -1.5% at later times (cf. Fig. 13 in PK20).

7 Discussion

7.1 Relevance of the current models

The models shown here represent a systematic upgrade of some of our previous SN-CSM interaction models in PK20 in the following sense: we systematically performed more sophisticated simulations of RHD models, including a proper treatment of radiation in light-curves, broad-band spectral temporal evolution, and the time evolution of a particular angle-viewed spectral line profile to provide more realistic predictions for observations. The presented study does not currently aim to demonstrate exact compliance with specific SN-CSM interaction events; rather, it represents a ‘work iteration’ of such a process in which it is challenging to find a configuration of the surrounding CSM in combination with the corresponding parametrisation of the progenitor star that would lead to observable characteristics fully consistent with some of the specific observed phenomena (e.g. Aretxaga et al. 1999; Stritzinger et al. 2012; Arcavi et al. 2017; Nyholm et al. 2017; Kangas et al. 2022; Kumar et al. 2025; Moriya 2026, among others). At the same time, our models were primarily aimed at the technically sophisticated handling of all numerically conceived computational processes, including RHD modelling, MC-RT calculations of LCs in individual relevant bands, and spectra—both broadband and of selected spectral lines—all in a multidimensional form (currently as axisymmetric 2.5D in SIROCCO code), allowing their study depending on the observer’s line of sight. This will allow us to estimate the situation for the next stage of the study, which will involve a better setting of such parameters and CSM configuration in order to achieve a more realistic match with the observed characteristics; for example, the initial decrease in the brightness of the LCs of interacting SNe, which is usually observed before irregular increases or re-brightenings occur, suggests that the dense CSM is probably located at a greater distance from the progenitor star. Only the calculation of relative polarisations has not been upgraded in this sense (see note in Sect. 6), but its semi-analytical implementation is more detailed than in PK20 in terms of time evolution, as it is calculated for each time step and not just for a few selected moments.

Given the 3D SN-asymmetric disc-like CSM interaction simulations (e.g. Blondin et al. 1996; Orlando et al. 2024, and a handful of others), we can anticipate differences or limitations inherent in current 2D models compared to potentially complete 3D models. Given that three dimensions imply more degrees of freedom, it is reasonable to expect a somewhat lower expansion velocity of the SN shock wave, as well as a more complex structure of instabilities, filaments, etc., particularly at the SN ejecta-CSM interface, but also within the expanding envelope. For example, 3D simulations show that Rayleigh-Taylor (RT) instabilities develop more rapidly and are more pronounced than in 2D models. Emission spectra may also depend more variably on the viewing direction in 3D simulations; whereas 2D models show that emission from the SN-disc interaction is clearly influenced by the viewing angle, 3D simulations illustrate in greater detail how emission spectra can be blurred or swallowed up when the ejecta surrounds the disc.

Fig. 7 Bolometric and visual LCs of selected models. The thick blue and brown lines show the numerically calculated pole-on viewed LCs of the models DMIN and LMOD (cf. Figs. 3 and 5). The red and green lines show the LCs of the SN-disc and SN-lobes interactions calculated semi-analytically in PK20. The graph compares the synthetic LCs from our models calculated as a ‘first iteration’ within the planned ‘chain’ of models with the observed LCs of strongly interacting SNe of type IIn, SN1988Z (Aretxaga et al. 1999), iPTF14hls (Arcavi et al. 2017), iPTF13z (Nyholm et al. 2017), and 2006jd (Stritzinger et al. 2012), shown with points.

7.2 Comparison with observed SNe

The LCs in the case of SN-disc interactions, calculated in our DMIN and DMAX models, are basically consistent with observations: both the bolometric and visual LCs show a relatively smooth and slow evolution over time; both curves in the DMIN model (left panel in Figs. 3 and 7) reach maxima around 10⁴³-10⁴⁴ergs⁻¹, while in the DMAX model (right panel in Fig. 3), with a disc density two orders of magnitude higher, the LCs rise steadily over the entire 2-year period, with the maxima occurring roughly one order of magnitude higher. Given the values calculated from the current model, it does not appear that the combination with the contribution of radioactive decay (black solid line in Fig. 3) significantly changes the resulting curve. The LCs of SN-lobe interactions (models LM) show more or less distinct drops and increases. Figure 7 compares the calculated LCs, derived only from the DMIN and LMOD models for an easier overview, with some actual observed curves of prominent interacting SNe IIn, as well as with our previous semi-analytical calculations performed in PK20. In particular, the LCs of these two models show a rough similarity to those determined semi-analytically in PK20. Since the initial densities of the aspherical CSM were estimated somewhat ad hoc, they may be, for example, too high compared to actual situations. The purpose of this and subsequent studies is to gradually determine, based on comparisons, the parameters and morphology of the interacting environment that best correspond to the observed reality. For further considerations related to the course, evolution, and comparison of LCs with actually observed data, see Sects. 4 and 7.1.

Within all models, we also calculate a synthetic profile of broad-band spectra, tracing the appearance or absence of spectral lines of a few selected important elements (see Sect. 5.1 and the top two rows of Fig. A.4) and a typical hydrogen spectral line from three different viewing polar angles (represented here by the Ly-α line; see Sect. 5.2 and the top two rows of Fig. A.5). Few examples of observations of broad-band spectra include, e.g. Bilinski et al. (2018, 2020) in SN2013L and SN2014ab for H lines, Mauerhan et al. (2017) in SN2013ej with the evolution of Ca and O lines, or Taddia et al. (2020) in SN2013L with the evolution of Ca and Fe lines; all of these are in the visible and IR range. Among the most recent are, for example, Kangas et al. (2022) for more H-rich SLSNe in the visible and near IR range, or Aamer et al. (2025) in the UV range of SN 2023taz, although the latter concerns hydrogen-poor SNe that show certain early evolutions of, for example, C, Mg, Ca, or Fe lines. The evolution of individual spectral lines (typically Hα) is documented, in addition to the articles already mentioned, for example, in Andrews & Smith (2018) for iPTF14hls, Kokubo et al. (2019) for KISS15s, Szalai et al. (2019); Weil et al. (2020) for SN2017eaw, and Andrews et al. (2019) for SN2017gmr; see PK20 for a more comprehensive list. One of the important features of the spectral lines in this context is the double-peaked spectral profile associated with the disc or torus-like morphology of the surrounding CSM (e.g. Gerardy et al. 2000; Jerkstrand 2017). Such a double-peaked profile is evident in our DM models, especially when viewed from the pole, and to a certain extent also in LM; see Fig. A.5. The DMAX model with a dense disc also shows significant central absorption when viewed equatorially, i.e. a typical ‘shell’ profile. The square profile (very approximate) with a flat peak and an indication of double-peaked evolution arises from bipolar lobes similar to those seen in η Car. An example may be the eruption of the impostor SN UGC 2773-OT (Smith et al. 2016). Our models LM in some evolutionary stages in Fig. A.5 clearly resemble such profiles, with the exact shape of the spectral profiles depending on the viewing angle.

All models show a clear blue-red asymmetry of the spectral line in Fig. A.5, at almost all times and from all viewing directions. The transition from the initially slightly predominant blue peak to the later dominant red peak is striking. This phenomenon is well observed and documented, for example, in Smith et al. (2015) for PTF11iqb, including a discussion of possible connections with the presence of dust in the vicinity of an SN. Our current models indicate that this typical bluered asymmetry in the Ly-α line should develop regardless of the dusty environment and irrespective of rotational asymmetry (which was an assumption discussed in PK20). We note that no dust is included in the spectral calculations; furthermore, the dust, together with the expanding, relatively dense envelope, should tend to block the red part of the spectrum, as seen in many SNe (cf. Smith et al. 2015). This developing asymmetry in the Ly-α line arises in the numerically calculated profile of spectral lines, even in completely axisymmetric models. The blue wing of the line may be stronger in the early phases because the hydrogen envelope has not yet expanded enough to develop a high column density or high-speed outflows; this allows ‘bluer’ photons to escape directly rather than being scattered into the red wing. At later times, the photons undergo repeated resonant scattering by neutral hydrogen within the expanded nebula with higher optical depth; they diffuse in frequency to the red side of the line. The red-shifted photons ‘see’ lower opacity, which allows them to escape more easily, creating a stronger red wing (e.g. Yamada et al. 2012; France et al. 2015). On the other hand, the blue-red asymmetry of the Hα spectral line, where the blue wing remains stronger than the red wing over an extended period, does not change significantly in Type II SNe, even in later stages; the emission originates from the inner regions of the ejected material, where the structure is predominantly formed by the initial explosion and remains relatively stable. The presence of the (however weak) line asymmetry indicates strong, ongoing interaction with a dense, asymmetric, H-rich CSM, providing direct evidence of extreme pre-SN mass loss (e.g. Blondin et al. 1996; Reynolds et al. 2025, etc.).

Our estimates of relative polarisation in the optically thin regime (Table 2 and Fig. A.7) are essentially similar to those observed (e.g. Dessart & Hillier 2011; Gal-Yam 2019). Aspherical CSM-discs and bipolar lobes-have a prolate ejecta shape and a maximum polarisation (in the equatorial line of sight) in the range of 0.5-3% for all density levels. In PK20, we determined the spherically symmetric radius of the He core as an optically thin boundary, while here we calculate the actual limit τ = 1 for the viewing direction θ = π/2 (see the explanation in Sect. 6).

8 Summary and future work

We have upgraded the two-dimensional (2.5 D) axisymmetric RHD simulations of SN-aspherical CSM interactions, incorporating a circumstellar disc and bipolar lobes with different mass levels. The typical CSM masses within our grid are 0.3 to 2.1 M_⊙ in the five presented models. Confrontation with observations will lead to corrections in the next stage of modelling. Snapshots of density, velocity magnitude, temperature, and radiative energy in the models DMIN and LMIN, which have the lowest masses, are shown in Figs. 1 and 2; cf. also Sect. 3. The higher mass RHD models are presented in Figs. A.1, 2, A.2, and A.3. Compared to purely adiabatic models in PK20, current RHD models show greater deceleration of ejecta (~60% of expansion velocity), while other features of gas-dynamics, including the deflection of ejecta material, the development of shock waves, and instabilities, may be regarded as essentially similar.

Following our RHD simulations, the most significant upgrade consists of using the MC-RT code to calculate the following observable quantities of CSM interaction in SNe, depending on the viewing direction: LCs, both bolometric and visual, in selected filters, as well as broadband spectra in the range of 900-9000 Â, with an evaluation of the presence of selected spectral lines, details of spectral line evolution (in this case Ly-α), and relative polarisation (calculated semi-analytically). These calculated observable characteristics are described in Sects. 4, 5, and 6; demonstrated in Figs. 3-6 and A.4-A.7. The observables show profiles and evolution that are generally consistent in the first ‘iteration’ with the observations of Type IIn SNe, as well as with the semi-analytically calculated results of our previous studies PK19 and PK20 (see the plotted comparison of LCs in Fig. 7). The evolution of the individual spectral line (Fig. A.5) shows the expected double-peaked profile in SN-disc interaction models and a flat-topped or multipeaked profile in SN-lobe interaction models; however, even in the spherically symmetric configuration, a remarkable blue-red asymmetry evolves (cf. e.g. the PTF11iqb Hα line evolution in Smith et al. 2015, or the 2018lqi and 2006gy Hα line evolution in Kangas et al. 2022).

In follow-up studies, we plan to further improve and refine the parametrisation and morphology of the SN progenitor and the surrounding CSM based on the agreement of observables with actual events. We also want to extend the computational domain to distances of the order of parsecs. In addition to more realistic luminous-blue-variable lobe scaling, we aim to calculate, for example, the interactions of SNe with nuclear accretion discs, jets, and other structures in the vicinity of galactic nuclei (cf. Czerny et al. 2023a; Kurfürst et al. 2025; Garland et al. 2026).

Data availability

The movies associated with Figs. 1,2, A.1-A.3 are available at https://www.aanda.org

Acknowledgements

We are grateful to the anonymous referee for the constructive comments that helped to improve the manuscript. We thank Filip Hroch for his technical support and useful improvements in the computational process. We greatly appreciate the advice on the Castro code from Michael Zingale and Max Katz. We appreciate the important discussions with Slah Abdellaoui and Jakub Podgorný about polarisation properties. We also thank Dan Kasen for useful information about the code SEDONA, as well as Stuart Sim and Christian Knigge for their advice on the code SIROCCO. PK received support from the OPUS-LAP/GACR-LA bilateral project (2021/43/I/ST9/01352/OPUS22 and GF23-04053L). JKs and BK were supported by grant GA CR 25-15910S. MZ acknowledges the GACR Junior Star grant no. GM24-10599M for support. Computational resources were supplied by the project ’e-Infrastruktura CZ’ (e-INFRA LM2018140) provided within the programme Projects of Large Research, Development, and Innovation Infrastructures.

References

Appendix A Additional graphs and figures

Appendix A.1 RHD simulation-model DMAX

Figure A.1 shows the RHD evolution snapshots of the DMAX model using the same structure as in Fig. 1; for further description, refer to Sect. 3.1. In this model, among other things, it is clearly visible how the dense inner part of the disc is pushed away by the expansion of the SN, creating a gradually thinning, propeller-shaped structure that roughly fills the angular range of the disc’s opening angle. Unlike the DMIN model (remember that the initial density of the aspherical component CSM - the disc - is two orders of magnitude higher than in the DMIN model), the expansion of the SN ejecta in the equatorial direction is almost completely blocked. The reason for this is the rapid loss of momentum in the direction of the disc and its deflection into areas of lower density. This evolution is also reflected analogously in the temperature structure and the structure of radiative energy, while the development of the velocity structure shows the opposite effect.

Fig. A.1 Model DMAX. Selected evolutionary stages of SN interaction with a ’denser’ circumstellar disk structure calculated by the RHD code. The individual columns and rows show the same quantities as in Figs. 1 and 2. The animated version of this image is available online.

Appendix A.2 RHD simulation-model LMOD

Figure A.2 shows the RHD evolution snapshots of the LMOD model, using the same structure as in Fig. 2; for further description, refer to Sect. 3.2. Unlike the LMIN model (remember that the initial density of the aspherical component CSM - the bipolar lobes - is two orders of magnitude higher on average than that in the LMIN model), where the maximum expansion velocity heads towards the equatorial direction, in the LMOD model, the expansion of the SN ejecta in both the equatorial and polar directions is significantly decelerated by the denser structure of the transversely positioned lobe slabs, especially by the much higher mass-concentrated narrow waist. As a result, the highest velocity magnitude appears towards the mid-positioned 45° angular direction, forming a significant protrusion arm there.

Fig. A.2 Model LMOD. Selected stages of SN interaction with a ’moderate density’ bipolar lobes structure calculated by the RHD code. The individual columns and rows show the same quantities as in Figs. 1, 2, or A.1. The animated version of this image is available online.

Appendix A.3 RHD simulation-model LMAX

Figure A.3 shows the RHD evolution snapshots of the LMAX model using the same structure as in Fig. 2; for further description, refer to Sect. 3.2. Since the initial density of the aspherical component CSM - the bipolar lobes - is two orders of magnitude higher than that in the LMOD model on average (four orders of magnitude higher than in the LMIN model), the SN expansion is almost entirely blocked by the lobe shells; this effect practically does not allow the central parts of the expanding SN configuration to rarefy. However, the relative expansion velocity is, in this case, highest towards the pole, where the efficiently deflected mass causes the highest concentration of momentum in this direction. The expansion of the SN ejecta in the equatorial and intermediate directions is almost entirely blocked. The expanding envelope is effectively slowed down by dense lobes where it accumulates for some time, thus emphasising their original structure, which is particularly visible in plots of density, temperature, and radiative energy.

In addition, Fig. A.8 shows the evolution of He and Fe densities in the models DMIN and LMIN at selected time snapshots. The structure of the plots is described in Sects. 3.1 and 3.2.

Fig. A.3 Model LMAX. Selected stages of SN interaction with a ’densest’ bipolar lobes structure calculated by the RHD code. The individual columns and rows show the same quantities as in Figs. 1, 2, A.1, or A.2. The animated version of this image is available online.

Appendix A.4 MC-RT simulations

Figures A.4 and A.5 show the temporal evolution of broadband spectra and the individual spectral line Ly-α of the presented models, respectively, calculated numerically using the MC-RT codes. The structure of these images and the nature of the physics behind them are sufficiently described in Sects. 5.1 and 5.2. In addition, Fig. A.7 shows the semi-analytically calculated relative polarisation evolution; for further description and explanation, refer to Sect. 6.

Fig. A.4 The overall behaviour of the spectral features calculated for all models marked on the right for individual rows, using the MC-RT code SEDONA (Kasen et al. 2006) under conditions described in Sects. 2 and 5. The three columns enhance the difference between the spectra if only H and He are included (left column), some intermediate mass elements C, O, and Si are added (middle column), and Fe with Ni are added (right column). The times of particular spectra snapshots are marked directly in the image next to the corresponding graph. The coloured graphs show the spectra calculated from the direction of the pole, while the paired black graphs are calculated from an equatorial perspective. For clarity, certain pairs of spectra are shifted upwards by the positive values shown next to them. Some significant spectral lines are marked in the graph. Blends of iron lines in the right column are not marked.

Fig. A.5 The time evolution of a hydrogen spectral line. In this case, the detailed calculations of Ly-α line were made using the MC-RT code SIROCCO (Matthews et al. 2025) under conditions described in Sects. 2 and 5. The individual rows show a particular model, labelled on the left; each column corresponds to a different viewing polar angle θ denoted within each panel. The times of shown simulations range from 50 to 730 days with the particular time-snapshots emphasised by different colours and labelled directly in the plot. Spectral line profiles are plotted on a linear scale and normalised for clarity of presentation; for this reason, we do not display units on the vertical axis (cf. fig. 12 in PK20). The originally slightly blue-shifted peak of emission transfers after about 300-500 days to a red-shifted peak, as corresponds to observations of the spectral lines of SNe (see, e.g. Smith et al. 2015).

Fig. A.6 As in Fig. A.5, now the detailed calculations of Hα line. The characteristic central emission peak is significantly narrow while the doublepeaked profile begins to appear only at later time, when the process enters the nebular phase (remarkable after some 200-300 days); the earlier profiles exhibit, except the central narrow peak, only an irregular noise within a continuum. See the further description in Sect. 5.2.

Fig. A.7 Graph of maximum (for θ = π/2) relative polarisation evolution for the five models, including only the optically thin region (τ ≤ 1). The values at selected times are listed in Tab. 2. (See also the polarisation calculations for models in PK20.)

Fig. A.8 Left panel: Model DMIN. Top row : Evolution of He density at two evolved times marked in the graph. The density scaling is the same as in Fig. 1. Bottom row : Evolution of Fe density at the same two times with different scaling. Right panel: Model LMIN. Top row : Example of map of the temporal evolution of helium density at two evolved times marked in the graph. The density scaling is the same as for the total density of this model in the top row of Fig. 2. Bottom row: Map of the temporal evolution of iron density at the same two times, with significantly different density scaling.

Appendix B Analytical scaling and parametrisation of initial states of the models

Appendix B.1 Initial scaling of spherically symmetric wind and circumstellar disc

We insert two components of the gaseous surroundings of an SN as an initial state: the spherically symmetric stellar wind, whose parameters correspond to a typical RSG stellar wind, and an aspherical component; that is, a circumstellar disc or a system of bipolar lobes, whose configuration may resemble the iconic Homunculus nebula around the star η Car (e.g. Gull & Damineli 2010; Smith 2013). Since we do not include the star’s rotation, we do not refer to the disc as equatorial; other designations of directions, such as equatorial or polar, are therefore more or less formal. We set the stellar wind density, velocity, and temperature structure as $${\rho }_{\text{w}}={\rho }_{0,\text{w}}{\left(\frac{R_{\star }}{r}\right)}^2,v_{\text{w}}=15\text{km}{\text{s}}^{-1},T_{\text{w}}=T_{\star ,\text{eff}}{\left(\frac{R_{\star }}{r}\right)}^{0.75}.$$(B.1)

Here, ρ_0,w ≈ 10⁻¹⁵ g cm⁻³ is the wind base density (the density at the star’s surface) that corresponds to the star-wind mass loss rate Ṁ_w = 10⁻⁶ M_⊙ yr⁻¹ (cf. van Marle et al. 2010, PK20), r is the spherical radial distance from the star’s centre, v_w corresponds to the asymptotic RSGs wind velocity (Goldman et al. 2017), and the stellar (surface) effective temperature T_{*, eff} ≈ 3300 K, which is typical for RSGs (all the star-subscripted quantities hereafter refer to the stellar surface values).

We include the initial gas-dynamic structure of the circumstellar α-disc (Shakura & Sunyaev 1973), using the wind temperature profile in Eq. (B.1) for the disc in the following way (cf. Lee et al. 1991; Kurfürst 2015; Kurfürst et al. 2014, 2018), $${\rho }_{\text{d}}={\rho }_{0,\text{d}}{\left(\frac{R_{\star }}{\varpi }\right)}^3\sqrt{\frac{r}{\varpi }}\text{exp}[-2\sqrt{\frac{R_{\star }}{\varpi }}{\left(\frac{R_{\star }}{H_{\star }}\right)}^2(1-\sqrt{\frac{\varpi }{r}})],$$(B.2)

where ϖ denotes the distance from the disc axis ($r=\sqrt{{\varpi }^2+z^2}$. Here, z is the vertical coordinate, H = c_s/Ω is the disc vertical scale height (Krtička et al. 2011),. c_s(r) being the isothermal speed of sound related to the disc temperature T_d(r) ∝ r^−0.75 (as in Eq. (B.1)) via the ideal gas EOS $c_{\text{s}}^2$ = kT/μm_u (where k is the Boltzmann constant, μ is the mean molecular weight, and m_u is the atomic mass unit), $\mathrm{\Omega }(\varpi )=\sqrt{G{M}_{\star }/{\varpi }^3}$ being the (Keplerian) angular velocity of the disc elements, and P_d(r) = ρ_d c_s²(r) is the isothermal ideal gas disc pressure. Since we assume a dominant external source of heating from the central star, we neglect the disc’s internal viscous-heating (Kurfürst et al. 2018). In the initial state setup, we neglect both velocity components, radial and azimuthal; the latter only formally defines Ω(ϖ) in the density prescription in Eq. (B.2). We insert two different initial values of the disc midplane base mass density ρ_0,d (i.e. the density at the disc mid-plane near the star’s surface), 5 × 10⁻¹³ g cm⁻³ in model DMIN (about 10 times higher than in PK20), and 5 × 10⁻¹¹ g cm⁻³ in model DMAX. The p_0,d value for model DMIN corresponds to the disc accretion (or decretion) rate Ṁ ≈ 10⁻¹⁰ M_⊙ yr⁻¹, while model DMAX represents the disc accretion rate two orders of magnitude higher Ṁ ≈ 10⁻⁸ M_⊙ yr⁻¹. The total mass of the spherical wind is in both cases M_w ≈ 10⁻³ M_⊙, while the superposed initial mass of the disc and the underlying wind is M_w+d ≈ 4 × 10⁻¹ M_⊙ in the DMIN model and M_w+d ≈ 2.4 M_⊙ in the DMAX model, which is about 2-3 orders of magnitude higher than that of the spherical wind alone. The parameters of both disc models were chosen to compare two adequately strong SN-CSM interactions with significantly different densities of a similar non-spherical formation.

We set the initial structure of the underlying stellar wind similar to models DM while using the polytropic approximation and the ideal gas law for the temperatures of the expanding lobes in models LM. We set the temperature in the unperturbed CSM to 20 K.

Appendix B.2 Initial scaling of bipolar lobes

For clarity, we list here Eqs. (1)-(7) (in a rather simplified form used in this work) and table 1 of González et al. (2010) used for the setup of bipolar lobe models, referred to in Sect. 2.2. The initial profile of bipolar lobes in models LM cannot be described simply analytically; we implement the semi-analytically pre-calculated structure that corresponds to five fundamental components (five ’colliding winds’ with different density and velocity parameters) of the Homunculus nebula, starting to evolve at five subsequent phases: pre-outburst wind (pew), giant eruption (ge), first post-outburst wind (pow3), minor eruption (me), and second post-outburst wind (pow5). We also paraphrase the models used in PK20 when we reduce the size of the nebula in a self-similar way to fit our current computational grid domain. For each of the five components j, a co-latitudinal velocity v_j dependence is scaled as $$v_j=v_{1,j}{F}_j(\theta )\text{with}{F}_j(\theta )=\frac{\left(v_{2,j}/v_{1,j}\right)+{\text{e}}^{2z}}{1+{\text{e}}^{2z}},$$(B.3)

where θ is the polar angle (θ = 0 at the pole, θ = π/2 at the equator), Z_j = λ_j cos (2θ), and the parameters v_1,j, v_2,j (velocities related to the polar and equatorial directions, respectively) and λ_j are given in Tab. B.1. The initial mass density of the components is $${\rho }_j={\rho }_0{\left(\frac{r_0}{r}\right)}^2\frac{1}{F_j(\theta )},$$(B.4)

where ${\rho }_0={\dot{m}}_0/(4\pi v_0{r}_0^{\text{2}})$ is the density at the reference radius r₀ (selected as ≈8 × 10¹⁴ cm ≈ 10R_*), the downstream velocity v₀ = 250 kms⁻¹ at r₀ is adopted as a free parameter, and Ṁ₀ = 10⁻⁵ M_⊙ yr⁻¹ is the total mass loss rate of the wind in the model LMIN. Besides that, we parameterise the mass loss rates of eruptive events in the model LMIN as Ṁ_ge = 5 × 10⁻³ M_⊙ yr⁻¹ and Ṁ_me = 10⁻⁴ M_⊙ yr⁻¹. The parameters Ṁ₀, Ṁ_ge, and Ṁ_me are in the models LMOD and LMAX increased by 10², respectively (see Tab. B.1). When the great eruption begins (at phase 2), the wind parameters increase to av₀ and bm₀, with a and b being constants. The flow parameters form a working surface at the base of the wind. It corresponds to the large Homunculus lobe (formed in the ’ge’ phase) that expands with a constant velocity $$v_{\text{ws}}=\sigma v_0,$$(B.5)

with σ = (a^1/2 + ab^1/2)/(a^1/2 + b^1/2). The wind velocity is intermediate between the low-velocity downstream wind v₀ and the faster upstream outflow av₀. We have also set five time limits: t₁ = 25 yr, t₂ = 50 yr, t₃ = 87.5 yr, t₄ = 100 yr, and t₅ = 112.5 yr, to separate five different courses of the nebula expansion behaviour. The two time intervals ∆t₂ = 30 yr and ∆t₄ = 10 yr represent the duration of the post-eruption phases. The mass loss rates Ṁ₃ and Ṁ₅ in the phases after the giant and minor eruptions, respectively, are $${\dot{m}}_3={\dot{m}}_0\varphi (t)$$(B.6)

and $${\dot{m}}_5={\dot{m}}_0\phi (t),$$(B.7)

where we adopt two functions $$\varphi (t)=\frac{{\dot{m}}_{\text{ge}}}{{\dot{m}}_0}+(1-\frac{{\dot{m}}_{\text{ge}}}{{\dot{m}}_0})\sin \left[\frac{\pi }{2}\left(\frac{t-t_2}{\mathrm{\Delta }{t}_2}\right)\right]$$(B.8)

and $$\phi (t)=\frac{{\dot{m}}_{\text{me}}}{{\dot{m}}_0}+(1-\frac{{\dot{m}}_{\text{me}}}{{\dot{m}}_0})\sin \left[\frac{\pi }{2}\left(\frac{t-t_4}{\mathrm{\Delta }{t}_4}\right)\right],$$(B.9)

which, due to the fact that different Ṁ form fractions, remain the same for all LM models (see Tab. B.1). The computational domain is initially filled with a homogeneous medium with a temperature of T_a = 10² K and a density of n_a = 10⁻²⁵ g cm⁻³. This is followed by the injection of a pre-outburst wind with a mass-loss rate of Ṁ₀ at a distance of r₀, with a temperature of T₀ = 10⁴ K. Afterwards, a gaseous distribution is gradually formed within the time limits according to the prescriptions and parameters introduced in this section.

All Tables

All Figures

Fig. 1 Model DMIN. Selected evolutionary stages of SN interaction with a ’thinner’ circumstellar disc calculated by the RHD code. The four columns show the snapshots at times t = 0, 240, 480, and 730 days. The individual rows shows four different quantities: from top to bottom it is density ρ, velocity magnitude v, temperature T, and radiative energy density Er. The resolution of the bitmaps had to be significantly reduced compared to resolution of the simulations in order to reduce the file size. The animated version of this image is available online.

In the text

Fig. 2 Model LMIN. Selected evolutionary stages of the SN interaction with the ’thinnest’ bipolar lobes structure calculated by the RHD code. The four columns show the snapshots at times t = 0, 240, 480, and 730 days. The individual rows show four different quantities: from top to bottom it is density ρ, velocity magnitude v, temperature T, and radiative energy density Er. The animated version of this image is available online.

In the text

Fig. 3 Pseudo-bolometric, directionally ‘observed’ LCs of the two models of the SN interaction with circumstellar disc-like structures, denoted as model DMIN (left panel) and model DMAX (right panel), up to approximately 2 years (see Sects. 3.1 and 4.1). The upper pair of lines shows the LCs calculated within the ’largest’ frequency range 1014-2 × 1016 Hz (see Sect. 2.4), while the bottom pair of lines shows the same within the range (3.875-7.825) × 1014 Hz, which roughly corresponds to the frequency range of visible light (for LCs corresponding to AB magnitudes in main filters, see Fig. 4). We demonstrate here the LC profiles for two different viewing angles in the graphs, the pole-on (polar) view marked in blue and the equatorial view marked in dashed black line for the bolometric LC (labelled Lbol in the legend), and the polar view marked in red dash-dotted line and the equatorial view marked in dotted black line for the visual LC (labelled jvis in the legend). The black solid line labelled ’no CSM’ corresponds to an SN without CSM interaction, heated with 0.28 M⊙ of radioactive Ni which roughly corresponds to the brightest normal H-rich SNe (Pejcha & Prieto 2015; Müller et al. 2017) or the population mean of stripped SNe (Anderson 2019).

In the text

	Fig. 4 Calculated AB magnitudes for the two models of the SN interaction with a circumstellar disk-like morphology. Each row in the figure corresponds to one of the models DMIN or DMAX labelled on the right. The three columns demonstrate the three various viewing angles marked at the top. The AB magnitudes are calculated for nine filters listed in the left top and bottom panels.
In the text

	Fig. 5 Pseudo-bolometric, directionally ‘observed’ LCs of the three models of the SN interaction with a bipolar lobes structure, denoted as model LMIN (left panel), model LMOD (middle panel), and model LMAX (right panel), up to approximately 2 years (see Sects. 3.2 and 4.2). The meaning and structure of the image, the order and types of lines, and the legend are the same as in Fig. 3.
In the text

	Fig. 6 Calculated AB magnitudes for the three models of the SN interaction with a bipolar lobes-like morphology. Each row in the figure corresponds to one of the models LM, labelled on the right. The other features of the figure are analogous to Fig. 4.
In the text

Fig. 7 Bolometric and visual LCs of selected models. The thick blue and brown lines show the numerically calculated pole-on viewed LCs of the models DMIN and LMOD (cf. Figs. 3 and 5). The red and green lines show the LCs of the SN-disc and SN-lobes interactions calculated semi-analytically in PK20. The graph compares the synthetic LCs from our models calculated as a ‘first iteration’ within the planned ‘chain’ of models with the observed LCs of strongly interacting SNe of type IIn, SN1988Z (Aretxaga et al. 1999), iPTF14hls (Arcavi et al. 2017), iPTF13z (Nyholm et al. 2017), and 2006jd (Stritzinger et al. 2012), shown with points.

In the text

	Fig. A.1 Model DMAX. Selected evolutionary stages of SN interaction with a ’denser’ circumstellar disk structure calculated by the RHD code. The individual columns and rows show the same quantities as in Figs. 1 and 2. The animated version of this image is available online.
In the text

	Fig. A.2 Model LMOD. Selected stages of SN interaction with a ’moderate density’ bipolar lobes structure calculated by the RHD code. The individual columns and rows show the same quantities as in Figs. 1, 2, or A.1. The animated version of this image is available online.
In the text

	Fig. A.3 Model LMAX. Selected stages of SN interaction with a ’densest’ bipolar lobes structure calculated by the RHD code. The individual columns and rows show the same quantities as in Figs. 1, 2, A.1, or A.2. The animated version of this image is available online.
In the text

Fig. A.4 The overall behaviour of the spectral features calculated for all models marked on the right for individual rows, using the MC-RT code SEDONA (Kasen et al. 2006) under conditions described in Sects. 2 and 5. The three columns enhance the difference between the spectra if only H and He are included (left column), some intermediate mass elements C, O, and Si are added (middle column), and Fe with Ni are added (right column). The times of particular spectra snapshots are marked directly in the image next to the corresponding graph. The coloured graphs show the spectra calculated from the direction of the pole, while the paired black graphs are calculated from an equatorial perspective. For clarity, certain pairs of spectra are shifted upwards by the positive values shown next to them. Some significant spectral lines are marked in the graph. Blends of iron lines in the right column are not marked.

In the text

Fig. A.5 The time evolution of a hydrogen spectral line. In this case, the detailed calculations of Ly-α line were made using the MC-RT code SIROCCO (Matthews et al. 2025) under conditions described in Sects. 2 and 5. The individual rows show a particular model, labelled on the left; each column corresponds to a different viewing polar angle θ denoted within each panel. The times of shown simulations range from 50 to 730 days with the particular time-snapshots emphasised by different colours and labelled directly in the plot. Spectral line profiles are plotted on a linear scale and normalised for clarity of presentation; for this reason, we do not display units on the vertical axis (cf. fig. 12 in PK20). The originally slightly blue-shifted peak of emission transfers after about 300-500 days to a red-shifted peak, as corresponds to observations of the spectral lines of SNe (see, e.g. Smith et al. 2015).

In the text

Fig. A.6 As in Fig. A.5, now the detailed calculations of Hα line. The characteristic central emission peak is significantly narrow while the doublepeaked profile begins to appear only at later time, when the process enters the nebular phase (remarkable after some 200-300 days); the earlier profiles exhibit, except the central narrow peak, only an irregular noise within a continuum. See the further description in Sect. 5.2.

In the text

	Fig. A.7 Graph of maximum (for θ = π/2) relative polarisation evolution for the five models, including only the optically thin region (τ ≤ 1). The values at selected times are listed in Tab. 2. (See also the polarisation calculations for models in PK20.)
In the text

Fig. A.8 Left panel: Model DMIN. Top row : Evolution of He density at two evolved times marked in the graph. The density scaling is the same as in Fig. 1. Bottom row : Evolution of Fe density at the same two times with different scaling. Right panel: Model LMIN. Top row : Example of map of the temporal evolution of helium density at two evolved times marked in the graph. The density scaling is the same as for the total density of this model in the top row of Fig. 2. Bottom row: Map of the temporal evolution of iron density at the same two times, with significantly different density scaling.

In the text