Diagnosis of Circumstellar Matter Structure in Interaction-powered Supernovae with Hydrogen Line Features

Some supernovae (SNe) are powered by the collision of the SN ejecta with dense circumstellar matter (CSM). Their emission spectra show characteristic line shapes of combined broad emission and narrow P Cygni lines, which should closely relate to the CSM structure and the mass-loss mechanism that creates the dense CSM. We quantitatively investigate the relationship between the line shape and the CSM structure by Monte Carlo radiative transfer simulations, considering two representative cases of dense CSM formed by steady and eruptive mass loss. Comparing the Hα emission between the two cases, we find that a narrow P Cygni line appears in the eruptive case but does not appear in the steady case due to the difference in the velocity gradient in the dense CSM. We also reproduce the blueshifted photon excess observed in some Type IIn SNe, which is formed by photon transport across the shock wave, and find the relationship between the velocity of the shocked matter and the amount of blueshift of the photon excess. We conclude that the presence or absence of narrow P Cygni lines can distinguish the mass-loss mechanism and suggest high-resolution spectroscopic observations with λ/Δλ ≳ 104 after the light-curve peak for applying this diagnostic method.

1. INTRODUCTION Supernovae (SNe), which are explosions following the deaths of massive stars, display a diversity in their light curves and spectra that reflects their diverse environment.A fraction of SNe exhibit narrow lines in their spectra, and such SNe are thought to be powered by the interaction between dense circumstellar material (CSM) and SN ejecta (Grasberg & Nadezhin 1986;Chugai 1991Chugai , 2001;;Chevalier & Irwin 2011).The CSM expands much slower than the SN ejecta, and is believed to produce the narrow lines in the spectra due to the Doppler effect of its line emission.
Estimations from light curves (e.g., Moriya et al. 2013Moriya et al. , 2014) ) and spectra (e.g., Kiewe et al. 2012;Taddia et al. 2013) of Type IIn SNe (SNe IIn), whose spectra exhibit narrow lines of hydrogen, indicate that their progenitors should have lost mass at a high mass-loss rate of Ṁ ≳ 10 −3 M ⊙ /yr to form the dense CSM.This is in stark contrast to late time radio observations of SNe with mass-loss rate estimates of Ṁ = 10 −7 − 10 −5 M ⊙ /yr (Weiler et al. 2002;Chevalier et al. 2006;Bietenholz et al. 2021), compatible with those observed in stellar winds of local evolved massive stars (e.g., de Jager et al. 1988;van Loon et al. 2005;Mauron & Josselin 2011;Beasor et al. 2020).Therefore, progenitors of SNe IIn should have undergone a different stage of enhanced mass loss just before the SN explosion 1 .
One of the plausible formation processes of such dense CSM is the mass ejection from a star prior to the SN explosion.Existence of such violent mass ejections was pointed out by Chugai et al. (2004), by modeling of Hα lines in the SN IIn 1994W using radiative transfer computations.Indeed, observations detected precursor outbursts in a significant fraction of SNe IIn (Ofek et al. 2014a;Strotjohann et al. 2021), the most famous one being SN 2009ip (e.g., Pastorello et al. 2013;Margutti et al. 2014;Smith et al. 2014).The outbursts should be strongly linked to the formation of the dense CSM (Ofek et al. 2014a;Matsumoto & Metzger 2022), since the precursor emission (10 40 -10 41 erg s −1 ) is much brighter than the Eddington limit of a massive star.However, the enhanced mass loss may also be reproduced by steady processes, such as constant near-surface energy deposition (Quataert et al. 2016), pulsation (e.g., Yoon & Cantiello 2010) or core neutrino emission (Moriya 2014).While observations are in great progress, we do not have a clear criterion to distinguish between steady and eruptive mass ejection from the observed features of the CSM.
A clue to distinguishing these mechanisms may be the structure of the resulting CSM.While the steady mass loss simply predicts a CSM with a density profile close to ρ(r) ∝ r −2 with constant velocity, the CSM from eruptive mass-loss has recently been investigated in detail.Recent numerical and analytical studies (Kuriyama & Shigeyama 2020, 2021;Leung & Fuller 2020;Linial et al. 2021;Tsuna et al. 2021b;Ko et al. 2022;Tsang et al. 2022;Tsuna & Takei 2023) have shown that the CSM consists of a bound and unbound part, whose masses strongly depend on the injected energy that triggers the eruption 2 .After several dynamical timescales of the progenitor, the inner region of the CSM gravitationally bound to the star has a shallower density profile with ρ(r) ∝ r −1.5 , and falls back to the star with a velocity of v(r) ≈ − 2GM * /r, where G denotes the gravitational constant and M * is the stellar mass.In contrast to steady-state mass loss with a nearly constant velocity, the CSM created by mass eruption is expected to have a positive velocity gradient, with the outermost region expanding homologously (v(r) ∝ r).
Such differences of the CSM structure in the two mechanisms may be observable in the line emission.Especially, in SNe IIn, the Hα lines have a characteristic shape, with a combination of a broad component and a narrow component with a P-Cyg profile (Salamanca et al. 1998;Pastorello et al. 2002;Kiewe et al. 2012;Zhang et al. 2012;Taddia et al. 2013;Fransson et al. 2014;Taddia et al. 2020).Some observed Hα lines in the late phase show a blue-shifted excess resulting in an asymmetric shape, which might originate from a region with strong radiative cooling immediately behind the shock front (cool dense shell; Dessart et al. 2015), or from the aspherical structures of the CSM (Smith et al. 2015;Kokubo et al. 2019).
The origin of these spectral components has been investigated in several previous works.Chugai (2001) and Huang & Chevalier (2018) have shown by Monte Carlo radiative transfer simulations that the broad component is formed as a result of multiple electron scattering for an optically-thick CSM (or by radiative acceleration at just after shock breakout; Kochanek 2019), while it is attributed to the bulk motion of the shocked region in the optically thin case.However they used a simple model for the shock propagation in the CSM, and also did not include line absorption in their opacity that characterizes the narrow P-Cyg feature.By Monte Carlo simulations including the narrow P-Cyg profile, Chugai et al. (2004) better reproduced line shapes as seen in the spectrum of SN 1994W assuming homologously expanding CSM rather than CSM formed by steady mass loss, and argued that the progenitor underwent an explosive mass loss event (see also Dessart et al. 2009).The simulations of SNe IIn by Dessart et al. (2015) have solved equations of radiative transfer by the Λ-iteration method, taking into account effects of the velocity gradient on the spectral shape of lines (see also Groh 2014;Boian & Groh 2019, 2020 for applications to early phase spectra).While the simulations generally reproduced the observations of a particular super-luminous Type IIn SN 2010jl, they were limited to CSM of a steady wind-line profile with a constant initial velocity.The relationship between the line shape and matter structure has also been investigated for Luminous Blue Variables (LBVs) in Groh & Vink (2011), which shows that abrupt variation in velocity profile is necessary to reproduce the Hα line shape characterising LBVs.
In this study, we quantitatively investigate the relationship between line shapes and CSM structures in order to establish a method to diagnose the CSM structure from the line shapes.We investigate the Hα line as a first step because it is prominently observed in a large fraction of interaction-powered SNe, including SN 1998S which we focus on as a reference for modeling.The propagation of a shock wave in the CSM is calculated by a radiation hydrodynamics code CHIPS (Takei et al. 2022(Takei et al. , 2023) ) for both cases of eruptive and steady mass loss, and the spectra are calculated using a Monte Carlo radiative transfer code as post-process.Then, the shapes of Hα lines in the spectra are compared between eruptive and steady mass-loss cases.We adopt a Monte Carlo approach, which investigates the shape of the Hα line that reaches the observer by following the trajectory of each packet of photons (e.g., Lucy 2005).The Monte Carlo approach is appropriate to diagnose the origin of each feature in the spectra, because one can trace back the history of the packets that form the feature.This paper is organized as follows.We present our models for the simulations and the detailed computational methods in Section 2. In Section 3 we present the spectra obtained by our simulations, with detailed discussion on the interplay of the processes of line emission, scattering and absorption.We furthermore compare the line shapes from the CSM formed by eruptive and steady mass-loss scenarios.We discuss our results in the context of observed Type IIn SNe and summarize this paper in Section 4.

METHODS AND MODELS
Our model for spectral calculations solves the radiative processes around the forward shock front that separates the shocked gas and the surrounding unshocked dense CSM, as depicted in the upper side of Figure 1.The structure of each region as a result of ejecta-CSM interaction is first obtained using the CHIPS code (Takei et al. 2022(Takei et al. , 2023)), which has been demonstrated to be able to simulate both the formation of the dense CSM and the subsequent interaction between the SN ejecta and the CSM.
Using CHIPS we perform radiation hydrodynamics simulations of the interaction between the ejecta and the dense CSM, with the latter considering two scenarios of enhanced mass-loss, eruptive or steady mass loss.Using the hydrodynamical outputs from CHIPS as background, we then perform post-process Monte-Carlo simulations to obtain the spectrum at each epoch, focusing on the features of the Hα emission.The model setup for each case is described in detail below.

Simulating Type IIn Supernovae by CHIPS
We first obtain a representative SN progenitor model with the stellar evolution code MESA version 12778 (Paxton et al. 2011(Paxton et al. , 2013(Paxton et al. , 2015(Paxton et al. , 2018(Paxton et al. , 2019;;Jermyn et al. 2023).We use the example_make_pre_ccsn test suite, and evolve until core-collapse a non-rotating star with zero-age main-sequence mass of M ZAMS = 20M ⊙ .The chemical compositions of hydrogen, helium and metals at ZAMS are set to be X = 0.718, Y = 0.268 and Z = Z ⊙ (= 0.014) respectively.The resultant progenitor is a red supergiant (RSG) that has a photospheric radius of R * = 843R ⊙ , luminosity of L * = 1.75 × 10 5 L ⊙ , and effective temperature of T eff = 4062 K.We have adopted the "Dutch" wind mass-loss rates of Nieuwenhuijzen & de Jager (1990), which has reduced the mass of the progenitor to M * ≈ 18.3M ⊙ upon core-collapse.As the structure of the hydrogen envelope hardly changes in the last decades of its life, we adopt the progenitor at core-collapse as the progenitor before enhanced mass loss.

Models of the Dense Circumstellar Material
We consider two scenarios of the enhanced mass loss from this progenitor that forms the dense CSM, via eruptive (model E) and steady mass loss (model S).The lower side of Figure 1 shows a schematic picture of the evolution of a SN in a dense CSM for both of the two models.A massive star generates the dense CSM through enhanced mass loss prior to an SN explosion.Once the massive star explodes as an SN, the SN ejecta collides with the CSM, generating shock-powered emission along with the narrow lines from reprocessing in the CSM.The remainder of this sub-subsection describes the settings shared by model E and model S, and the settings that differ for each model, respectively.
• Shared parameters for both models Some model parameters characterizing the CSM and SN ejecta share the same values in both models for ease of comparison.The main parameters that characterize the dynamics of the CSM interaction are the mass and radial extent of the CSM.As described in detail below, we adopt the same extent of r CSM ≈ 3 × 10 15 cm and total mass of M CSM ≈ 1.8M ⊙ .The subsequent SN is assumed to leave a neutron star remnant of M rem = 1.4 M ⊙ .The SN ejecta is assumed to have a kinetic energy of 10 51 erg, with ejecta mass of M * −M CSM −M rem ≈ 15.1M ⊙ .Following common prescriptions (Matzner & McKee 1999), the ejecta is set to be homologous with a double power-law density profile (see equation 6 of Takei et al. 2022), with outer and inner power-law indices of n = 11.43 and δ = 1 respectively.

• Eruptive mass loss (model E)
We generate the CSM with the mass eruption module of the CHIPS code (for details see Kuriyama & Shigeyama 2020), by injecting energy into the base of the star's outer envelope and following its response.We inject an energy of 6.0 × 10 47 erg, which is 80% of the envelope's binding energy.The SN is assumed to happen 10 years after the energy injection, which results in a dense CSM of M CSM ≈ 1.8 M ⊙ with a dense part extending out to r CSM ≈ 3 × 10 15 cm.The values chosen generally reproduce the light curve morphology (rise time and peak luminosity) typically seen in SNe IIn (Takei et al. 2022).The velocity of the resulting CSM monotonically increases with radius, with the velocity in the outermost part of the dense CSM of ≈ 100 km s −1 which is comparable to the progenitor's escape velocity.Before the eruptive mass loss, the star should have been losing mass by blowing a steady wind.To include this, we smoothly connect a stellar wind component out to r out = 3 × 10 16 cm the same as Takei et al. (2022), with a mass-loss rate of Ṁ ≈ 1.1 × 10 −5 M ⊙ /yr and velocity of ≈ 24 km/s.These values are obtained from the formulae in Takei et al. (2022) and are rather arbitrary, but the density in this region is so small that it would not affect the results of the subsequent radiative transfer simulations.

• Steady mass loss (model S)
The progenitor is supposed to suddenly start blowing a much denser wind than before, and continue blowing at a constant mass-loss rate till the progenitor explodes.We adopt a double power-law density profile for the CSM, similar to the shape of model E but with the inner region of a wind profile of ρ ∝ r −2 , where the inner and outer power-law indices are set to n in = 2 and n out = 12, respectively.Similar to model E, we smoothly connect a tenuous stellar wind of mass-loss rate Ṁ ≈ 1.1 × 10 −5 M ⊙ /yr.The normalization for ρ is chosen so that the total CSM mass is the same as model E. The largest difference from model E is the initial velocity profile; the CSM velocity is assumed to be 24 km s −1 independent of radius.This is slower than the outermost part of the dense CSM in model E by a factor of 4, and thus corresponds to the enhanced mass loss starting ≈ 40 years before core-collapse3 .

Modelling the Circumstellar Interaction
With these CSM profiles as inputs, we model collisions with the SN ejecta using the SN module of the CHIPS code (for details see Takei & Shigeyama 2020;Takei et al. 2022).The code is split into two steps, the shock dynamics and radiative transfer in the unshocked CSM.The former solves the propagation of the forward and reverse shocks formed by the collision, as well as the detailed structure of the shocked region taking into account radiative transfer in the flux-limited diffusion approximation (Levermore & Pomraning 1981).The radiative flux at the forward shock, emitted outwards in the unshocked CSM, is simultaneously obtained.
With the time-evolving forward shock luminosity as input, the second step solves radiation transfer in the unshocked CSM that is assumed to be a two-temperature fluid, where the gas and radiation temperatures separately evolve with time.The number of computational cells in the unshocked CSM is set to 1000 at the beginning of the interaction of the SN ejecta with the CSM, and decreases with time as the CSM is engulfed by the shock.The width of the narrow lines, which is characterized by the velocity of the unshocked CSM, can be affected by acceleration due to radiation produced from shock interaction (Chevalier & Irwin 2011;Tsuna et al. 2023).We thus add radiative acceleration a of the unshocked CSM at each radii r as a(r) = κ(r)L(r)/4πr 2 c, where κ(r) is the total (scattering + absorption) opacity adopted in the CHIPS code, L(r) is the incoming luminosity, and c is the speed of light.The opacities in the calculation are referenced from opacity tables, generated for (mass-averaged) hydrogen, helium, and metal abundances in the CSM of X = 0.6402, Y = 0.3457, and Z = 0.014 respectively.
Figure 2 shows the temporal evolution of the forward shock generated by the interaction: the shock radius r fs , the effective temperature at r = r fs , and the shock velocity together with the fluid velocity at the immediate downstream.The shock in model E has a higher velocity than in model S due to the slightly lower density in the innermost region of the CSM.In the bottom-left panel of Figure 2, shock velocity increases at the earliest time for both models and then begins to monotonically decrease, which is due to the initial setting of the CHIPS simulations.
The bottom-right panel of Figure 2 shows the bolometric light curves of models E and S generated by CHIPS, measured at r out = 3 × 10 16 cm.The CSM interaction generally dominates the light curve when the shock is running Figure 2. The time-dependent forward shock radius (top-left), effective temperature at the shock front (top-right), the fluid velocity inside the shock wave accompanied with the forward shock velocity (bottom-left), and bolometric light curves (bottomright) for models E and S. We also show for reference the light curves of two well-known Type IIn SNe, 1998S and 2009ip.The x-axis of the light curves are shifted so that the peak of luminosity is at t = 0.
through the dense CSM.As the dense CSM in these two models has similar mass and extent, the resulting light curves are similar through day 80.The light curves from our models have rise time (≈ 15 days) and peak luminosity (≈ 10 43 erg s −1 ) that are representative in the parameter space of SN IIn (Ofek et al. 2014b;Nyholm et al. 2020), and are similar to the well-known SNe like 1998S (Fassia et al. 2000) and 2009ip (Margutti et al. 2014) overlaid in the plot.
Figure 3 shows distributions in the shocked region and the CSM of gas density, velocity, gas temperature, and the number density ratio n Hα /n H of excited hydrogen that causes Hα absorption among the total hydrogen, at days 10/12 (start of the spectrum simulation for model E/S), 30, and 60 after the SN explosion.The ratio n Hα /n H is calculated assuming local thermodynamical equilibrium, with temperature equal to the gas temperature obtained from the CHIPS simulation.The number density n Hα in the fluid comoving frame of excited hydrogen atoms with n = 2 that can cause Hα is given as where n HI is the number density of neutral hydrogen atoms in the fluid comoving frame, g n is the statistical weight of the level n, ν n is the frequency corresponding to the energy difference between level n and the ground level, and T is the gas temperature.We actually take the summation in the denominator of the right hand side only from n = 1 to n = 3.As explained in detail in the end of Section 2.2, the degree of ionization is referenced from a pre-computed table.This creates small fluctuations of n Hα /n H as can be seen in the figure, but this is found to have a negligible effect on the subsequent radiative transfer calculations.For model E at day 10, the unshocked CSM at r ≈ 5 × 10 14 cm is falling back with a negative velocity, hence does not appear in the bottom-left panel.
The shock breaks out at 10 days from explosion in model E, while it breaks out on 12 days in model S because of its slightly larger optical depth of the CSM.As Figure 3 shows, while the unshocked CSM in both of the models have similar density structures, they exhibit different velocity distributions.Such different velocity profiles are preserved even after the CSM is radiatively accelerated by photons from the interaction.As demonstrated in Section 3, this difference results in distinct features of the Hα spectra in the two models.

Post-process Monte-Carlo Radiative Transfer Simulations
The CHIPS code computes light curves from interaction-powered SNe, but it cannot compute detailed spectra.To obtain them, we perform post-process radiative transfer simulations assuming steady-state background flow-field, using the snapshots of the density, velocity and temperature profiles in the shocked region and the unshocked CSM from Section 2.1.The steady-state assumption is valid during the post-breakout interaction phase, as the hydrodynamical properties do not significantly change within the diffusion timescale of photons in the relevant regions, except for the early phase soon after the interaction of SN with CSM occurs.
The computational domain and radiation processes included in the simulation are shown in the upper part of Figure 1.We illustrate what kind of radiation process is important in which part of the computational domain, and how each process affects the line shape of the resulting spectrum (details are described in Section 3).Photons from the shock front are assumed to be emitted as blackbody radiation, since the photons in the optical (and UV) wavelength range are expected to be immediately thermalized in the shocked CSM with a large mass-loss rate of Ṁ ≳ 0.1M ⊙ yr −1 (v CSM /100 km s −1 ) (Chevalier & Irwin 2012;Svirski et al. 2012;Tsuna et al. 2021a;Margalit et al. 2022).We further incorporate hydrogen recombination emission and free-free emission in the unshocked CSM.Photons experience the collisional processes such as electron scattering, Hα absorption and re-emission, and free-free absorption.We use a Monte Carlo method for radiation transfer simulations (Lucy 1999(Lucy , 2005;;Maeda et al. 2014;Ishii et al. 2017).The fundamental part of photon transport and electron scattering process are following Ishii et al. (2017), in which transport in a relativistic flow-field is considered by taking into account Compton scattering, while thermal motion of electrons is ignored.In this study, the Doppler effect by the thermal motion of electrons following a Maxwell distribution with gas temperature is taken into account for scattering because this effect is important in determining the line shape.In addition the recombination emission, free-free emission/absorption, and line absorption/re-emission processes are installed in this study.A "packet" of photons with the same frequency is treated as a sample particle.We use a total of 10 8 packets for calculating the spectrum at each epoch, which is verified to be enough for convergence.
In this study we focus only on the line shape of Hα, and thus the initial photon frequency distribution is limited to the vicinity of the Hα transition frequency.Packets for the continuum photons are initially assigned to a frequency in the range of 0.9ν Hα ≦ ν 0 ≦ 1.2ν Hα in the fluid comoving frame, where ν Hα = c/(656.28nm) is the Hα photon frequency, while packets for Hα photons via recombination emission process are assigned to ν 0 = ν Hα .Even if the energy distribution in the fluid comoving frame is limited to around Hα, it is broadened in the observer frame depending on the photon direction and the fluid velocity.Since a large number of photon packets is required to resolve the characteristic narrow lines in interaction-powered SNe as well as the broadened (often called "intermediate-width") lines, we need to focus on such narrow energy range to save the computational cost.
The computational region consists of the shocked and unshocked CSM.As the profiles of the physical quantities in the shocked region are found to be nearly constant in radius, we treat this region as one zone for simplicity.At a given time, we determine quantities of this zone by averaging the radial profile of the shock dynamics calculations in CHIPS.First, the average density ρs and internal energy ūs are calculated by a volume weighted average in the shocked region.The average temperature Ts is then calculated from where µ(ρ s , Ts ) is the mean molecular weight obtained by solving Saha equations with hydrogen and helium, m u is the atomic mass unit, k B is the Boltzmann constant, and a r is the radiation constant, respectively.The average velocity vs is obtained by dividing the total momentum of the shocked region by its total mass.The width of this zone ∆r s is determined so that ∆r s = ∆τ s /n e σ T , where σ T is the Thomson scattering cross section and ∆τ s is the Thomson scattering optical depth in the shocked region derived by integration over the region using the density and temperature profiles obtained from CHIPS.The one-zone electron number density ne (ρ s , Ts ) is calculated by solving Saha equations for mixture of hydrogen and helium gas assuming thermal equilibrium.We obtain the ionization degrees of hydrogen and helium by referencing from a pre-computed table for given values of n HI , the number density for neutral helium atoms n HeI in the fluid comoving frame, and gas temperature T .The values of n HI , n HeI , and T in the table cover a broad range of 10 − 10 15 cm −3 , 0.1 − 5 × 10 13 cm −3 , and 10 3 − 10 5 K, respectively.Here, the mass fractions for hydrogen and helium are set to be the same as the CHIPS calculations, and the contribution to n e from heavier elements is neglected for simplicity.

Emission processes
We assume that the radiation emitted from the forward shock front is a blackbody, and thus emit a packet of photons of frequency ν with a luminosity where B ν and h are respectively the Planck function and the Planck constant, and T eff,s is the effective temperature at the shock front.Thermal photons from the shocked region are initially located just 10 −6 % outside of r fs and are emitted outwards.For each Monte-Carlo simulation, both the duration of emitting photons and sampling photons exiting the edge of the computational region at r out = 3 × 10 16 cm are set to 2 × 10 6 s.We set this by the sum of the diffusion time for our CSM models (≈ 10 6 s) and the light-crossing time in the computational region r out /c = 10 6 s.4 Soon after CSM interaction begins, hydrogen atoms in the CSM are ionized by radiation from the shocked region, and recombines with free electrons to emit the Hα line as well as other hydrogen lines.The transition from n > 1 to n = 1 can be neglected, because the emission from such transition is immediately absorbed (case B recombination; Osterbrock & Ferland 2006) for a very dense environment like the CSM in SNe IIn.The case B emission coefficient j Hα (n e , T ) [erg cm −3 ] of the recombination emission in Hα is obtained from Table 4.4 of Osterbrock & Ferland (2006).Since the original table is coarse, we interpolate the values of j Hα as a function of n e and T .The ranges of n e (10 2 -10 6 cm −3 ) and T (5000-20000 K) in the table do not fully cover the ranges required in this study, so we linearly extrapolate this table to the values of n e and T requested5 .The Hα photons are isotropically emitted, with their initial positions determined from a probability distribution proportional to the product of the emission coefficient and the volume at each computational cell.
The frequency ν of the photon packet with free-free emission process in the CSM is determined by the spectral emission coefficient due to bremsstrahlung (e.g.Zel'dovich & Raizer 1967).We consider hydrogen and helium for the ions and neglect the contribution from heavier elements.Similar to the recombination emission, photons are emitted isotropically with their initial positions following the product of the frequency-integrated emission coefficient and the volume at each computational cell.

Transport process
To investigate the spectral shape of the Hα line from the CSM with some velocity gradients, we use the Sobolev approximation as introduced in Sobolev (1960) and used in many previous works (e.g., Tanaka & Hotokezaka 2013;Maeda et al. 2014;Kerzendorf & Sim 2014;Dessart et al. 2015).While the velocity gradient in the unshocked CSM is much smaller than the SN ejecta, the part of the CSM that contributes to Hα absorption is also radiatively accelerated (see Figure 3), which generally creates a large enough velocity gradient to justify this assumption.
The optical depth τ Hα for Hα absorption in the observer frame is obtained as (for details see Appendix A), where f is the oscillator strength of the transition, dν/ds denotes the frequency change along the photon path in the observer frame.If we denote the angle between the path and the radial direction θ and the velocity distribution v(r), the term of dν/ds is written as Since the velocity v(r) in the relevant region is smaller than a few thousand km/s, the terms of the order of O(v/c) 2 or higher are negligible.In actual calculations, (dν/ds) −1 is obtained from equation (A3) along a finite path δs the photon travels.
To simulate the trajectory of each photon packet, we consider the mean free path of each scattering/absorption process.The mean free path of photons for Hα transition, electron scattering, and free-free absorption are given as respectively, where κ ν is the free-free absorption coefficient, which can be obtained from the spectral emission coefficient, j ff,ν , and Kirchhoff's law j ff,ν = κ ν B ν .For l Hα , the line profile function is given by the δ-function in the Sobolev approximation where line broadening by the thermal motion of atoms is not taken into account.This can result in δν = 0, which is unphysical.To avoid that, δν is assumed to have a minimum value of 3 × 10 −5 ν Hα s −1 , which corresponds to the velocity of thermal motion for hydrogen atoms at T ∼ 10 4 K. Since the floor value of δν might change the depth of the absorption lines seen in the spectrum only if the velocity gradient is negative everywhere, we have to note what assumptions are used to determine the floor value.We will discuss this issue in section 3.1.For a single step of a photon packet, we determine which reaction occurs using a random number R in the range of 0 ≤ R ≤ 1 so that the probability of each reaction is inversely proportional to its mean free path.The actual distance l the packet travels is determined by another random number so that R = exp(−l/l sc − l/l ff − l/l Hα ).If a packet experiences Hα absorption, it is re-emitted at the same position with an energy of ν Hα in a random direction in the fluid comoving frame.When a packet is determined to be absorbed by a free-free transition of an ion, the packet is erased at the position.Our treatment of line radiative transfer is validated in Appendix B, with a test problem of modeling the P-Cyg line emerging from a homologously expanding medium (Jeffery & Branch 1990).

RESULTS
In this Section we present spectra calculated for models E and S introduced in section 2.1 at several epochs up to t = 80 days, where the time is defined as t = 0 at explosion.

Spectra of Hα line profiles
Figure 4 shows snapshots of spectra around the Hα wavelength for model E. The spectra generally show a combination of a narrow P-Cyg line and a broad emission line, with their shapes determined by the CSM structure immediate outside the shock wave.The CSM in the vicinity of the shock front is significantly accelerated due to the effect of radiation in the early phase (t ≲ 30 d in Fig. 3), while the CSM velocity in the vicinity of the shock front does not significantly change as the radiation flux decreases in the later phase (t ≳ 60 d).The wavelength of the absorption minimum in the P-Cyg profile likely traces the minimum velocity of the CSM affected by the radiation from the shock, as seen from the right panel of Figure 4.
The Hα emission line is broadened by a large number of electron scattering in the CSM, as previously claimed in Chugai (2001) and Huang & Chevalier (2018).The wavelength shift ∆λ when a line photon incident from a direction ⃗ n is scattered in a direction ⃗ n 1 is given as where λ Hα is the wavelength of Hα emission line in the rest frame of the hydrogen atom, ⃗ p denotes the momentum vector of an electron (Rybicki & Lightman 1979), and p = |⃗ p|.Here θ and θ 1 are the angles between the momentum vector of the electron and the vectors ⃗ n and ⃗ n 1 , respectively.Since photons travel a distance s = τ 2 es /(n e σ T ) on average in a CSM with an optical depth τ es (> 1) for electron scattering before diffusing out of the CSM, the line width can be estimated as The right hand side can be expressed as y/2 in terms of the Compton y parameter (Rybicki & Lightman 1979).
Comparison of equation ( 9) and the line widths in the left panel of Figure 4 indicates that Hα photons are emitted in a region with an optical depth τ es of at most a few, even in the early phases such as day 30 when the total Thomson optical depth of the unshocked CSM is ≈ 20.This is because photons emitted by Hα transitions in deeper regions are scattered by electrons and mostly change back close to the Hα frequency to be absorbed by the same Hα transitions, and hence cannot contribute to the observed Hα emission line.As the shock expands, the temperature and ionization degree in the CSM also decrease.Then the number of electron scatterings and the Compton y parameter decrease, resulting in a narrower emission line.
Figure 5 shows the time evolution of the spectrum for model S. At the initial phase, the line profile is composed of a combination of broad and narrow emission lines as in model E. The time evolution of the broad component is almost the same as model E, as this only depends on the time-varying τ es of the CSM that is similar between the two models.However, the shape of the narrow lines appear differently, with no clear absorption lines as seen in model E.
As seen in the right panel of Figure 5, a small absorption feature appears at the velocity of ∼ −30 km s −1 in the broad emission component after day 30.This absorption line is expected to appear from a medium with a negative velocity gradient, and the similar issue was already noted in Dessart et al. (2015).In model S the velocity gradient is negative everywhere, in contrast to model E in which the velocity increases with radius at the outer part of the dense CSM (the bottom left panel of Fig. 3), keeping the memory of the eruption.If the velocity gradient is negative, there is a direction for which the frequency shift (δν in Eq. ( 6)) vanishes.In terms of the angle θ from the radial direction, this condition can be expressed as because dν ∝ d(v cos θ) = (dv/ds cos θ − v sin θdθ/ds) ds = 0, dθ/ds = − sin θ/r, and dr/ds = cos θ along the path ds making an angle θ with the radial direction.This means that a photon emitted in this direction by the Hα transition is immediately absorbed and re-emitted in a different direction by the Hα transitions and results in a strong absorption feature in the spectrum.In reality, the thermal motion of hydrogen atoms prevents δν from vanishing and instead increases δν of the photon emitted in the direction satisfying equation ( 10) up to a value of the order of ν Hα 2k B T /m p c 2 at minimum.Thus we set the lower limit of 3 × 10 −5 ν Hα s −1 to δν, which corresponds to ν Hα 2k B T /m p c 2 at T ∼ 10 4 K, to mimic this situation as already mentioned in Section 2.2.2.The wavelength of the absorption feature traces the wind velocity as seen in Figure 5 and does not significantly change with time.This is in contrast with model E, in which the wavelength of the absorption feature moves to shorter wavelength with time.Figure 6 shows a comparison of the the spectra calculated for the two models at day 30.As mentioned in the previous section, while the broad lines are similar in the two models (left panel) the narrow lines appear very different (right panel).Model E shows a P-Cyg like profile with a narrow absorption line at v ≈ −100 km s −1 , while for model S no such deep absorption line appears.To understand what causes this difference, we overplot in Figure 7 the "initial" spectra of photon packets emitted at the start of calculation.Here, all emitted photons are sampled for the "initial" spectra, whereas the observed spectra (denoted as "final" in Figure 7) are composed only of photons escaping from the outer boundary.

Comparison of spectra for different CSM models
In both models, the "initial" spectra have a symmetrically broadened structure with respect to the wavelength of Hα in the rest frame.The Hα emission line is produced by recombination of hydrogen, and the broadening of the line up to ∼ 300 km s −1 indicates that Hα photons are emitted in the vicinity of the shock front where the CSM is accelerated by radiation from the shocked region.The width of the plateau in the "initial" spectra is determined by the minimum velocity in the ionized region, which is ∼ 80 km s −1 in model E and ∼ 30 km s −1 in model S. Most of these photons experience electron scattering, free-free absorption, and absorption and re-emission by the Hα transitions.The luminosity of the emergent emission line component decreases by a factor of ∼10 due to these processes, and the line is further broadened due to electron scattering as discussed above.The emergent spectrum has a blue-shifted peak in both models.This is because a significant fraction of photons entering the inner boundary of the computational region (i.e.shocked region), which are red-shifted in the observer frame, do not contribute to the emergent spectrum.The lack of these redder photons in the observed spectrum creates the blue-shifted peak.
As we describe below, the difference in the detailed structure of the emission line for the two models is due to the difference in the velocity gradient of the unshocked CSM (Figure 3).In model S, a narrow emission line appears at around −30 km s −1 , which is the bluest edge of the plateau in the "initial" spectra.This narrow emission line is composed of photons that are emitted in nearly radial directions, and escape from the CSM with no further reactions.In model S, the velocity gradient is negative at all locations, due to a combination of a flat initial velocity profile and the radiative acceleration being stronger at smaller radii.With such a velocity distribution, when a photon is initially emitted with ν Hα in the fluid comoving frame and travels in the radial direction, it would never have the frequency ν Hα again unless it experiences electron scattering.Since L λ of this emission line is reduced by a factor of ∼ 6 from the "initial" value, we infer that the emitting region of this line has an electron scattering optical depth of τ e ∼ ln(6) ≈ 1.8.This value is roughly consistent with τ es estimated by equation (9).
To further demonstrate this relationship between the velocity structure and the narrow line feature, we perform Monte-Carlo calculations of several models in which only the velocity distribution is artificially modified from model E while the density and temperature distributions are fixed (Fig. 8).For these modified models we assume that the CSM velocity from each location in the outside of the shock wave to ∼ 5 × 10 15 cm has a power law distribution v ∝ r α , with α > 0 (α < 0) corresponding to a positive (negative) velocity gradient as seen in the left panel of Figure 8.The right panel of Figure 8 shows the corresponding spectra for the assumed velocity distributions in the left panel.When α is positive, a larger α makes the absorption line deeper.On the other hand, when α takes a negative value,  the absorption line disappears and an emission line appears.These results confirm that the velocity gradient explains the difference in the line profiles in Figure 7.A similar trend of distinct absorption features in Hα line produced by CSM with different velocity distributions was pointed out by Chugai et al. (2004), which modeled a particular SN IIn 1994W.They claimed that a homologously expanding CSM produces prominent absorption lines, while what we show here is that CSM with a positive velocity gradient produces prominent absorption lines and the part of the CSM contributing to the Hα line formation does not expand homolgously in our model.

Blue-shifted excess component
Returning to the time-dependent spectra in Figures 4 and 5, an excess from the continuum emission is seen on the blue side at around 650 nm.The blue-shifted excess component emerges as a bump at around day 40, while from around day 60 a smaller and redder component appears between the main emission line and the first blue component.From inspection of the photon trajectory, we find that the bluer (redder) excess is formed by Hα photons from recombination in the unshocked CSM, which enters the shocked region and is scattered by electrons (absorbed and re-emitted by the Hα transition) there, and exits the shocked region towards the observer without further reprocessing.Figure 9 shows the spectrum and its components decomposed by the number of Hα absorption/re-emission in the shocked region for model E at day 60.The second redder component appears if the number of Hα absorption/re-emission in the shocked   region is larger than 1, while the component does not appear when the Hα process does not occur in the shocked region.Thus, we see that the second redder component is caused by photons that experience the Hα process in the shocked region.If a photon in the shock downstream encounters an electron moving with the downstream bulk flow and changes its direction closer towards the radial direction by scattering, the photon acquires energy from the electron in the observer frame as bulk Compton scattering process and results in the wavelength shift according to Equation ( 8).Such photons produce the first bluer excess, with a maximal wavelength shift of where v down is the downstream bulk velocity.If a photon in the downstream is instead absorbed and re-emitted by Hα transitions, the photon re-emitted in a nearly radial direction would also have a higher energy in the observer frame due to the bulk motion of the shock downstream.The resultant wavelength shift is at most, half of the wavelength shift due to electron scattering.This is why two excess components appear at the blue shifted wavelengths.
To quantitatively see the effect of shock crossing on the photon spectra, we show in Figure 10 the spectra at day 60, grouping the photons with the number of round-trips across the shock wave.The position of the blue-shifted excess seen in the total spectrum (e.g.left panel of Fig. 4) is mainly reproduced by photons with one round-trip.Photons experiencing two or more round trips make an even bluer and broader component, with the amount of blue-shift proportional to the number of the round trips.However the wavelength of the redder component hardly changes with the number of round-trips.This is because photons with multiple round-trips via electron scattering in the shocked region would gain energy each time they are scattered, whereas photons with round-trips via Hα absorption/re-emission would always have a frequency of ν Hα in the comoving frame upon re-emission.Since the number of photons significantly decrease for higher numbers of round trips, the component with one round trip mostly determines the total spectrum.Nevertheless, one may be able to see multiple blue-shifted excess components due to multiple round-trips, for a spectrum of a bright SN IIn with a good resolution and photon count.
We note that the second redder excess is sensitive to the temperature in the shocked region as it involves Hα transitions.In this particular model, the temperature of the shocked region between days 60 and 70 is T ≈ 3 × 10 4 K, just around the temperature where hydrogen could cause Hα transitions.If the shock wave were stronger and the temperature of the shocked region were higher, the second excess component would disappear.Thus while the second component may be sensitive to the assumed shock downstream temperature, the main bluer excess seen in our simulations should be robust.
The blue-shifted excess is evident only at late times, as the electron scattering optical depth in the unshocked CSM has to be sufficiently low for the excess to be notable.At early times, the large τ es spreads photons in this blue-shifted excess component to a wider wavelength range, which smooths out this component and makes it buried either in the broad line or the continuum.
Table 1 shows the velocity inside the shock wave and the corresponding minimum wavelengths for the photon with electron scattering or Hα absorption/re-emission inside the shock wave obtained with equations ( 11) and ( 12). Figure 11 shows each spectrum from day 40 to day 70 with vertical lines indicating the minimum wavelengths.Since these lines successfully trace the bluer edges of the two blue-shifted components, one can estimate the velocity of the shocked matter and its evolution, if either of the two components is clearly seen in the spectra.6

DISCUSSION AND CONCLUSION
In this study, we modelled spectra emitted from interaction of SN ejecta with an extended dense CSM, formed by either eruptive or steady mass-loss prior to the SN, using our newly developed Monte Carlo radiative transfer code.We focused on the Hα line, and investigated in detail the physical processes that cause the broad component and narrow P-Cyg profile, and the blue-shifted excess component which have been observed in some SNe IIn.
As summarized in the top panel of Figure 1, our model predicts three features in the spectra: the broad component, narrow component, and the blue-shifted excess component at late phases.The broad component is formed by Hα photons emitted by the recombination experiencing a large number of electron scattering in the CSM, consistent with the claims of Chugai (2001) and Huang & Chevalier (2018).The narrow P-Cyg profile is formed by Hα photons experiencing a small number of electron scattering and the subsequent Hα absorption and re-emission in the CSM.The blue-shifted excess is caused by Hα photons interacting with the matter in the shocked region, which then travel almost radially to the unshocked CSM.The excess can have two components, with interaction in the shocked region being electron scattering or Hα absorption/re-emission corresponding to a bluer or redder component, respectively.This excess does not appear until the optical depth for electron scattering in the unshocked CSM becomes small enough, which happens at around day 40 for our particular model.
Comparing the spectral shapes of the Hα line between the CSM models formed in the two scenarios, the narrow P-Cyg line is only observed for the CSM modelled by eruptive mass loss.We find that this is due to the difference in their velocity distributions of the CSM in the pre-SN phase.For CSM formed by steady mass-loss with constant velocity profile, the velocity after SN monotonically decreases with radius due to radiative acceleration.A photon initially emitted from the Hα transition cannot undergo Hα absorption unless the traveling direction satisfies Equation (10) or it changes its direction and wavelength via electron scattering.Otherwise the photon could never again have the wavelength of Hα in the comoving frame of any other hydrogen in its path.Therefore, the absorption line is absent in the steady mass-loss case, while it appears in the eruptive mass-loss case in which the velocity distribution includes a weak shock-like structure with positive velocity gradient.We confirm the relation of the absorption line feature and velocity gradient through test calculations using artificially modified velocity distributions and find that the absorption line appears when there is a positive velocity gradient, and that a steeper positive velocity gradient deepens the absorption line.
For both CSM models a blue-shifted excess appears in the spectra between days 40 and 70.This is caused by electron scattering and Hα absorption/re-emission in the shock downstream, corresponding to the first bluer excess and second redder one, respectively.From this finding, we analytically estimated the wavelength changes of photons by electron scattering or Hα transition using the flow velocity at the immediate downstream.The analytical estimates of the wavelengths at the bluer edges of the two blue-shifted excesses agree with the numerical results at all times.
From our relationship between the Hα line shape and the CSM structure, we claim that the CSM structure and its origin can be diagnosed by high-resolution spectroscopic observations of Hα lines.Since the presence of a narrow absorption line is important as evidence for the eruptive mass-loss, a resolution of λ/∆λ > c/v CSM ∼ 10 4 (v CSM /30 km s −1 ) −1 is required to distinguish the CSM-forming mechanisms.Here v CSM denotes the typical velocity of the unshocked CSM in the interaction region, and is dependent on the progenitor.While we considered here an RSG progenitor, the requirement on resolution would likely be less stringent for CSM formed from more compact stars, such as LBVs with expected v CSM ≳ 100 km s −1 .

Applications to Observed SNe
Our work is aimed to be an in-depth study of the various spectral features we expect from interaction-powered SNe, and we thus did not aim to fit a particular observed spectrum.Nevertheless, our findings can lead to some qualitative understanding of the CSM environment for some of the observed SNe.As a preliminary demonstration, we show a comparison between the observed spectrum of SN 1998S and our simulation results.We also apply our findings to two SNe IIn, SN 1997ab (Salamanca et al. 1998) and KISS-15s (Kokubo et al. 2019).
Since our models E and S give rising times and peak luminosities consistent with those observed in SN 1998S as shown in Figure 2, we also compare the spectra.We obtain the spectral data from WISeREP (Yaron & Gal-Yam 2012) and focus on the spectrum observed 9 days after the luminosity peak by Keck1/LRIS.Spectra from our simulations are constructed by binning photon packets to an equal bin-width of ∆λ = 0.8 nm and each bin is shifted by 0.1 nm for rebinning to be consistent with the observation.Figure 12 shows a comparison of the observed spectrum of SN 1998S and the simulation results 40 days after the SN explosion for both models E and S.This period is chosen because it roughly corresponds to 9 days after the luminosity peaks in models E and S. As shown in the top panel of this figure, the original simulation results in a too strong Hα emission line.On the other hand, the bottom-left panel shows the simulation spectra from models with the same CSM structures but with a reduced density distribution by a factor of 0.06, in which the density only in the unshocked region is reduced and that in the shocked region is unchanged, roughly reproduces the strength of the observed line.To create a less dense CSM at 40 days after the SN explosion without changing the shape of the light curve, there can be the following suggestion: shorten the time duration between eruption and SN explosion while maintaining the density of the inner dense part of the CSM that powers the light curve around peak, which creates a model where the shock radius is further outward and CSM density is lower at day 40.Since some different parameter sets result in similar light curves in CHIPS calculations, a detailed parameter survey is required to create a variety of flow-field profiles without changing the light curve, which will be clarified in the future work.
As shown in the bottom-right panel of Figure 12, the difference in line shapes between models E and S is not pronounced in the low-resolution spectrum while it is more significant in the high-resolution spectrum.Thus, high- resolution spectra are necessary to distinguish CSM structures formed from different mass-loss histories.The difference in the width of the narrow absorption line between the observation and simulations may be partly due to the limited wavelength resolution of the observation, which is indicated that the full width at half-maximum of the narrow line is less than 300 km/s (Leonard et al. 2000).Furthermore, it seems difficult to simultaneously explain the depth of the absorption line and the height of the emission line seen in the observed spectrum with a spherically symmetric CSM model, because deeper absorption lines require higher CSM density, which leads to too strong recombination emission to be compatible with the observation.To match both absorption and emission lines between the observation and simulations, a non-spherically symmetric CSM model with high density only around the line of sight of the observer is required.
As the bottom-left panel shows, the dip around 640 -650 nm in the observed spectrum is qualitatively reproduced by our simulations.As described in Section 3.3, the dip is located between the two photon-excess components created by electron scattering and Hα absorption/re-emission inside the shock wave, and the characteristic wavelengths of these components can be used to estimate the fluid velocity inside the shock wave.The observed dip is slightly bluer than the simulations, suggesting a faster shock velocity.A detailed parameter survey is still needed to produce a model with larger shock velocity without changing the light curve shape.On the other hand, it seems to be difficult for our model to reproduce the photon excess around 657 -660 nm in the observation, because this excess is thought to be composed of Hα photons that are emitted in the CSM in front of the shock, traveling toward the shock, and redirected to the observer by electron scattering or Hα absorption/re-emission with a positive line-of-sight velocity of several thousand km s −1 in the shocked CSM.However, due to the high optical depth inside the shock wave, such a photon would immediately be randomly oriented due to electron scattering.Therefore, a non-spherically symmetric CSM (as suggested by Leonard et al. 2000 from spectropolarimetry) is probably needed to reproduce this red-side photon excess.
SN 1997ab is one of the few SN IIn where high-resolution echelle spectra had been taken.Inspecting the spectra taken around 1 year after discovery (Salamanca et al. 1998), the co-existence of the P-Cyg line overlaid with the broad emission line suggests an eruptive origin for the dense CSM generating the lines.However, a notable difference from our resultant spectra is that the peak of the broad line in SN 1997ab is shifted to the blue side relative to that of the narrow P-Cyg emission line, whereas both of the peak positions are coincident in our simulations.Salamanca et al. (1998) explain that the red side of the broad emission line is reduced because Hα photons are self-absorbed in the red side resulting in the blue-shifted broad emission line.Based on our simulations, the presence of the narrow P-Cyg line suggests that eruptive mass-loss occurs.Time from the start of eruption to the present is written as t erup = t obs × (v sh /v CSM ), with the observation time of t obs ∼ 1 yr and v CSM ∼ 90 km/s, which is estimated from the width of narrow P-Cyg line in Salamanca et al. (1998).Although we cannot estimate the shock wave velocity, v sh , because the blue-shifted excess is not observed in SN1997ab, we obtain t erup ∼ several decades for v sh of several thousand km/s.The spectral features of this SN will be studied in future work, with a refined modelling of the CSM that reproduces the luminosity evolution.
KISS15s is a long-lasting SN IIn with a clear blue-shifted excess in the spectra from day 124.2 to day 431.4 (their Figure 20), which is also seen in our model at late phase.We compare the low-wavelength edge of the observed excess with the expected maximal wavelength shift from our model (Equation 11), and find v down ≈ 3500 km s −1 .From the Rankine-Hugoniot relation under the strong shock limit, this translates to a shock velocity of v sh ≈ 4100-4700 km s −1 for an adiabatic index of γ = 4/3-5/3.This is a factor of ≈ 2 larger than v sh ≈ 2000 km s −1 adopted by Kokubo et al. (2019) from the width of the intermediate-width lines.The change in v sh drastically affects the inferred mass-loss rate and total shocked CSM mass within a duration t duration , based on equations ( 19) and (20) of Kokubo et al. (2019), where L bol is the bolometric luminosity which is nearly constant (0.8-0.9×10 43 erg s −1 ) from day 100 to 600, v w ≲ 47.2 km s −1 is the CSM velocity constrained from high-resolution spectra, and ϵ is the radiation conversion efficiency.With other parameters fixed, the inferred Ṁ and M CSM,600d for our updated v sh are drastically reduced to Albeit the possible uncertainty in ϵ the inferred CSM mass is found to be more compatible with the envelopes of RSGs, and thus may be a more natural explanation given by the low v w that may favor an RSG over an LBV.The inferred mass-loss rate indicates a Thomson scattering optical depth in the unshocked CSM of at most where κ es is the Thomson scattering opacity.We note that this is an upper limit, under the assumption that the entire unshocked CSM is ionized.From equation ( 9), we expect that scattering in this CSM creates an emission line with half-width ∼ 1000 km s −1 , which is roughly consistent with the observed intermediate-width line.

Caveats and Avenues for Future Work
In this work we have done the Monte-Carlo simulations as post-process, assuming that the structure of the CSM does not change while the emitted Monte Carlo packets move in the computational region.We have tested the convergence of our results by changing the duration of sampling the packets by half shorter and an order longer, finding similar results.However, the post-process approach would not work when the hydrodynamical profiles of the shocked region and CSM significantly change during the diffusion time in the CSM, such as before the SN has reached its luminosity peak.Future work will consider time-dependence of the hydrodynamical profiles of the CSM and the shocked region, which may help when modelling the very early phase of the SNe uncovered by e.g.flash spectroscopy (Khazov et al. 2016;Yaron et al. 2017;Bruch et al. 2021).
The fractions of the ionization and excitation states are calculated assuming the local thermodynamic equilibrium with the gas temperature.This is a rough approximation in the CSM outside the photosphere.Thus we have discussed the information that can be deduced from the shape of the Hα line by an analytical method based on results of the Monte Carlo simulations.As a consequence, we find that the velocity structure of the CSM can be inferred from the line shape.
There are claims that a significant fraction of SNe IIn show evidence for asymmetric CSM (e.g., Leonard et al. 2000;Bilinski et al. 2018;Soumagnac et al. 2020).Our simulations using Monte-Carlo radiation transport can be straightforwardly extended to asymmetric CSM, by coupling to multi-dimensional hydrodynamical simulations of SNe interacting with aspherical CSM (Vlasis et al. 2016;McDowell et al. 2018;Suzuki et al. 2019).
We have discussed Hα spectra from a limited set of models, composed of a specific SN explosion with two CSM structures having a similar density distribution.It is known that SNe IIn have a diversity in their light curves and spectral evolution (see e.g., Taddia et al. 2013, for one such categorization), probably due to the diversity in the extent and mass of the CSM.While we adopted parameters of SN ejecta and CSM that lie close to the representative light curve properties of SN IIn, we did not attempt to fit our results with actual observations of SN IIn.In future work we plan to conduct a more comprehensive study on comparison of our model with observed SN IIn events, including the two events SN 1997ab and KISS-15s mentioned in this paper.
We thank Akihiro Suzuki for helpful comments throughout the advancement of this work.This work is supported by JSPS KAKENHI Grant Numbers 22K03688, 22K03671, 21J13957, 20H05639, 19H00693 MEXT, Japan, and by the Sasakawa Scientific Research Grant from The Japan Science Society.DT is supported by the Sherman Fairchild Postdoctoral Fellowship at the California Institute of Technology.Numerical computations were in part carried out on Cray XC50 at Center for Computational Astrophysics, National Astronomical Observatory of Japan.
Software: MESA (Paxton et al. 2011(Paxton et al. , 2013(Paxton et al. , 2015(Paxton et al. , 2018(Paxton et al. , 2019;;Jermyn et al. 2023), CHIPS (Takei et al. 2022(Takei et al. , 2023) )  The subscript 0 indicates that the quantities are defined in the fluid comoving frame, and σ Hα has a non-zero value when the frequency of a photon ν 0 in the fluid comoving frame is equal to ν Hα .Using this, the optical depth for Hα is as follows, τ Hα = σ Hα n Hα ds Using s 1 , we rewrite r 1 as r 1 = r p + s 1 n = [x p + s 1 n x , y p + s 1 n y , z p + s 1 n z ].
Therefore, cos θ 1 = r 1 • n r 1 = (x p + s 1 n x )n x + (y p + s 1 n y )n y + (z p + s 1 n z )n z r 1 .
Here, we interpolate the flow velocity between the values of v 1 and v 2 to obtain v p at the current photon position.
(ii) The case of cos θ p < − 1 − (r 2 /r p ) 2 The photon reaches the inner boundary of a cell by traveling a path s 2 , which corresponds to δs.Let r 2 = [x 2 , y 2 , z 2 ] be the position vector of the inner boundary where the photon reaches, we obtain from the cosine formula as in (i) Solving this for s 2 , we obtain δs = s 2 = −r p cos θ p − r 2 p cos 2 θ p − (r 2 p − r 2 2 ). (A5) Then we can rewrite r 2 as r 2 = r p + s 2 n = [x p + s 2 n x , y p + s 2 n y , z p + s 2 n z ].
Therefore, the cosine of the angle θ 2 between the photon's path and the radial direction (see Fig. 13) can be obtained as cos θ 2 = r 2 • n r 2 = (x p + s 2 n x )n x + (y p + s 2 n y )n y + (z p + s 2 n z )n z r 2 .
Thus we calculate the optical depth τ Hα for the Hα transition from equation (A4) by substituting δs and the frequencies in each case.

B. TEST CALCULATION OF THE P-CYG LINE
In order to verify our Monte Carlo calculations, we performed test calculations to reproduce the P-Cyg line from a spherically symmetric flow-field undergoing homologous expansion.The simulated Hα spectra are compared with those obtained from the semi-analytical model in Jeffery & Branch (1990).
The input of the Jeffery & Branch (1990) model is a homologously expanding medium, spanning from an arbitrarily defined "photosphere" to the outer edge, which reprocesses the continuum emission from the photosphere.Here, we choose the parameters of this medium to be similar to the CSM that are relevant to this paper.Specifically, the radius of the photosphere is set to 10 15 cm and that of the edge of the computational domain is set to 5 × 10 15 cm.The flow velocities at the photosphere and the outer edge are set to 20 km s −1 and 100 km s −1 , respectively.Photons are injected from the photosphere with a blackbody spectrum, and are initially assigned to a frequency range 0.999ν Hα ≦ ν 0 ≦ 1.002ν Hα in the fluid comoving frame.
The number density of hydrogen atoms that cause Hα absorption is assumed to follow a power-law profile of n Hα ∝ r −2 .The optical depth τ ph of Hα absorption at the photosphere is defined in the semi-analytic model as We perform simulations with τ ph = 10 and 100, under the assumption that each photon experiences at most only one Hα absorption/re-emission during its transport.To compute a single spectrum we use 10 7 photon packets, which are initially emitted from the photosphere and sampled at the outer edge of the computational domain.semi-ana τ ph =10 semi-ana τ ph =100 Figure 14.Comparison of P-Cyg lines between simulations and the semi-analytic calculations of Jeffery & Branch (1990).
Figure 14 compares P-Cyg lines calculated by the Monte-Carlo simulations and the semi-analytic model.For both cases of τ ph = 10 and 100, the simulated spectra are in good agreement with the semi-analytical results.Our radiative transfer code reproduces the P-Cyg line, as well as results of previous works (Roth & Kasen 2015;Wagle et al. 2023).

Figure 1 .
Figure 1.Schematic picture of the emission from interacting SNe that we model in this work.The upper side shows the line emission computed by Monte Carlo radiative transfer simulations (this work), and lower side shows the radiation hydrodynamical part computed by the CHIPS code.

Figure 3 .
Figure 3.The distributions of gas density (top-left), gas temperature (top-right), velocity (bottom-left), and the ratio of number densities of excited hydrogen that causes Hα absorption among total hydrogen (bottom-right) of the shocked region and the CSM at 10/12 (model E/S, respectively), 30, and 60 days after the the SN explosion.For model E at day 10, the unshocked CSM at r ≈ 5 × 10 14 cm is falling back with a negative velocity, hence does not appear in the bottom-left panel.

Figure 4 .
Figure 4. Time-dependent spectra in the eruptive mass-loss model (model E).For better visibility, the spectra after 20 days are multiplied by constants as indicated on the legend.The right panel is an enlarged view of the left panel to show the narrow lines around the Hα.

Figure 5 .
Figure 5. Same as Figure 4, but for the steady mass-loss model (model S).

Figure 6 .
Figure 6.Comparison of the Hα line shapes between models E and S at day 30.

Figure 7 .Figure 8 .
Figure 7.Comparison of the initial and final spectra at the day 30.Left:model E. Right: model S.

Figure 9 .
Figure 9. Spectrum and the components decomposed by the number of Hα absorption and the subsequent re-emission processes in the shocked region for model E at day 60.

Figure 10 .
Figure10.Spectrum and the components decomposed by the number of photon round-trips across the shock wave, for model E at day 60.Here a round-trip means that photons initially emitted in the unshocked CSM cross the shock wave, are scattered by electrons or absorbed and re-emitted by Hα transition in the shocked region, and return to the unshocked CSM.The "5 round trips" are photons that experience 5 or more round trips.

Figure 11 .
Figure 11.Spectra including the blue-shifted photon excess.The minimum wavelengths of the photon with electron scattering or Hα transition inside the shock wave are shown with the vertical lines.Top-left: the spectrum at 40 day.Top-right: the spectrum at 50 day.Bottom-left: the spectrum at 60 day.Bottom-right: the spectrum at 70 day.

Figure 12 .
Figure12.Comparison of the observed spectrum of SN 1998S and the simulation results 40 days after the SN explosion for both models E and S. Spectra from our simulations are constructed by binning photon packets to an equal bin-width of ∆λ = 0.8 nm and shifting each of the bins by 0.1 nm to be consistent with the observed spectrum.Top: Comparison with the simulation results of the original models of E and S. Bottom-left: Comparison with the models of decreased CSM density by a factor of 0.06.Bottom-right: Comparison of observation with simulated spectra with different bin widths.'r006' means the same as the bottom-left panel and 'r006HR' means high-resolution spectra with ∆λ = 0.1 nm.

APPENDIXA.
CALCULATION METHOD OF OPTICAL DEPTH FOR Hα ABSORPTIONHere we summarize the procedure to calculate the optical depth for the Hα transition using the Sobolev approximation.The cross section for Hα absorption is(Rybicki & Lightman 1979)σ Hα = πe 2 m e c f 23 ϕ(ν),where f 23 is the oscillator strength for Hα and ϕ(ν) is the line profile function of the photon frequency ν.In the Sobolev approximation, the line profile function is replaced with the delta function as follows,σ Hα = πe 2 m e c f 23 δ(ν 0 − ν Hα ).

Table 1 .
The fluid velocity inside the shock wave and the minimum wavelengths of the photon with electron scattering or Hα absorption/re-emission inside the shock wave.Time from Explosion (days) velocity [km/s] minimum λ for scattering [nm] minimum λ for Hα transition[nm]