{'fraction_non_alphanumeric': 0.034627639088267215, 'fraction_numerical': 0.021812961731031828, 'mean_word_length': 4.199144443433148, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 7, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 1, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
CHIPS: Complete History of Interaction-Powered Supernovae

We present the public release of the Complete History of Interaction-Powered Supernovae (CHIPS) code, suited to model a variety of transients that arise from interaction with a dense circumstellar medium (CSM). Contrary to existing modellings which mostly attach the CSM by hand, CHIPS self-consistently simulates both the creation of the CSM from mass eruption of massive stars prior to core-collapse, and the subsequent supernova light curve. We demonstrate the performance of CHIPS by presenting examples of the density profiles of the CSM and the light curves. We show that the gross light curve properties of putative interaction-powered transients, such as Type IIn supernovae, rapidly evolving transients and recently discovered fast blue optical transients, can be comprehensively explained with the output of CHIPS.


INTRODUCTION
Core-collapse supernovae (SNe) are explosive events which occur at the end of massive stars' lives. They show a variety of features in their spectra, mainly depending on the chemical abundances of the progenitor. Though emission lines are usually broad, some SNe show very narrow hydrogen emission lines in their spectra, which are classified as Type IIn (Schlegel 1990). The narrow width strongly suggests that the progenitor is surrounded by a slow-moving dense circumstellar medium (CSM). A dense CSM also efficiently converts the kinetic energy of the SN ejecta to radiation. Such "interactionpowered" events often have brighter light curves (LCs) than other types of SNe powered by the internal energy and/or radioactive decay of 56 Ni and 56 Co in the ejecta (Richardson et al. 2014).
Recent high cadence transient surveys have found that SNe IIn are not the entirety of interaction-powered transients. These surveys have discovered luminous transients that display rapid rise and decline (of within Corresponding author: Yuki Takei takei@resceu.s.u-tokyo.ac.jp 10 days) in the LC (e.g. Drout et al. 2014;Arcavi et al. 2016;Pursiainen et al. 2018), which are sometimes classed as rapidly evolving transients (RETs). Furthermore, in the past few years there have been a handful of discoveries of transients with even faster evolution on the order of days and luminosity exceeding 10 44 erg s −1 , named as fast blue optical transients (FBOTs; e.g., Prentice et al. 2018;Perley et al. 2019;Ho et al. 2020;Perley et al. 2021). While multiple interpretations exist for these transients, a fraction of these transients also show signs of the presence of dense CSM, that can potentially explain the brightness of these transients as well.
The key questions for these interaction-powered SNe are how the dense CSM has formed, and how the CSM is distributed around the SN progenitor when it explodes. In general, line emissions of metal induce steady mass loss with a rate up to ∼ 10 −3 M yr −1 (e.g., Castor et al. 1975;Vink et al. 2001). On the other hand, Kiewe et al. (2012) estimate from the strength of the Hα emission that the mass-loss rates of the progenitor of SNe IIn reach ∼ 0.1M yr −1 , which cannot be achieved by the line-driven mass loss. While theoretical models suggest some mechanisms to transport energy which induce arXiv:2109.05871v2 [astro-ph.HE] 24 Mar 2022 eruptive mass loss from massive stars (e.g., Woosley et al. 2007;Quataert & Shiode 2012;Moriya 2014;Smith & Arnett 2014;Woosley & Heger 2015;Soker & Gilkis 2017), there is no observational constraint on how they lose their envelope intensely. By contrast, there are recent analyses of SN observations which show that some massive stars experience eruptive mass loss before the actual SN event to form CSM, the famous example being SN 2009ip (e.g., Chugai et al. 2004;Pastorello et al. 2007Pastorello et al. , 2013Ofek et al. 2014a;Elias-Rosa et al. 2018).
Since interaction-powered SNe shine by converting the kinetic energy into thermal energy through interaction between SN ejecta and dense CSM, their LCs strongly depend on the density structure of CSM. Thus, the theoretical investigation into the LCs of interacting supernovae are direct probes for the surrounding CSM as well as the progenitor activity just before the core-collapse, which is too faint compared to an SN event to be observed.
Many works modelling LCs of interaction-powered SNe concentrate on the CSM formed by steady massloss ρ ∝ r −2 (where ρ, r denote the density and radius), or CSM which follows a single power-law (e.g., Moriya et al. 2011;Ginzburg & Balberg 2012;Chatzopoulos et al. 2012;Moriya et al. 2013;Dessart et al. 2015;Tsuna et al. 2019;Takei & Shigeyama 2020;Suzuki et al. 2020). However, it is unclear whether an eruptive mass loss event makes such a simple density structure around a star.
Recent numerical simulations indicate that the energy deposition at the stellar envelope likely triggers the mass eruption, although the precise energy transport mechanism is unknown (Dessart et al. 2010;Owocki et al. 2019;. In particular,  conducted one-dimensional radiation hydrodynamics simulations by injecting energy comparable to the binding energy of the envelope. They not only successfully reproduced the mass of the CSM inferred from Type IIn SNe, but also found that the density profile is quite different from a simple power-law. Hence, numerical simulations for eruptive mass loss is necessary in order to accurately predict the observations of interaction-powered SNe from LC modelling. For better modelling of interaction-powered SNe, we have developed the open source code "Complete History of Interaction-Powered Supernovae (CHIPS)" 1 that self-consistently calculates the LC of SNe interacting with the mass eruption from the progenitor. We show that the output of CHIPS comprehensively reproduces the LC properties of the observed interaction-powered transients. Here we introduce the code in detail and show some representative results. Should an interesting interaction-powered transient be discovered in the future, CHIPS would serve as a code for theoretical modelling of the LC and CSM.
This paper is constructed as follows. In Section 2 we briefly outline the calculations that are done in the CHIPS code, and what can be obtained as outputs of this code. In Section 3 we describe the technical details of the CHIPS code, focusing on the formulations and underlying assumptions of the mass eruption and LC calculations. In Section 4 we show results that we obtained from a series of calculations performed with the CHIPS code, and explore the variety of CSM density profiles and LCs obtained from our parameter study. We discuss the planned future updates of the code in Section 5, and conclude in Section 6. Appendix A summarizes how to install and execute CHIPS code and what kind of results are obtained.

OVERVIEW
We show the schematic picture of our code in Figure  1. The CHIPS code is composed of the following three components: 1. The code first loads a stellar model obtained by the MESA stellar evolution code (Paxton et al. 2011(Paxton et al. , 2013(Paxton et al. , 2015(Paxton et al. , 2018(Paxton et al. , 2019. The user can generate their own model using MESA, or use a sample stellar model at core-collapse generated by the authors using the example_make_pre_ccsn test suite in MESA version 12778. We have verified that the code is compatible with output files generated by the latest revision of MESA (r15140) at the time of writing.
2. Then using the input stellar model, the mass eruption and the resulting CSM are calculated using a radiation hydrodynamical code based on the code developed in , as detailed in Section 3.1.
3. Finally the bolometric and multi-wavelength LCs powered by the CSM interaction are numerically calculated using the methods detailed in Section 3.2. This calculation is based on the code developed in Takei & Shigeyama (2020).
For the first step, the stellar model can be specified with an appropriate path using the argument --stellar-model in the execution script run.py. Our sample consists of stars in a ZAMS mass range of Figure 1. Schematic overview and flowchart of the CHIPS code. The code first simulates the outburst just prior to core-collapse that creates a dense CSM. It then simulates the light curve from a supernova explosion powered by interaction with this CSM.
13 − 26M and with ZAMS metallicity of solar, and is used for the demonstration of our code in Section 4. For historical reasons the code for calculating the mass eruption is written in Fortran, and that for obtaining the LC is written in C. These codes can be compiled and executed with the standard gfortran and gcc compilers. The execution script run.py also internally does computationally inexpensive calculations to connect these steps. These include remeshing radial profiles of the star and CSM, generating opacity tables, and obtaining the gross property of the emission (e.g. peak luminosity, rise time) from the LC data.
Apart from various physical parameters that can be set in the stellar evolution done in MESA, our code consists of four parameters as shown in Table 1. We also show the range of values recommended by the authors to get a physical and reliable result. The detailed procedures to install and run the CHIPS code is shown in Appendix A.

Mass Eruption Calculation
The dynamics of the mass eruption and the formation of CSM are calculated by our 1-D Lagrangian radiation hydrodynamical simulation code developed by  and Kuriyama & Shigeyama (2021). Our code adopts the stellar model generated by the MESA code as the initial condition. The hydrogen-rich envelope is extracted from this file as the calculation target, and remeshed in order to enhance the resolution and to reduce the differences of the cell masses between neighboring cells prior to our hydrodynamical calculation. The total number of cells is specified by a parameter hydroNumMesh, with a default value of 10000.
After a short (50 seconds) initial relaxation of the remeshed initial models, thermal energy is injected into the inner boundary region, and the subsequent dynamical eruption is simulated. If the progenitor has a sufficiently massive (> 0.2M ) hydrogen envelope, the default mass coordinate of the inner boundary is set to be 0.2M outside its helium core to avoid a severe Courant condition caused by a sudden increase of the density around the boundary between the helium core and hydrogen envelope. Otherwise, a custom inner boundary must be set by the argument --eruption-innerMr (in units of solar masses) when executing run.py. This code is agnostic of the energy source, but the amount of injected energy and the duration of the injection are set by the parameters f inj and injectDuration respectively. The parameter f inj is defined as the ratio of the injected energy and the absolute value of the total energy inside Explosion energy of supernova ejecta 3 × 10 50 -10 52 ergs b a This parameter space reproduces the general characteristics (mass, radii) of the CSM observed in Type IIn supernovae . b Our assumption of homologous ejecta limits the explosion energy to sufficiently high energies (see also Section 3.2.1).
the envelope where r, v, ρ, e int , M r are respectively the distance from the center, velocity, density, specific internal energy and enclosed mass, G is the gravitational constant, and the three terms in the parenthesis are respectively the kinetic energy, internal energy and gravitational energy. Since the progenitor's envelope is gravitationally bound, E tot is negative for all progenitors.  finds that f inj of a few 10% results in partial ejection of the envelope and that the total ejected mass is sensitive to the value of f inj . The parameter injectDuration is by default set to 1000 seconds, which is much shorter than the dynamical timescale in the envelope where R * and M * are respectively the radius and mass of the progenitor. A longer timescale comparable to t dyn would make the gravitational pull from the central core important, thereby reducing the mass and energy of the erupted material .
A part of the envelope is erupted on the dynamical timescale and becomes a dense CSM. After enough time ( 2t dyn ) has passed since the mass eruption, the ejected envelope expands and radiation pressure gradient can no longer affect its dynamics. Therefore, at this stage the radiative transfer calculation scheme is turned off to save computational cost. In addition, because our interest is on the CSM profile at core collapse, the innermost region (< 10 13 cm) of the envelope is excluded from the computational range. The hydrodynamical simulation is continued until core-collapse, with duration given by the parameter t inj (see Table 1).
The density profile created after turning off radiative transfer is found to contain artificial shocks that propagate to the CSM (see green dashed line of Figure 2). This is because part of the bound CSM falls back from the star, creating shocks that do not disappear in the adiabatic approximation. Since the cooling and diffusion timescales in the shock-heated region are much shorter than the age, we expect that the radiative cooling will suppress these discontinuities from propagating to the CSM if correctly taking into account radiative transfer. In other words the discontinuity is artificial, and the profile should instead smoothly connect to the region unaffected by the shock. This is verified from long-term radiation hydrodynamical simulations done in Tsuna et al. (2021c) as seen in the red dotted line of the same figure.
where r * , ρ * , y are the fitting parameters. The code fits the outer region unaffected by the shock to the analytical profile, and extrapolates the fit inwards to r in = 2R * , where r in denotes the innermost radius of the CSM. This choice of the inner boundary is rather arbitrary; in our case it is mainly from the assumption of homologous supernova ejecta, which is valid only after the ejecta expands roughly twice its initial radius. Although r in is dependent upon M * , f inj , it does not affect the calculations of the CSM mass and thus the subsequent LCs (see next section). Outside the erupted material, the code adds by hand CSM composed of a stellar wind emitted by the progenitor before the eruptive mass loss. The CSM is assumed to have a steady wind profile ρ =Ṁ /4πr 2 υ w that extends out to 3 × 10 16 cm. The choice of the outer radius does not affect the final LC since the optical depth is only determined by the inner CSM. The mass loss ratė M and wind velocity υ w are assumed to scale with the progenitor's luminosity L * and effective temperature T * as (Nieuwenhuijzen & de Jager 1990; Mauron & Josselin The values L * and T * are taken from the stellar model generated by MESA. This prescription is purely for extending the computational region. Since the mass carried by the stellar wind is much smaller than the mass of the erupted material, the details of this do not affect the major characteristics of the final LC. A demonstration of this remeshing process is shown in Figure 2. The correction of artificial shocks in the inner region successfully reproduces the profile calculated from long-term radiation hydrodynamical simulations.

Light Curve Calculation
Several years to decades after the mass eruption, a supernova explosion occurs and the ejecta collide with the CSM. The shocked region formed from circumstellar interaction is enclosed with the reverse and forward shocks, which propagates into the SN ejecta and the unshocked CSM, respectively.
The radiative flux emerging from the forward shock heats the unshocked CSM. The CHIPS code calculates the bolometric and multi-band LCs following the methodology in Takei & Shigeyama (2020). The details of the LC model are mostly the same as Takei & Shigeyama (2020), but we have made several important updates for performance and compatibility with CHIPS. Thus we describe the detailed setup and the numerical procedures, as well as various updates from Takei & Shigeyama (2020).

Parameter setup
For the SN ejecta, we adopt homologous ejecta with a double power-law density profile (Matzner & McKee 1999) which is valid roughly after the ejecta expand to r ≈ 2R * and kinetic energy dominates over thermal energy. The constants g and υ t are functions of the ejecta mass M ej and energy E ej as The two exponents n and δ are constant. The explosion energy E ej is a free parameter in the code. The code calculates three light curves for the three different energies {1, 3, 10} × 10 51 erg by default. The ejecta mass M ej is calculated from where M cc is the mass of the star at core-collapse obtained from the MESA stellar model, M CSM is the mass of the CSM measured from r = 2R * , and M rem is the mass of the remnant compact object formed in the center. Since we cannot find the true radius of the progenitor due to the artificial shock as shown in Figure 2, we assume that the progenitor expands to r = 2R * at core-collapse 3 . Thus the mass enclosed between the true radius and r = 2R * is contained in the ejecta. From equation (3), the mass of the CSM is estimated to be M CSM r * 2R * 4πr 2 ρ CSM (r)dr ∝ r 1.5 * −(2R * ) 1.5 and very weakly depends on the choice of the inner boundary as long as 2R * r * . This is generally satisfied for t inj of a few years or longer. The value of M rem is estimated from the mass of the progenitor's carbon-oxygen core M CO (also obtained from the MESA model), using a simple formula obtained for single stars (Schneider et al. 2021) where x CO = M CO /M . The intermediate range is actually expected to form black holes as remnants, possibly with much weaker explosion energies as low as E ej = 10 46 -10 48 ergs and smaller ejecta mass (Lovegrove & Woosley 2013;Fernández et al. 2018;Tsuna et al. 2020;Ivanov & Fernández 2021). The code is nonetheless only capable with values of E ej around 10 51 ergs seen in canonical supernovae (see Table 1), because our assumption of homologous ejecta (υ = r/t) breaks down for such low-energy explosions. This is especially the case for red supergiants, which takes months for the ejecta to expand to twice the progenitor's radius.
The inner exponent δ typically takes a value between 0 and 1 (Matzner & McKee 1999), and here δ is fixed to be 1. The value of the outer exponent n is determined by first obtaining a "global" polytropic index N pol inside the envelope at core-collapse that best fits the relation p = Kρ 1+1/N pol , where constants K and N pol are the fitting parameters. Then we use a formula connecting N pol and n as (Matzner & McKee 1999;equation 25) where β ≈ 0.19 only weakly depends on N pol . The abundance is assumed to be uniform throughout the computational region, and the code takes the representative helium mass fraction Y as the average of that inside the CSM. We fix the metal mass fraction Z, and adopt the solar abundance of Z = 0.014 (Asplund et al. 2009). The hydrogen mass fraction is then obtained as

Shock Structure
The CHIPS code calculates the bolometric and multiband LCs at an observer in two steps. The code first calculates the shock structure as a function of radius r for each of the shocks, by the following equations at the shock rest frame that assumes steady state.
For simplicity we assume gas and radiation are in thermal equilibrium, and set the pressure p and specific internal energy e as a function of ρ and temperature T (= T g = T r ) as where µ, m u , k B , a, T g , T r are respectively the mean molecular weight, the atomic mass unit, the Boltzmann constant, the radiation constant, the gas temperature, and the radiation temperature. Though the radiative flux F was calculated by assuming a diffusion approximation in Takei & Shigeyama (2020), we refine this and adopt a flux-limited diffusion (Levermore & Pomraning 1981) to extend our model to the optically thin regime, where E = aT 4 , c is the speed of light, and κ R is the Rosseland mean opacity. For the Rosseland mean opacity κ R , we use the TOPS opacity table (Magee et al. 1995), which tabulates κ R as a function of density and temperature.
To solve the hydrodynamical equations (12)- (14), we first set boundary conditions at the immediate downstream of each shock by the jump conditions, which can be written in the rest frame of the shock as where the subscripts down (up) denote variables at the downstream (upstream) of the shock. The upstream density ρ up for the forward (reverse) shock is given from the density profiles in the CSM (ejecta), and p up , e up are assumed to be 0 for both shocks.
There are three unknown parameters in the equations: the velocities of the two shocks and the radiative flux at the forward shock. The radiative flux at the reverse shock front is set to 0 since the reverse shock weakly dissipates the kinetic energy into radiation compared to the forward shock. We derive these three values so that the velocity, pressure, and flux become continuous at the contact discontinuity by integrating the equations iteratively. The detailed procedure to calculate the temporal evolution of the shock waves is given in Takei & Shigeyama (2020).
One of the most important updates from Takei & Shigeyama (2020) is on the initial condition. In order to start the calculation, we need to set the initial position of the reverse shock r 0 and initial time t 0 . Takei & Shigeyama (2020) fixed t 0 = 1 day for simplicity. The assumption of a steady-state shock can be violated at early phases when the photon diffusion is slower than the expansion of the shocked shell, which we have seen to occur for a large f inj or small t inj . We will see the situation where the photon diffusion time becomes comparable to the expansion time scale of the shocked shell by using a simple model and explain a method to deal with this situation below.
As seen in equation (3), the inner density profile of the CSM closely follows ρ ∝ r −1.5 . Initially the shocked region is well approximated to be adiabatic with an index γ = 4/3, and a self-similar solution for the shock propagation in a CSM with ρ = qr −1.5 (Chevalier 1982) gives the radius and the velocity of the shock front as functions of time as where A is a constant dependent on n and γ. The value of q is determined by the total mass M CSM of the CSM distributed within a radius r out as The optical depth of the unshocked CSM is given as where the average opacityκ is assumed to be given by the Thomson opacity. The latter equation assumes that the profile ρ ∝ r −1.5 continues up to a certain radius r * > r sh , and that the contribution beyond r * is negligible 4 .
We further assume r * r sh so that the second term in the parenthesis can be neglected. Therefore the diffusion and dynamical timescales become comparable when t and r sh satisfy the following condition Substituting this into equation (22) we solve for r sh and obtain We plot r sh diff as a function of q in Figure 3 for three values of n that are realistic for stellar progenitors. We see that the steady-state assumption is invalid at r = 2R * for q 10 9 cgs. If the CSM expands for a time t inj with a velocity comparable to the escape velocity at the stellar surface υ w ∼ 2GM * /R * , then the CSM will extend out to r out ∼ υ w t inj . Thus the radiative equilibrium is not realized at around the forward shock in the CSM with the corresponding mass until the forward shock reaches the radius r sh, diff . We thus set r 0 to be the larger of 2R * and r sh, diff . Then following Takei & Shigeyama (2020) we adopt an analytical model of shock propagation (Moriya et al. 2013) to estimate the corresponding t 0 .
In the case of r sh, diff > 2R * , the shock breaks out at the radius r sh,diff inside the CSM. The internal energy E int, bo that has been stored as radiation in the shocked region before the shock breakout is released on roughly the dynamical time. To take this effect into account, we inject this internal energy as a source term at the inner boundary (r = r fs ) when solving the radiative transfer equation in the unshocked CSM outlined in the next section. The energy is injected depending on time in the form of an energy flux as (Arnett 1980 The internal energy stored in the shocked region at time t is proportional to the ejecta's kinetic energy dissipated up to time t, and the proportionality factor depends on the CSM density profile within the breakout radius, γ, and n. We fix the first two as ρ ∝ r −1.5 and γ = 4/3, and tabulate the value of this factor for various values of n using the self-similar solution (Chevalier 1982). For a given value of n we obtain the energy E int, bo through interpolation.

Radiative Transfer in Unshocked CSM
By solving equations for the shock structure at each time step, we obtain the flux F at the forward shock as a function of time from equation (17). To obtain the thermal structure inside the unshocked CSM with this flux and the LCs observed by a distant observer, the CHIPS code solves the following two-temperature radiative transfer and energy equations for E(= aT 4 r ), and U (= 3k B T g /(2µm u )), where U denotes the specific internal energy of the gas, p g the gas pressure, and κ a the absorption opacity. The boundary condition for the radiative flux is obtained as follows. Integrating equation (30) over the volume enclosed between r = r fs and r = r out yields rout r fs where the source term is omitted. From this equation, we obtain the time derivative of the total radiation energy in this volume as ∂ ∂t rout r fs 4πr 2 Edr = rout r fs ∂(4πr 2 E) ∂t dr − 4πr 2 fs E(r fs )u fs −4πr 2 out F (r out ) + 4πr 2 fs (F (r fs ) − E(r fs )u fs ), Since the last term F (r fs ) − E(r fs )u fs denotes the flux at the shock front, this should be equal to F at the forward shock (Otherwise the radiation energy would be lost or generated at the shock). Thus we obtain as the boundary condition.
Takei & Shigeyama (2020) calculated the absorption opacity κ a by subtracting the scattering opacity from the OPAL Rosseland mean opacity (OPAL opacities: Iglesias & Rogers 1996). However, the above opacity is not realistic as the absorption opacity for the source term in equations (30), (31) since the Rosseland mean opacity is derived assuming that the region is optically thick and used for the diffusion approximation. Therefore, the Planck mean opacity is newly implemented in CHIPS as the absorption opacity from the TOPS opacity table (Magee et al. 1995) in order to more precisely calculate the source term. We generated the Rosseland and Planck mean opacities for helium fractions of Y = 0, 0.1, · · · , 0.7 with a fixed metallicity of Z = Z . For other helium abundances, opacities are generated by linear interpolation.
We also implement multigroup radiative transfer calculation using post-process ray-tracing method to obtain multi-color LCs at U, B, V, R, I bands. For ray-tracing we follow the methods in Kutsuna & Shigeyama (2015) and Suzuki et al. (2021). We first integrate the radiative transfer equation along each ray passing through the CSM with a different impact parameter b, where I ν , S ν are the intensity and the source function. κ ν , σ ν are the absorption and scattering opacity at frequency ν. We can integrate this equation with respect to the optical depth dτ ν = (κ ν + σ ν )ρds and obtain the formal solution, where we assume that S ν is the Planck function B ν (T g ). Again, we use the TOPS opacity for the frequencydependent absorption opacity to keep the consistency with the Planck mean opacity. We obtain I ν (b, 0) from the flux emerging from the forward shock by assuming that the spectral shape is blackbody for simplicity if a ray penetrates the shocked region. On the other hand, along rays with b > r fs , the boundary value of I ν is 0 since the ray emanates from infinity. We can summarize the boundary conditions for I ν as where T fs satisfies F (r fs ) = (ac/4)T 4 fs . Then, equation (36) is integrated from τ ν = 0. We note that the exact spectral shape should depend on the density of the CSM. A less dense CSM , for instance, that of a stellar wind with a mass-loss rate of 10 −2 M yr −1 for a wind velocity of υ w = 100 km s −1 , is expected to produce a harder spectrum that extends to X-rays (Chevalier & Irwin 2012;Svirski et al. 2012;Tsuna et al. 2021a).
We calculate the luminosity per unit frequency, L ν by integrating the surface intensity I ν (b) = I ν (b, τ ν (b)) at each impact parameter b over the surface as, where b max = r out . Finally we can calculate the absolute magnitude at each band, using a filter function in Bessell (1990).

CHIPS RESULTS
Here we present some results from a series of simulations performed with the CHIPS code and explore the dependence of LCs of interaction-powered SNe on four key parameters: the mass M * of the progenitor, the injected energy f inj for the precursor mass eruption and its timing t inj , and the explosion energy E ej of the SN. First, we present and discuss results of mass eruptions with various parameter sets. Then we present resultant LCs to summarize the dependence of some observable quantities on the key parameters, which will be useful to discuss the origin of an observed interaction-powered transient. We compare CHIPS results with other works in Appendix B.
In what follows, we present our results for massive stars evolved with the example_make_pre_ccsn test suite in MESA version 12778. To reduce the dimensionality of our parameter space, we have restricted the metallicity to be solar and neglected stellar rotation. The properties of the progenitors that we have created are summarized in Figure 5.

Mass Eruption: Dependence on the Progenitor
and Injected Energy Figure 6 shows the dependence of erupted mass on the injected energy f inj . As shown in this figure, the erupted mass increases with f inj for all the models we create. The erupted mass of model M * = 19M is lower than those of models M * = 17, 18M . The photospheric radius plotted in Figure 5 is large at M * = 19M , which means that the shock breakout occurs at more inner part of the envelope compared to other progenitor models. Due to the inefficient energy transport towards the outer envelope, the erupted mass becomes small. This tendency is prominent for lower f inj .
The density profiles of the remeshed CSM for some parameters f inj , t inj are plotted in Figure 7. While the density profile at inner region well follows ρ ∝ r −1.5 at t inj = 10, 30 yrs, the exponent slightly deviates from −1.5 at t inj = 3 yrs. This is because the gravity of the progenitor does not control the motion of the inner region yet or the matter does not start to fall back yet. It takes significantly longer time than the free-fall time scale given by ∼ r 3 /GM * ∼ 0.9 yr(r/10 14 cm) 3/2 (M * /10 M ) −1/2 for the inner region to enter this fall-back phase (Tsuna et al. 2021c).

Supernova Light Curves
Bolometric LCs are a useful probe of SNe because i) bolometric LCs can be obtained from a simplified numerical model and ii) the bolometric LC constructed from photometric observations include important information about SN events such as the heating source, the size of the progenitor, the explosion energy, the CSM structure, etc. The CHIPS code provides bolometric LCs for interaction-powered SNe, as shown for several example sets of parameters in Figure 8. We focus on two observables that are often used to characterize the LC: Figure 5. Properties of the progenitors we adopt in this work. The left panel shows the masses of the hydrogen-rich envelope, helium core and the carbon-oxygen core, and the remnant mass Mrem obtained from equation (10). The right panel shows photospheric radius (in R ) of each star and the effective temperature (in Kelvins). These progenitors are available when downloading the CHIPS code, and one can use them without running the MESA calculation. Figure 6. The erupted mass as a function of finj at tinj = 10 yrs. The CSM mass decreases with increasing tinj due to the fallback of the bound mass elements towards the progenitor.
peak luminosity L peak and rise time t rise . In this subsection we summarize the dependence of these observables on the input parameters of CHIPS. Figure 9 shows the dependence of L peak and t rise on t inj . We have defined t rise to be the elapsed time from the luminosity being 0.01L peak to L peak . As seen from the figure, L peak becomes larger for shorter t inj , while the dependence of t rise on t inj is more complex. The rise time depends on the transition radius r * and the shock breakout radius R d ≡ (2κqυ fs /c) 2 , which is equal to r sh,diff in the limit of large n. For convenience we define a ratio of these two radii  where we have assumed for simplicity that the CSM extends as a power law ρ ∝ r −1.5 and truncates at r = r * . Following the calculations of Chevalier & Irwin (2011), we can analytically obtain the rise time for the profile ρ ∝ r −1.5 as

Dependence on the Timing of Eruption tinj
As these two are discontinuous, we instead adopt their harmonic meañ which peaks at R d /r * ≈ 1 when r * is fixed. To see if the numerical results from CHIPS match this expectation, we test two cases with different t inj of 10 years and 3 years, and the other parameters fixed as M * = 15M , f inj = 0.8, E ej = 10 51 erg. For t inj = 10 years, we obtain from the analytical model R d /r * ≈ 0.21 < 1 andt rise ≈ 16 days. The numerical results from CHIPS show a good agreement with this, with t rise ≈ 16 days. On the other hand, changing t inj to 3 years gives R d /r * = 80 1 andt rise ≈ 8.6 days. The numerical simulation nicely matches as well in this case, with t rise ≈ 9.0 days.
The analytical model in equation (40) implies that for an optically thick CSM the rise time peaks at a certain t inj for fixed f inj , M * , E ej . We can derive the relation ρ * ∝ t −3 inj , r * ∝ t inj since around r ∼ r * the CSM expands homologously 5 . Substituting it into equation (39), we get R d /r * ∝ t −4 inj υ 2 fs , which is sensitive to t inj . For the same parameter set M * = 15M , f inj = 0.8, E ej = 10 51 erg, we find thatt rise reaches a maximum of ∼ 30 days at t inj ∼ 7 years.
From Figure 8, the LC is found to decline more slowly for larger t inj . This is simply because the dense part of the CSM can be more extended for larger t inj , and the SN ejecta takes longer time to sweep up the entire CSM. The duration of this phase is roughly which is clearly seen as a plateau with a cutoff in the case where r * /R d 1 is achieved, e.g. for the t inj = 10 yr and 30 yr models. Otherwise the LC instead shows a monotonic decrease because the shock breakout component, with time dependence in equation (29), would be dominant throughout the time the forward shock runs through the dense part of the CSM.

Dependence on the Injected Energy finj
In Figure 10, we plot the t rise -L peak relation for f inj . The interpretation is rather simple compared to section 4.2.1. More mass is erupted from larger injected energy f inj . Except for the case of M = 19M , t inj = 3 yrs, E ej = 10 51 erg, L peak for all models plotted in this figure monotonically increase with f inj . This is roughly consistent with the previous result that the ratio of the radiated energy to the explosion energy is The rise time t rise also monotonically increases with increasing f inj because the erupted mass increases with f inj and forms a CSM with a large optical depth.

Dependence on the Explosion Energy Eej
The t rise -L peak relations as a function of E ej are plotted in Figure 11. Since υ fs ∝ E (n−3)/2(n−1.5) ej for given f inj , t inj , M * (Moriya et al. 2013), it takes a shorter time to reach the peak luminosity for a larger explosion energy. The larger E ej gives higher peak to the LC. It is expected from these results that the lower explosion energy gives lower peak and longer rise time.  Since the t rise -L peak relation does not largely depend on M * except for M * = 26M as seen in Figures 9, 10, 11, we cannot solve the degeneracy only by the relation. The difference of the M * = 26M model is caused by the difference of the ejecta mass. The calculated M ej for the M * = 26M model is ≈ 5.4M , which is much smaller than M ej ≈ 14M for the M * = 17M model. This difference arises from M rem , given in equation (10) and plotted in Figure 5. The carbon-oxygen core mass of the 26M model exceeds the threshold of black hole formation, corresponding to the second case of equation (10) where the value of M rem jumps. Since the ejecta with smaller mass moves fast for fixed E ej (υ ej ∝ M −1/2 ej ), the LC evolves rapidly and the peak luminosity becomes higher.

Multi-band Light Curves
We plot in Figure 12 absolute magnitudes in U, B, V, R, I bands. As seen from absolute magni-tudes of the model shown in the left panel, U > B > V R I at around their peaks. Since the color V − R ∼ 0.1 around the peak, the temperature of the model is estimated to be ∼ 4, 500 K. Changing E ej from 10 51 erg to 10 52 erg with fixed other parameters, multiband LCs are plotted in the right panel. Compared to the model with lower explosion energy, the color becomes bluer (corresponding to the color temperature of ∼ 14, 000 K.). This is caused by the higher radiative flux than that of the lower explosion energy model.
Multi-band LC modelling can be a tool to additionally constraint the parameters jointly with the t rise -L peak relation.

Comparing with Observations
We compare CHIPS results with photometric data of SN IIn 1998S. This SN was discovered on 1998 March 2.68 UT in NGC 3877 by Z. Wan (Li et al. 1998). The nickel mass is estimated to be 0.15 ± 0.05M from the Figure 11. The same as Figure 9, but the dependence on Eej. tail of the bolometric light curve (Fassia et al. 2000). The left panel of Figure 13 shows the bolometric LC of SN 1998S constructed by Fassia et al. (2000) and the CHIPS result for parameters E ej = 2.5 × 10 51 erg, f inj = 0.7, t inj = 11 yrs, M * = 20M . While we can successfully reproduce the bolometric LC for the first ∼ 60 days, the CHIPS model LC is fainter than the observed LC in the later epochs. This is because the energy input due to the radioactive decays of 56 Ni and 56 Co is not implemented in the CHIPS code (see also section 5.2). The large explosion energy of E ej = 2.5×10 51 erg, which is necessary to reproduce the observed LC, is consistent with the production of the large nickel mass of 0.15M (e.g., Hamuy 2003;Müller et al. 2017). The multi-band LCs obtained by Fassia et al. (2000); Liu et al. (2000) are plotted in the right panel, together with those of the same CHIPS model. As can be seen from the panel, each light curve in the early-middle epochs can be reproduced by CHIPS code as well.
In the following, we discuss the dependence of the characteristic observable quantities L peak , t rise on our results.
• FBOTs: L peak ∼ (1 − a few) × 10 44 erg/s, t rise ∼ (1 − a few) days (Arcavi et al. 2016) We can successfully obtain parameter spaces that are consistent the above observed regions in Figures 9, 10, and 11. The parameter spaces that give t inj , L peak of SNe IIn can be inferred as f inj > 0.8, E ej < 10 51 erg, t inj ∼ 10 yrs, for example. The indicated lower explosion energies are consistent with the small amount of 56 Ni compared to that of typical SNe (e.g., Smith 2014;Elias-Rosa et al. 2018). t inj has to be selected so that t rise takes a maximum for given other parameters. Meanwhile, there is a degeneracy in parameters which give a shorter rise time and higher peak luminosity, typical observables of RETs and FBOTs. Fitting the multi-color LCs to these observed quantities may break this degeneracy, which we will explore in a forthcoming study.

FUTURE IMPROVEMENTS
While we have shown that the current version of the CHIPS code successfully reproduces the diverse parameter spaces observed in interaction-powered SNe, there remain various simplifications and limitations. In this section we list our planned future improvements to the CHIPS code. We plan to release a new version of CHIPS when any of these updates have been done.

Relaxing the Steady-state Assumption in the Shocked Region
We resolve the shocked region at the shock rest frame with the assumption that the region is in a steady state. However, with this assumption we cannot take into account the effect of the diffusion of photons in the shocked region at early epoch, which may affect the shock structure itself. Moreover, we cannot resolve the shocked region at early times for some models in which the initial radius r 0 becomes larger than 2R * .
We plan to update our code by following the temporal evolution of the shocked region with time-dependent radiative transfer equation from the very early times. Then we can relax the approximations made when estimating the luminosity around CSM breakout (see Section 3.2), and we can predict the rise time and the peak luminosity more precisely.

Inclusion of Emission from the Ejecta
In the current version of the code we have taken into account only the emission powered by circumstellar interaction. While the interaction component is dominant around the peak of the LC, energy sources from the SN ejecta may become important when the CSM density drops at the late phase. In the future we plan to incorporate two robust energy sources into our LC model: internal energy deposited by the blast wave in the ejecta and radioactive decays of 56 Ni and 56 Co.
Though most of the internal energy in the ejecta is converted into kinetic energy as a star expands after the blast wave passes the stellar surface, a part of the internal energy is emitted from the expanding ejecta. The duration of this emission is determined by the expansion timescale, thus depends on the stellar radius. SNe originating from red supergiants shine for a few months, while those from Wolf-Rayet stars shine only for a few minutes by this energy source. If these SNe originate from red supergiants, then this source may contribute to the luminosity as much as another heating source, i.e., the radioactive decays of 56 Ni and 56 Co in the ejecta that emit energy with rates (Arnett 1982) Ni ≈ 4.8 × 10 10 erg g −1 s −1 exp[−t/8.8 days], (43) Co ≈ 2.6 × 10 8 erg g −1 s −1 exp[−t/110 days], (44) mainly in the form of γ-ray. The energy from these two sources, once stored in the optically thick ejecta, is gradually released into space as hydrogen ions recombine from the outside to the inside of the ejecta in Type II-P SNe. This displays a plateau in the LC with luminosity L p lasting for t p , given by the following formulae (Popov 1993) respectively, where T ion is the recombination temperature of hydrogen. The ejecta are subject to engulfment and photoionization from the reverse shock, which can modify the resulting emission. Nonetheless this component more or less exists in SNe IIn as well, and a simple comparison implies it could be important for many of our LCs in the late phase. In fact SNe with co-existence of narrow hydrogen lines and plateau LC are observed (Kankare et al. 2012;Mauerhan et al. 2013), and sometimes classified as SN IIn-P (Smith 2014). We note that the plateau would not affect our argument on the rise time and peak, since the former is generally much shorter than t p .
There are a handful of inferences and upper limits on the synthesized nickel mass from late phase observations of SNe IIn. Most have values around M Ni 0.04M (e.g. Smith 2014;Elias-Rosa et al. 2018), while brighter SNe may have larger values (e.g. SN 1998S in the previous section having M Ni ≈ 0.15M ). The input from radioactive decay at the early phase is expected to be 10 48 -10 49 erg. Therefore unlike Type I SNe, the contribution from radioactive decay is not expected to be dominant in the early phase.

Extension to Black Hole Formation
As mentioned in Section 3, not all the stars are expected to explode as canonical SNe with energy around 10 51 erg. There are multiple tentative evidence from observations that a non-negligible fraction of massive stars die without being visible as canonical SNe and form BHs (e.g. Smartt et al. 2009;Horiuchi et al. 2011;Kochanek 2014).
However the mass ejection of these kinds of corecollapse is poorly known, and can have a diversity. A weak "explosion" by reduction of gravity at the core due to neutrino emission (Nadezhin 1980) has been studied by recent simulations, which found explosions of energy 10 46 -10 48 ergs (Lovegrove & Woosley 2013;Fernández et al. 2018;Tsuna et al. 2020;Ivanov & Fernández 2021). This can be significantly modified if an accretion disk can form around the nascent black hole, which are predicted to launch outflows and/or jets that can carry much more energy (Woosley 1993;MacFadyen & Woosley 1999;Kashiyama & Quataert 2015;Quataert et al. 2019;Tsuna et al. 2021b).
The mass eruption studied here can naturally occur for massive stars that form black holes as well. There is a recent claim that circumstellar interaction after black hole formation can explain some peculiar transients found in recent optical surveys (Tsuna et al. 2020). Modelling these kinds of emission would be an important extension to our code.

CONCLUSION
We have developed the open source code CHIPS for modelling the LC of interaction-powered transients including the simulation of the mass eruption prior to a SN event. CHIPS successfully calculates the LC for four key parameters, the ZAMS mass of the progenitor M * , the injected energy normalized by the envelope's binding energy f inj , the time from energy injection to corecollapse t inj , and the explosion energy E ej . After showing the methodologies of CHIPS, we explored the dependence of the resultant CSM and LCs on these parameters, f inj = 0.3, 0.4, 0.5, 0.6, 0.8, t inj = 3, 10, 30 yrs, E ej = 10 51 , 3 × 10 51 , 10 52 erg, M * = 13, 14, · · · , 26M . It is found that the peak luminosity L peak becomes high for shorter t inj and larger E ej . Additionally, the rise time of a LC, t rise , is longer for larger f inj , E ej . In contrast to the simple relation of L peak with t inj , E ej , the dependence of t rise on t inj is complicated; a larger t inj can make the extent of the CSM larger but also reduces the diffusion time (see Section 4.2.1 for details). For an optically thick CSM, there is a maximum value of t rise when t inj is varied, for given f inj , E ej , M * . We find that the parameter space we explore using CHIPS covers a wide range of t rise -L peak space, including that of SNe IIn, RETs, and FBOTs. We expect that CHIPS will serve as the tool for observers to utilize in obtaining physical quantities of interaction-powered transients.
CHIPS will be updated and released when we achieve any of the planned updates: i) to follow the temporal evolution of the shocked region between SN ejecta and CSM including radiation transport, ii) to solve radiative transfer equation in the unshocked ejecta in order to reproduce the LC at late phase, iii) to extend to the cases of black hole formation, or iv) any other improvements in the modelling that would be beneficial to the transient community.
We deeply thank the anonymous referee for helpful comments that greatly improved the manuscript, and Viktoriya Morozova for helpful guidance on using the SNEC code. Y.T. is supported by the RIKEN Junior Research Associate Program. D.T. is supported by the Advanced Leading Graduate Course for Photon Science (ALPS) at the University of Tokyo, and by the JSPS Overseas Challenge Program for Young Researchers. This work is also supported by JSPS KAKENHI Grant Numbers 21J13957, JP19J21578, JP20H05639, MEXT, Japan.
Software: MESA (v12778; Paxton et al. 2011Paxton et al. , 2013Paxton et al. , 2015Paxton et al. , 2018Paxton et al. , 2019 Here we present the procedures to run the CHIPS code and describe the output files generated by the code. The name of output files and physical quantities listed in each file are shown in Table 2, and Figure 14 explains the directory tree of CHIPS.
1. First check whether gcc, gfortran, python, numpy and scipy can be called/imported on the command line. If any of these are not installed, install them.
2. Then install mesa reader, a module that lets us easily extract data from the mesa output. This is used in some of our Python scripts. The detailed description of the module is in http://mesa.sourceforge.net/output.html.
3. Clone the CHIPS repository using the command git clone https :// github . com / DTsuna / CHIPS . git 4. The sample pre-SN models generated by the authors are compressed in a zip file in the directory input/mesa_models. If one plans to use them, unzip this zip file.
5. At the top directory, compile the scripts for the mass eruption and light curve (LC) calculations using the makefile, with the command make.
6. Execute the script run.py with the parameters in Table 1 as arguments. For example, to simulate a star with ZAMS mass 15 Msun and solar metallicity (assumed to be 0.014) which experiences mass eruption with parameters f inj = 0.3 and t inj = 5yr, and finally explodes with energy E ej = 10 51 erg, run the following command python run . py --tinj 5 --finj 0.3 --Eej 1 e51 --analytical -CSM --stellar -model input / mesa_models /15 Msun_Z0 .014 _preccsn . data The argument --stellar-model specifies the MESA stellar model file to be used as input for the mass eruption calculation. The option --analytical-CSM corrects the artificial shock compressions that arise from the mass eruption code before obtaining the LCs. For details of this procedure, see Section 3.1 and Figure 2.
The mass eruption calculation reads in the MESA stellar model file given by the argument --stellar-model, and outputs files with the CSM profile under the directory EruptionFiles/. The calculation also provides profiles at the intermediate stages, outputting one file per year with the name intermediateXXyr.txt, where XX are integer coded as two digit from 01 to t inj . Another set of files resultYY.txt (YY being integers from 01 to 99), which finely resolves the period before and around the mass ejection event, is also created. In these files, the radius, velocity, density, hydrogen and helium fraction, and pressure of CSM XX years after the energy injection are listed as functions of the enclosed mass. The final CSM profile at t inj is outputted with the name atCCSN.txt, and is used to calculate the LCs.
The LCs of interaction-powered SNe are then simulated, with a set of explosion energies given by the argument --Eej. This argument can be given multiple times to calculate LCs for multiple values of E ej , e.g., --Eej 1e51 --Eej 3e51. When nothing is given, E ej = {1, 3, 10} × 10 51 erg are used as default values. The LC calculation yields files named shock_output_...erg.txt and IIn_lightcurve_...erg.txt as outputs under the directory LCFiles/. In LCFiles/shock_output_...erg.txt, the radii of the reverse and forward shocks, the velocities of the both shocks, the radiative flux at the forward shock are listed as functions of time since explosion, while LCFiles/IIn_lightcurve_...erg.txt list the luminosity at the edge of CSM r = r out = 3×10 16 cm, the photospheric radius, and the color temperature. Multi-band LCs (see examples in Figure 12) can also be requested by the argument --calc-multiband; in this case, another set of files called LCFiles/IIn_lightcurve_...erg_mag.txt are created. These files contain LCs in U, B, V, R and I bands, with 1 day interval.
As done in the above example command, one can use the sample MESA models generated by the authors for the argument --stellar-model. We have 11 models for stars of a solar metallicity in the ZAMS mass range 13 − 26M , with 8 models in 1M increments for masses up to 20M and 3 models in 2M increments for masses greater than 20 M . The pre-SN models are in a zip file in the directory input/mesa models/. After un-zipping this file, you will find MESA data files with the naming showing the mass and metallicity at ZAMS.
If a mass eruption calculation has already been done and the CSM profiles EruptionFiles/intermediateXXyrs.txt exist, one can use those files as CSM to obtain the LCs of interaction-powered SNe. This is done with the code after_eruption.py, with a command like python after_eruption . py --Eej 1 e51 --analytical -CSM --stellar -model input / mesa_models /15 Msun_Z0 .014 _preccsn . data --profile -at -cc EruptionFiles / intermediate05yr . txt By using this code one can avoid re-doing the mass eruption calculation, which is computationally costly compared to the LC calculation.

B. COMPARING WITH OTHER WORKS
We compare our CHIPS code with past works on similar LC modelling, and adopt SNEC as an example as it is one of the most widely used open-source codes. SNEC 6 is an open-source 1-D Lagrangian code for the hydrodynamics and equilibrium-diffusion radiation transport, aiming at modeling LCs of core-collapse SNe (Morozova et al. 2015). More recently SNEC has also been used to simulate LCs of Type II-P/L SNe interacting with dense CSM (e.g., Morozova et al. 2017;Das & Ray 2017;Morozova et al. 2020;Li & Morozova 2021).
We simulate an explosion and the resulting LC of our 15 M progenitor using SNEC, using the CSM corresponding to a mass eruption with parameters f inj = 0.3, t inj = 5 yrs. We first stitch the profile of the inner helium core (obtained by MESA) to that of the outer hydrogen envelope obtained from the hydrodynamical simulation in the CHIPS code. Then we remove the innermost ≈ 1.4 M from the computational region, corresponding to the value of M rem in equation (10) for our 15 M progenitor. To simulate the explosion, SNEC injects energy by hand at the innermost computational region. Using the thermal bomb formalism, we set the value of the injected energy so that the final explosion energy becomes E ej = 10 51 erg.
The red solid line in Figure 15 shows the resulting LC obtained by extracting the luminosity at the outermost cell, which is always optically thin and travels much slower than the speed of light. The LC consists of two main phases: CSM breakout and subsequent CSM interaction (t 10 days), and the plateau phase powered by cooling and radioactive decay of 56 Ni in the ejecta. In the same figure we also show in blue dashed lines the LC obtained by CHIPS code, which at present only simulates the contribution from the former CSM interaction. While the timescale of the emission is roughly in agreement, the peak luminosity in our calculation is dimmer by a factor of few. We presume this to be due to our simplified assumption of the ejecta density profile as a double power-law (equation 6), where the outer ejecta follow n ≈ 11 for this parameter set (equation 11). This simplified profile is often used in the literature, and in our case is required for our semi-analytical evaluation of the shock breakout emission in and after equation (27).
The actual ejecta density profile is predicted to have a smooth transition from the inner and outer asymptotic power-law indices (Matzner & McKee 1999, Figure 7). For interaction-powered SNe with massive CSM, the outer unshocked ejecta at CSM breakout will have a flatter profile than the asymptotic profile at the highest velocity. Thus a flatter density profile of n = 7 for the outer ejecta was assumed in some of the previous works (Chevalier & Irwin 2011Svirski et al. 2012). This is indeed seen from inspection of the density profile of the unshocked ejecta in the SNEC calculations. With this in mind, we additionally calculate another model where the exponent of the outer density profile is set by hand to be n = 7, with all the other parameters fixed. We show this LC as green dash-dotted lines in Figure 15. We find that a model with n = 7 agrees well with the LC from SNEC. mIn essence, assumptions in the density profile of the outer ejecta can somewhat affect the LC, and the simple double power-law profile should eventually be updated to a more realistic one including curvature. Since the only limitation is from the semi-analytic framework of the breakout emission (which requires a power-law ejecta profile for self-similarity), this can be realized once the early-phase LC modelling is replaced to a more rigorous radiation transfer simulation, as outlined in Section 5.1.
We note that past works using SNEC have mainly focused on modelling multi-color/bolometric LCs of Type IIP/L SNe interacting with a relatively compact CSM, that formed by mass eruptions which occur within months to a few years before the terminal explosion. Although this situation is similar to our model, at present it is difficult to conduct our calculations under their setting. This is because at such short timescales (i) the CSM extends only to several times the progenitor radius, so the ejecta is not homologous during the interaction phase, (ii) the LC code suffers numerical instabilities at high densities during the interaction, and possibly (iii) the CSM profile deviates from the analytical profile assumed throughout this work. These will be addressed in future updates of CHIPS to consistently model the early-phase LCs.