Simulations of light curves and spectra for superluminous Type Ic supernovae powered by magnetars

TL;DR

This study uses detailed radiative transfer simulations to explore how magnetar-powered superluminous Type Ic supernovae evolve in light curves and spectra, highlighting the effects of magnetar properties and ejecta clumping.

Contribution

It provides the first comprehensive non-LTE time-dependent radiative transfer models of magnetar-powered SLSNe Ic, including effects of clumping and spectral diversity.

Findings

Magnetar power boosts ejecta temperature and ionization.

Clumping influences spectral line formation and diversity.

Most observed SLSNe Ic require ejecta clumping from early to nebular phases.

Abstract

Numerous superluminous supernovae (SLSNe) of Type Ic have been discovered and monitored in the last decade. The favored mechanism at their origin is a sustained power injection from a magnetar. This study presents non-local thermodynamic equilibrium time-dependent radiative transfer simulations of various single carbon-rich Wolf-Rayet star explosions influenced by magnetars of diverse properties and covering from a few days to one or two years after explosion. Nonthermal processes are treated; the magnetar-power deposition profile is prescribed; dynamical effects are ignored. In this context, the main influence of the magnetar power is to boost the internal energy of the ejecta on week-long time scales, enhancing the ejecta temperature and ionization, shifting the spectral energy distribution to the near-UV (even for the adopted solar metallicity), creating blue optical colors. Varying the ejecta and magnetar properties introduces various stretches and shifts to the light curve (rise time, peak or nebular luminosity, light curve width). At maximum, all models show the presence of OII and CII lines in the optical, and more rarely OIII and CIII lines. Non-thermal effects are found to be negligible during the high-brightness phase. After maximum, higher energy explosions are hotter and more ionized, and produce spectra that are optically bluer. Clumping is a source of spectral diversity after maximum. Clumping is essential to trigger ejecta recombination and yield the presence of OI, CaII, and FeII lines from a few weeks after maximum until nebular times. The UV and optical spectrum of Gaia16apd at maximum or the nebular spectrum of LSQ14an at +410d are compatible with some models that assume no clumping. However, most observed SLSNe Ic seem to require clumping from early post-maximum to nebular times (e.g., SN2007bi at +46 and +367d; Gaia16apd at +43d).

Tables

Table 1. Table 1: Summary of properties for the progenitor, the ejecta, the magnetar, as well as results at bolometric maximum obtained with the cmfgen  simulations. Numbers in parenthesis correspond to powers of ten. Each group of models (separated by a horizontal line) has a specific purpose. Models 5p11Bx2 and 5p11Bx2th address the influence of non-thermal processes. The reference magnetar model r0e2 is confronted to model r0e2n (no magnetar), to model r0e2s (narrow magnetar-power deposition profile), to models r0e1 and r0e4 (different explosion energies), to model r0e2ba (different magnetar field strength), and to models r0e2ea, r0e2eb, and r0e2ec (different magnetar spins). Also discussed is the influence of ejecta clumping for model r0e2 (Section  9 ).

Model	${}_{}{}^{56}N$i	$M_{\mathrm{ej}}$	$E_{\mathrm{kin}}$	$E_{\mathrm{pm}}$	$P_{\mathrm{pm}}$	$B_{\mathrm{pm}}$	$V_0$	$dV$	$t_{\mathrm{pm}}$	$t_{\mathrm{peak}}$	$L_{\mathrm{peak}}$
	[M⊙]	[M⊙]	[erg]	[erg]	[ms]	[G]	[km s-1]	[km s-1]	[d]	[d]	[erg s-1]
5p11Bx2th	0.19	3.63	2.49(51)	0.4(51)	7.0	3.5(14)	6000	3000	19.1	26.0	4.83(43)
5p11Bx2	0.19	3.63	2.49(51)	0.4(51)	7.0	3.5(14)	6000	3000	19.1	26.7	4.86(43)
r0e2	0.13	9.69	4.12(51)	0.4(51)	7.0	3.5(14)	5700	2850	19.1	31.7	3.75(43)
r0e2n	0.13	9.69	4.12(51)	…	…	…	…	…	…	44.4	2.25(42)
r0e2	0.13	9.69	4.12(51)	0.4(51)	7.0	3.5(14)	5700	2850	19.1	31.7	3.75(43)
r0e2s	0.13	9.69	4.12(51)	0.4(51)	7.0	3.5(14)	1425	713	19.1	51.7	2.81(43)
r0e1	0.09	9.52	1.14(51)	0.4(51)	7.0	3.5(14)	2700	1350	19.1	43.0	2.53(43)
r0e2	0.13	9.69	4.12(51)	0.4(51)	7.0	3.5(14)	5700	2850	19.1	31.7	3.75(43)
r0e4	0.17	9.86	1.23(52)	0.4(51)	7.0	3.5(14)	10000	5000	19.1	25.8	4.86(43)
r0e2	0.13	9.69	4.12(51)	0.4(51)	7.0	3.5(14)	5700	2850	19.1	31.7	3.75(43)
r0e2ba	0.13	9.69	4.12(51)	0.4(51)	7.0	1.0(14)	5700	2850	234.4	48.9	1.61(43)
r0e2	0.13	9.69	4.12(51)	0.4(51)	7.0	3.5(14)	5700	2850	19.1	31.7	3.75(43)
r0e2ea	0.13	9.69	4.12(51)	0.8(51)	5.0	3.5(14)	6800	3400	9.6	25.0	7.20(43)
r0e2eb	0.13	9.69	4.12(51)	1.2(51)	4.1	3.5(14)	8000	4000	6.4	21.6	1.03(44)
r0e2ec	0.13	9.69	4.12(51)	5.0(51)	2.0	2.0(14)	8000	4000	4.7	24.6	3.21(44)

Table 2. Table 2: Characteristics of the observed SLSNe Ic used in this paper, including the time of r 𝑟 r -band maximum, the redshift, the distance, the reddening, and the reference from where these quantities and observational data were taken.

	$t_{\max }\left(r\right)$	$z$	$d$	$E\left(B-V\right)$	Ref.
	[MJD]		[Mpc]	[mag]
SN 2007bi	54160.0	0.1279	591.6	0.03	a
SN 2010gx	55283.0	0.2297	1181.3	0.04	b
PTF12dam	56095.0	0.1073	462.8	0.0	c
PTF12gty	56143.4	0.176	750.0	0.058	d
LSQ14an	56595.0	0.1637	766.0	0.07	e,f
SN 2015bn	57108.0	0.1136	513.0	0.0	g,f
Gaia16apd	57554.0	0.1018	481.93	0.013	h

Equations

$${\dot{e}}_{\mathrm{pm}}=(E_{\mathrm{pm}}/t_{\mathrm{pm}})/{\left(1+t/t_{\mathrm{pm}}\right)}^2,t_{\mathrm{pm}}=\frac{6I_{\mathrm{pm}}{c}^3}{B_{\mathrm{pm}}^2{R}_{\mathrm{pm}}^6{\omega }_{\mathrm{pm}}^2},$$

Full text

1 Introduction

Over the last decade, numerous superluminous Type Ic supernovae (SNe) have been discovered and monitored extensively (see e.g., Quimby et al. 2011, Inserra et al. 2013, Nicholl et al. 2013, De Cia et al. 2018, Lunnan et al. 2018). These superluminous SNe (SLSNe) are generally found in low-mass low-metallicity host galaxies (Lunnan et al., 2014; Leloudas et al., 2015; Perley et al., 2016; Schulze et al., 2018). They reach peak luminosities of up to a few 1044 erg s*-1* and visual magnitudes of about $-21$ mag (Quimby et al., 2011; De Cia et al., 2018; Nicholl et al., 2017; Lunnan et al., 2018). Their rise times to bolometric maximum ranges from 20 to 125 d, while their post-maximum brightness decline rates cover from 0.01 to 0.06 mag d*-1*. At maximum, their optical spectra have a blue color. Apart from a few weak features associated with O ii (Pastorello et al., 2010; Quimby et al., 2011), these maximum-light optical spectra are nearly featureless. After maximum, the optical spectra show lines attributed to O i, Ca ii, and Fe ii (e.g., Gal-Yam et al. 2009; Nicholl et al. 2013), with similar properties to broad-line SNe Ic (Liu et al., 2017). At late times, their optical spectra exhibit forbidden lines of O i and Ca ii, as well as Fe ii lines, and are similar to the contemporaneous spectra of GRB SNe such as 1998bw (Nicholl et al., 2016a; Jerkstrand et al., 2017). However, SLSNe Ic also show O ii and O iii lines, which are strong in LSQ14an (Inserra et al., 2017) or PS1-14bj (Lunnan et al., 2016), and weak in SN 2007bi or SN 2015bn (Jerkstrand et al., 2017)). In rare instances, SLSNe Ic exhibit a strong and broad H$\alpha$ line at $100-200$ d after maximum (Yan et al., 2015, 2017a).

The current consensus is that SLSNe Ic owe their exceptional luminosities, color evolution, and spectral properties to the power contribution from a compact remnant (Bodenheimer & Ostriker, 1974; Kasen & Bildsten, 2010; Woosley, 2010). In standard core-collapse SNe, the explosion is thought to result from complex circumstances involving neutrino absorption and multi-dimensional fluid instabilities (for recent results from numerical simulations, see e.g., Lentz et al. 2015; Müller et al. 2017; Glas et al. 2018; O’Connor & Couch 2018; Vartanyan et al. 2019). The resulting SN luminosity then corresponds to the release of previously stored energy (left by shock passage) and the continuous release of energy by the decay of unstable isotopes, in particular ${}^{56}\rm{N}$i. Although expected to be rare, the iron core of a massive star may be fast spinning at the time of collapse (Hirschi et al., 2004; Yoon & Langer, 2005; Woosley & Heger, 2006; Georgy et al., 2009), establishing suitable conditions for the enhancement of the magnetic field during the collapse and the early post-bounce phase. Depending on the initial spin period and the magnetic field amplification, the proto-neutron star may lead to the formation of a proto-magnetar with a range of magnetic field strengths. For the strongest fields, the magneto-rotational effects may lead to a highly energetic explosion, perhaps associated with the formation of a baryon-free relativistic jet and a $\gamma$-ray burst (Usov, 1992; Wheeler et al., 2000; Akiyama et al., 2003; Thompson et al., 2004). For weaker fields, the magnetar spin-down timescale is longer, the dynamical impact of the magnetar is weaker, but its influence on the SN luminosity can then be large (Kasen & Bildsten, 2010; Woosley, 2010), and completely swamp the contribution from ${}^{56}\rm{N}$i. Figure 1 illustrates the predicted outcome for SN peak luminosity and rise time to peak as a function of magnetar spin period and field strength, and for an ejecta of 5 M*⊙* and a kinetic energy of 1051 erg (and no ${}^{56}\rm{N}$i). The maximum luminosity is reached for millisecond spin periods (i.e., large rotational energy budget) and moderate field strengths of 1014 G (corresponding to week-long spin-down timescales). In the context of the magnetar model, the current sample of SLSNe Ic yields inferred magnetar properties spanning the range $1.2-4$ ms initial spin periods, 0.2 to $1.8\times 10^{14}$ G field strengths, while the associated ejecta spans masses from 2.2 to 12.9 M*⊙, and kinetic energies from 1.9 to $9.8\times 10^{51}$ erg (Nicholl et al., 2017). Sophisticated numerical simulations of magnetar-powered SNe will probably not change these numbers much. Slower rotators and less magnetized magnetars can also appreciably alter SN radiation because the time-integrated luminosity of a SN is only 1049* erg, which is comparable to the rotational energy of a $\sim$ 40 ms period magnetar. Magnetars may also influence the properties of some He-rich explosions (e.g., SN 2005bf; Maeda et al. 2007) as well as some H-rich explosions (e.g., iPTF14hls; Dessart 2018). Future observations will probably extend the distribution of events lying between superluminous SNe and standard-luminous SNe (there may be a continuum of values for initial magnetar spins and magnetic field).

One-dimensional radiation-hydrodynamics modeling has been done to characterize the influence of the magnetar on the SN radiation and the ejecta structure (Kasen & Bildsten, 2010; Woosley, 2010; Bersten et al., 2016; Moriya et al., 2016; Tolstov et al., 2017). These simulations yield results in agreement with analytical predictions (Kasen & Bildsten, 2010). When extended to 2-D, hydrodynamical simulations exhibit strong differences with key features obtained in 1-D. In 2-D, the magnetar power leads to strong turbulence in the inner ejecta, which alters the shock properties and prevents the formation in the inner ejecta of a narrow dense and fast moving shell. The inner ejecta thus contains material down to low velocities rather than a low-density high temperature cocoon bounded by a fast-moving dense shell. In 3-D, the effect should persist but will probably be quantitatively different. The power released from the magnetar is likely aspherical and may not follow the smooth evolution of a dipole, which is generally assumed throughout the early magnetar life. Another concern is the efficiency with which the radiation from the magnetar is thermalized by the ejecta (see., e.g., Kasen et al. 2016). This power may eventually leak out from the ejecta, and indeed, leakage has been invoked to explain the late-time light curve (for the case of SN 2015bn, see Nicholl et al. 2018b).

Spectral modeling at the time of maximum provided strong evidence for energy injection from a compact remnant compared to a scenario in which the ejecta is influenced by an exceptionally large mass of ${}^{56}\rm{N}$i (Dessart et al., 2012b), later confirmed by similar modeling at the time of maximum (Howell et al., 2013; Mazzali et al., 2016) and at nebular times (Dessart et al., 2013; Jerkstrand et al., 2017).

This work is a follow-up to Dessart et al. (2012b), but this time using a time-dependent rather than a steady-state approach; the magnetar-powered SN can thus be followed from about a day until one or two years after explosion in a continuous fashion. This evolution is followed with cmfgen (Hillier & Dessart, 2012), using the approach of Dessart (2018) for the treatment of the magnetar power. The next section presents the carbon-rich Wolf-Rayet progenitors and Type Ic SNe models used, as well as the grid of magnetar properties considered. Section 3 discusses the results for the magnetar powered SN model r0e2. This model is used as a reference for additional comparisons and investigations. The subsequent sections investigate how the SN properties vary when various ejecta or magnetar characteristics are modified. Section 4 considers the influence of the adopted deposition profile for the magnetar power. Section 5 presents how a magnetar-powered ejecta is influenced by a change in ejecta kinetic energy. Sections 6 and 7 present the impact of varying the magnetar initial spin and field strength on the SN radiation and ejecta properties. Section 8 discusses the impact of the SN radiation when non-thermal processes are ignored (here, all simulations treat non-thermal processes unless otherwise stated). Section 9 quantifies the influence of clumping on the SN radiation, both at photospheric and nebular epochs. Section 10 presents a comparison between these models and a few observed SLSNe Ic. Section 11 presents the conclusions.

2 Numerical setup

2.1 Progenitor and explosion models

The simulations presented here are based on Type Ic SN ejecta that have been used in previous studies on SNe Ibc and GRB SNe (Dessart et al., 2015, 2016, 2017b). Two progenitor models and a variety of ejecta properties were used. One progenitor is model 5p11 from Yoon et al. (2010, model 31 in Table 1 of that paper) and the other progenitor is model r0 (Dessart et al., 2017b). These two models correspond to solar metallicity carbon-rich Wolf-Rayet stars that die with a final mass of 5.11 M*⊙* (rather than 4.95 M*⊙* quoted in Yoon et al. (2010); see Dessart et al. 2015 for discussion) and 11.4 M*⊙*. The influence of metallicity, and the use of a low metallicity progenitor more suitable for SLSNe Ic will be explored later (this is not essential given the large metal content of such Wolf-Rayet stars). These were exploded with v1d (Livne, 1993; Dessart et al., 2010b, a) by means of a piston to produce ejecta with different kinetic energies and ${}^{56}\rm{N}$i masses.

Model 5p11Bx2 (see Dessart et al. 2016 for details, including the mixing procedure, which is applied to all species) corresponds to an ejecta with a mass of 3.63 M*⊙, a ${}^{56}\rm{N}$i mass of 0.19 M⊙, and a kinetic energy of $2.49\times 10^{51}$ erg. The ejecta contains 0.34 M⊙* of helium, 0.9 M*⊙* of carbon, 1.38 M*⊙* of oxygen, 0.43 M*⊙* of neon, and 0.14 M*⊙* of magnesium.

Models r0e1/r0e2/r0e4 derive from progenitor model r0 and correspond to ejecta with a total mass $M_{\rm ej}$ of about 9.52/9.69/9.86 M*⊙, a ${}^{56}\rm{N}$i mass of 0.09/0.13/0.17 M⊙, and a kinetic energy $E_{\rm kin}$ of $1.14/4.12/12.3\times 10^{51}$ erg. Ejecta model r0 and its variants have a composition dominated by oxygen (a cumulative amount of about 5.4 M⊙), neon (about 1.9 M⊙), carbon (about 1.3 M⊙), and magnesium (about 0.35 M⊙). The helium content rises from 0.18 (r0e1) to 0.22 M⊙* (r0e4) because of the extra helium produced for higher energies in the inner envelope during explosive nucleosynthesis – most of the helium is however present in the outermost layers (see Fig. 2 for an illustration of the composition in the progenitor model r0 and of the explosion models r0e2 and r0e4; the left part of Table 1 summarizes the model properties). In these ejecta models based on the progenitor model r0, the mixing was applied only to ${}^{56}\rm{N}$i – the composition profile for other species was left as in the progenitor model (see Section 2 of Dessart et al. 2017b for discussion).

2.2 Magnetar power in cmfgen simulations

At 1 d after explosion, models 5p11Bx2, r0e1, r0e2, and r0e4 are remapped into cmfgen. These simulations use the same model atoms as in Dessart et al. (2015, 2016, 2017b). The grid typically consists of 100 points to cover evenly in optical-depth scale from 500 to $30000-50000$ km s*-1*. The simulations treat radioactive decay from the ${}^{56}\rm{N}$i chain. The energy deposition from radioactive decay is computed using a gray-absorption (opacity $\kappa_{\gamma}=0.06Y_{\rm e}$ cm2 g*-1*, where $Y_{\rm e}$ is the local electron fraction) solution to the radiative transfer of $\gamma$-rays (Swartz et al., 1995). Positrons are absorbed locally. For the magnetar power, the approach is similar to that presented in Dessart (2018) for the study of the Type II SN iPTF14hls.

Starting at day one (the magnetar power released prior to day one is ignored), a magnetar power $\dot{e}_{\rm pm}$ is introduced within the ejecta in homologous expansion, and is given by the expression

[TABLE]

where $E_{\rm pm}$, $B_{\rm pm}$, $R_{\rm pm}$, $I_{\rm pm}$ and $\omega_{\rm pm}$ are the initial rotational energy, magnetic field, radius, moment of inertia, and angular velocity of the magnetar; $c$ is the speed of light. In all cases, $I_{\rm pm}=10^{45}$ g cm2 and $R_{\rm pm}=10^{6}$ cm (see Kasen & Bildsten 2010 for details).

To determine how the magnetar power is deposited within the ejecta is difficult because the process is dynamical and multi-dimensional. In contrast, cmfgen does not treat dynamics (the ejecta structure is frozen) and is 1-D. In 1-D, a hydrodynamical simulation of a magnetar powered SN naturally leads to the formation of a dense shell if the energy is deposited in the innermost layers of the SN ejecta (Kasen & Bildsten, 2010; Woosley, 2010; Bersten et al., 2016; Dessart & Audit, 2018; Orellana et al., 2018). In cmfgen, adopting this deposition procedure would be unpractical since a huge temperature spike would form in the cmfgen simulation, producing large gradients that would be hard to resolve and would likely impede the convergence. It would also be physically inconsistent since the swept-up mass in the dense shell is accelerated as a result of magnetar power, while cmfgen cannot treat this dynamical effect (neither the snow-plow effect nor the acceleration).

In 2-D, hydrodynamical simulations indicate that this swept-up dense shell does not form (Chen et al., 2016; Suzuki & Maeda, 2017). Hence, the dense shell is an artifact of the 1-D treatment. Without an imposed spherical symmetry, the magnetar power triggers convection in the inner ejecta. This turbulent medium also causes advection of the magnetar energy faster than achieved by diffusion. An angle-dependent and radial-dependent (i.e., over a broad range of radii rather than in the innermost shells) energy deposition might also prevent the formation of a shock. The impact of the magnetar energy injection on the density structure is much weaker in 2-D (and probably also in 3-D) than in 1-D.

In the 1-D radiation hydrodynamical simulations with heracles, Dessart & Audit (2018) used this result to adopt a broad profile for the deposition of magnetar energy. A consequence is that then, even in 1-D, no dense shell forms. Furthermore, the broad deposition causes an earlier brightening of the SN luminosity, capturing in a simplistic way the energy advection seen in 2-D hydrodynamical simulations. This approach is therefore qualitatively sound but quantitatively uncertain. For example, in the process of advection of magnetar energy deposited at depth, the heated material would cool by stretching out in radius. By adopting a broad deposition profile, the cooling from expansion is not properly accounted for since the energy is deposited directly at large velocities.

The present simulations with cmfgen use a relatively extended energy deposition as parametrized in Dessart & Audit (2018) and Dessart (2018). Numerically, the deposition profile follows the density $\rho$ for $V<V_{0}$, and $\rho\exp\left(-[(V-V_{0})/dV]^{2}\right)$ for $V>V_{0}$. This profile is frozen in time. Figure. 3 shows two deposition profiles, one broad and one narrow. A normalization is applied so that the volume integral of this deposition profile yields the instantaneous magnetar power at that time.

The impact of the choice of $V_{0}$ and $dV$ is explored in one model (r0e2s; Section 4 and Table 1). In other simulations, $V_{0}$ and $dV$ scale with the magnetar energy or power since this seems plausible from the point of view of hydrodynamics. With this choice, the energy deposition profile influences the model luminosity mostly before maximum (Dessart & Audit, 2018). After maximum, the conditions are close to nebular so that it is the total magnetar power that primarily matters rather than its precise distribution within the ejecta.

While the $\gamma$-rays from radioactive decay can escape at late times, the entire magnetar power is designed to remain trapped at all times and is deposited according to the prescription given above. This approach is not optimal but a proper handing of the process would require transport modeling for the deposition of the high-energy radiation and particles released by the magnetar together with the multi-dimensional hydrodynamical modeling of its impact on the 3D structure and dynamics of the ejecta. At late times, one may adjust the trapping efficiency of the magnetar power so that the power absorbed equals the observed luminosity. This could be considered when comparing models to specific observations.

The approach in cmfgen currently requires that the composition is homogeneous at a given ejecta depth (or velocity). There is no chemical segregation. In the context of core-collapse SNe exploded through the neutrino-driven mechanism, this choice is not optimal (see e.g., Wongwathanarat et al. 2015). However, in the context of magnetar-powered SNe, the inner ejecta becomes very turbulent (Chen et al., 2016; Suzuki & Maeda, 2017) and chemical mixing may take place down to small scales. So, the assumption of complete microscopic mixing in cmfgen may in fact be more suitable than assuming no microscopic mixing at all. The issue is most relevant in Type II SNe because of the strong composition stratification of the progenitor massive stars. In Type Ic SN progenitor like model r0, the progenitor envelope is essentially a large ONeMg rich envelope on top of an iron core.

2.3 Grid of models

In this work, magnetars with a range of initial spin and magnetar field strengths are considered. The initial magnetar energy $E_{\rm pm}$ covers from 0.4 to $5.0\times 10^{51}$ erg. The magnetar field strength $B_{\rm pm}$ covers a small range from 1 to $3.5\times 10^{14}$ G. Larger values are a problem because the main effect from magnetar power would be dynamical, which is not treated in cmfgen. Lower values would produce brighter light curves at nebular times, but with fainter maxima so the SN would no longer be superluminous (this is still of interest but out of the scope of the present study). These choices of $E_{\rm pm}$ and $B_{\rm pm}$ correspond to magnetar spin-down timescales in the range $4.7-234.4$ d. The energy released during the first day, which is neglected, is $0.4-17.6$ % of the total magnetar energy, which is quite substantial for the highly magnetized fast-spinning models. Because the dynamical influence of magnetar power is ignored, all the energy liberated by the magnetar goes into ejecta internal energy (and eventually radiation). Throughout the time sequence, our magnetar-powered model will retain the same ejecta structure (mass, density, velocity, and radius) and thus the same kinetic energy.

In practice, the code treats the magnetar power the same way as radioactive decay. Energy is injected as 1 keV electrons for which the degradation spectrum is computed. The contribution to the heating of thermal electrons, and to non-thermal excitation and ionization is computed explicitly in all but one model – one simulation is run with non-thermal processes turned off (i.e. all the energy is in this case entirely channeled into heat).

The model nomenclature is the following. The same magnetar properties ($B_{\rm pm}=3.5\times 10^{14}$ G and $E_{\rm pm}=4\times 10^{50}$ erg) are used in the “reference” models 5p11Bx2, r0e1, r0e2, and r0e4. There is no specific suffix to characterize the magnetar in these 4 models. The model r0e2 is extensively discussed and used for various comparisons. The properties of its magnetar are the same as those used in the study of Dessart & Audit (2018) — this particular choice was in part arbitrary but also driven by the need to use a magnetar with a modest dynamical influence (i.e., with $E_{\rm pm}$ significantly smaller than $E_{\rm kin}$, which is $\sim$ 10% in the case of model r0e2). Variants of these models are then explored. Model 5p11Bx2th is the same model as 5p11Bx2 except that non-thermal processes are switched off (all magnetar and decay power goes into heat). Model r0e2n is equivalent to model r0e2 except that the magnetar power is turned off. Model r0e2s is equivalent to model r0e2 except that the magnetar energy deposition profile is narrower. Model r0e2ba is equivalent to model r0e2 except that the magnetar field strength is lower. Models r0e2ea, r0e2eb, and r0e2ec are equivalent to model r0e2 except that the initial magnetar rotational energy is enhanced (model r0e2ec also uses a weaker field). Model r0e2cl is equivalent to model r0e2 except that the material is clumped. This simulation is started around the time of maximum. In addition, in models r0e1, r0e2, and r0e4, various levels of clumping are enforced to explore the effects at one nebular epoch (details and results are given in Section 9).

This set of models includes only one magnetar with a large rotational energy because in that regime, the neglect of dynamical effects in cmfgen  becomes a severe limitation. Consequently, most of the present simulations are limited in maximum luminosity and optical brightness. And the quantitative results are less robust than the trends.

Figure 4 shows the bolometric light curve for the full set of time-dependent cmfgen  simulations. The peak luminosity covers from about $10^{43}$ up to few times $10^{44}$ erg s*-1*, with rise times from about 25 up to about 50 d after explosion. The slight offset at $100-200$ d between the model curves and the magnetar power is related to the small contribution from radioactive decay and residual optical-depth effects (i.e., there is still stored energy released in addition to the instantaneous power absorbed by the ejecta). Eventually, the bolometric luminosity equals the instantaneous magnetar power. The contribution from decay power depends on the ${}^{56}\rm{N}$i mass of each model (see Table 1), as well as the amount of $\gamma$-ray escape (function of $E_{\rm kin}/M_{\rm ej}$).

3 Detailed discussion for magnetar-powered SN model r0e2 – the reference case

Figure 5 illustrates the photometric properties of the model r0e2 with a magnetar ($B_{\rm pm}=3.5\times 10^{14}$ G and $E_{\rm pm}=4\times 10^{50}$ erg) relative to its counterpart r0e2n without a magnetar. Model r0e2 reaches a peak luminosity of $3.75\times 10^{43}$ erg s*-1* at 31.7 d after explosion while model r0e2n peaks later at 44.4 d at a smaller peak luminosity of $2.26\times 10^{42}$ erg s*-1*. The $V$-band light curve is similar to the bolometric light curve (middle panel of Figure 5). The $V$-band maximum occurs in both at 43.3 d but it is $-19.4$ mag in model r0e2 and $-17.5$ mag in model r0e2n. The width of the light curve (bolometric or visual) is much broader in the magnetar-powered model, because the magnetar power in model r0e2 (black thin line Fig. 5) overwhelms the decay power from 0.13 M*⊙* of ${}^{56}\rm{N}$i (dashed thin line in Figure 5; $\gamma$-ray escape is allowed for) at all times. The boost in luminosity or visual brightness in model r0e2 is accompanied by a shift to blue optical colors, which is maintained throughout the photospheric phase (bottom panel of Fig. 5).

Figure 6 shows the evolution of the photospheric properties (the electron scattering opacity is used for the computation of the optical depth and the location of the photosphere). Although models r0e2 and r0e2n have the same ejecta mass, kinetic energy, density structure, the magnetar-powered model r0e2 has a photosphere at larger velocities $V_{\rm ph}$, (which extends to larger radii, with a maximum of about $4\times 10^{15}$ cm at 120 d) and larger temperatures $T_{\rm ph}$. Around the time of maximum, $T_{\rm ph}$ is $\sim$ 12,000 K in model r0e2, compared to $\sim$ 6000 K in model r0e2n. The higher temperature in model r0e2 causes a greater ionization of the photospheric layers, so that oxygen is once ionized from 4 d onwards in model r0e2, while it is neutral up to 20 d and only half ionized afterwards in model r0e2n (bottom panel of Fig. 6).

Since the spectrum formation region tracks approximately the photospheric layers111This is true up to bolometric maximum. After that, the spectrum forms over an increasingly extended region in the present simulations, in part because the ionization remains large above the photosphere so that no recombination front forms. as long as the ejecta is optically thick (hence up until $\sim$ 200 d after explosion in model r0e2, and $\sim$ 140 d in model r0e2n), the contrast in photospheric properties shown in Fig. 6 should lead to a modest change in line profile width and strength but a much bluer optical spectrum with lines from once ionized elements (if one considers that O is a representative element). Figure 7 confirms this. The magnetar powered model shows a rapid shift from a red to a blue optical spectrum after a few days past explosion, with a spectral energy distribution (SED) that peaks around 2500 Å (Fig. 7 shows the quantity $\lambda F_{\lambda}$, which gives more importance to the long-wavelength flux) from $\sim$ 10 to beyond 100 d. The onset of that shift to the blue (and to greater luminosity and visual brightness is function, for example, of the magnetar power and deposition profiles — see next sections for various dependencies). The contrast with the model without magnetar is striking (right panel of Fig. 7). In model r0e2n, the SED stays red at all times, with strong signs of blanketing in the UV and in the optical (model 5p11Bx2 and variants have been extensively discussed in Dessart et al. 2015, 2016). Spectral differences between models r0e2 and r0e2n remain throughout the photospheric and nebular phases.

The left panel of Figure 7 also shows the predicted continuum flux (thin line). At the earliest epoch (here, at 5 d), the continuum flux follows the optical flux closely (outside strong-line regions), but it is strongly blanketed in the UV. As time proceeds, the continuum peaks more and more in the blue and the continuum flux weakens relative to the total flux. By 70 d, the continuum flux is essentially zero. The difference between the total flux and the continuum flux arises from the multitude of broad and weak overlapping lines, which cannot be clearly identified in the full spectrum. These lines contribute both in emission and absorption, although as time passes, the emission part dominates. A consequence is that the color temperature obtained by fitting the full spectrum with a blackbody overestimates the temperature of the thermalization layer. A lower gas temperature can accommodate the observed flux because of the line flux contribution (which a blackbody ignores). In these conditions, blackbody fits, which are often used, overestimate the temperature. The weakness of the continuum flux arises from the metal-dominated composition and the weakness of continuum processes relative to line processes. A similar effect holds in SNe Ia (see e.g., Hillier et al. 2013).

Figure 8 shows the evolution of the spectral region centered on the C  ii 6579.7 Å (doublet) transition. This line is strong (to assess the strength of the line, see Figs. 7 and 9) and appears clearly as a P-Cygni profile as soon as the photospheric conditions are hot and ionized (this line is weak or absent in the model r0e2n at the corresponding epochs). The Doppler velocity at the maximum absorption in this feature tracks closely the velocity of the electron-scattering photosphere for up to 70 d (for this illustration, this C  ii doublet is a better line to use than those in the $4000-5000$ Å region, where O  ii lines are numerous and overlap). The location of maximum absorption stays the same from about 60 d until 100 d because the photosphere then recedes to inner ejecta layers where C is under-abundant: the C-rich shell is located beyond 6000 km s*-1* in model r0e2 (analogous to the effect on the H$\alpha$ line profile caused by the H stratification in SNe IIb; see Fig. 2). After $\sim$ 100 d, the C  ii lines disappear. In our WR star models, a C-rich shell is always present in the outermost layers of the ejecta and is the cause of the strong C  ii lines in our models.

Figure 9 uses two different ways to illustrate the main line contributions in the spectra of models r0e2 at 32.0 and 240.0 d after explosion. The top panel shows the contribution of individual ions, obtained by comparing the full spectrum ($F_{\lambda,{\rm All}}$) with the spectrum computed by omitting the bound-bound transitions of a given ion ($F_{\lambda,{\rm less}}$). The bottom panel marks the strongest transitions with a Sobolev equivalent width (EW) greater than 7000 Å (at nebular times, the very weak continuum flux and the presence of strong forbidden emission lines explain these large EWs). The top panel of Fig. 9 shows the contributions from a large number of lines associated with once and twice ionized states of C, O, Mg, and Fe. At 32.0 d, the spectrum forms around 10000 km s*-1*, in the C-rich part of the progenitor CO core (see middle panel of Fig. 2 and Fig. 6), hence the spectrum naturally shows signatures of these intermediate mass elements, although the number of line transitions is large and line overlap is strong, often preventing an association between a feature and a given transition. The photospheric temperature and ionization favor the presence of these elements in their once or twice ionized state. The region $4000-5000$ Å is primarily shaped by O ii  which has about one hundred lines with an EW greater than 100 Å (3p–3s, 3d–3p, and 4f–3d transitions). At longer wavelength, the strongest features are due to C ii, mostly doublets or triplets at 5889.7, 6579.7, 7234.8 Å. He i 5875.6 Å contributes weakly (about 1/5 of the total Sobolev EW) to the feature at 5900 Å. O i 7773.4 Å is the main contributor of the feature at 7800 Å, with a weak contribution at 7896.4 Å by Mg ii – the main Mg ii contribution is from multiple lines around 2800 Å, which also overlap with C ii and C iii transitions. There are weaker Mg ii transitions around 2930 and 4480 Å. In the UV range, most lines are caused by similar ions, together with a few regions of strong blanketing by iron-group elements (Ti iii around 1500 Å; Fe iii in the range $1700-2200$ Å). The ups and downs in the UV range are not actual lines, not even blends of lines. These peaks and dips are regions of reduced and enhanced opacity, and the entire UV flux is subject to blanketing.

The bottom panel of Fig. 9 shows the spectrum of model r0e2 at 240 d. The total ejecta optical depth to electron scattering at that time is 0.48, and the conditions for the spectrum formation are hybrid, with the presence of permitted transitions (similar to those seen at 24.5 d, as well as a forest of Fe ii lines, including strong lines at 6149.3, 6247.6, and 6456.4 Å — there are numerous weaker components overlapping with these transitions) and forbidden lines, most notably of O  ii 7324.3 Å (multiplet) and O iii at 5006.8 Å. At later times, the permitted transitions weaken and the spectrum becomes more and more dominated by forbidden lines (which ions contribute depend on how the ionization state evolves). Oxygen absorbs most of the magnetar power (the power absorbed scales roughly linearly with the mass of each element), although oxygen lines are not always the strongest lines (this depends on the composition mixture and on atomic properties of the ions and atoms present).

4 Influence of the magnetar energy deposition profile

This section discusses the influence of the adopted deposition profile. Model r0e2s uses $V_{0}=1425$ km s*-1* and $dV=713$ km s*-1*, compared to $V_{0}=5700$ km s*-1* and $dV=2850$ km s*-1* in model r0e2 (see Fig. 3 and Section 2.2). The resulting SN ejecta and radiation properties computed with cmfgen are qualitatively similar but there are quantitative differences (Fig. 10). At early times, the magnetar influence is delayed in model r0e2s, so that the increase in luminosity and visual brightness and the shift of the optical color to the blue are all delayed. The sudden rise in luminosity occurs at $\sim$ 3 d in model r0e2 but is delayed until nearly 20 d in model r0e2s. From model r0e2 to r0e2s, the rise time increases from 31.7 to 51.7 d, the bolometric maximum drops by 25%, the $V$-band magnitude rises by 0.2 mag. During the rise to maximum, model r0e2s stays red longer but it has roughly the same $U-V$ and $V-I$ colors at maximum as model r0e2. In model r0e2, more energy is radiated (i.e., the time-integrated luminosity is greater) because the magnetar power is released over a broader range of velocities. By depositing more power at smaller optical depth, radiation is allowed to escape earlier with a reduced degradation by expansion. As discussed in Section 2.2, this approach is not fully consistent since magnetar energy advected out would be subject to expansion cooling. So, this neglect tends to overestimate the influence of magnetar power. However, our neglect of dynamics means that the inner ejecta layers are not accelerated, which would contribute to bringing them faster to a lower optical depth where they could radiate their stored energy (mass conservation also implies that this swept-up shell becomes very dense; a radial compression, however, does not change the radial column density).

Figure 11 illustrates the impact of the deposition profile on the photospheric properties. With a narrower profile, the photosphere is located deeper in (smaller radii and velocities), is cooler and less ionized on the rise to maximum. Around the time of maximum (which differs by 20 d between the two models), the difference in photospheric properties is small, mostly limited to the temperature ($\sim$ 11000 K in r0e2 compared to $\sim$ 9000 K in model r0e2s). This offset likely results from the longer rise time in model r0e2s, which implies that the volume occupied by the ejecta in model r0e2s is greater at maximum (smaller equilibrium radiation temperature in the optically-thick volume) and the magnetar power at that later time is smaller. After maximum, the deposition profile has less impact on the SN properties because the ejecta is not very optically thick, the spectrum forms over a more extended volume, and the outer ejecta layers contribute less.

The difference in spectral evolution reflects the delay in the magnetar influence, especially the shift to a hotter and more ionized photosphere radiating a bluer optical color (top panel of Fig. 12). Model r0e2s exhibits a similar spectrum to model r0e2n (without magnetar) for about 30 d. By 40 d, the spectra of models r0e2 and r0e2s show similar lines in the optical, while model r0e2 has a greater flux in the UV (its photosphere is hotter). Line widths are comparable around maximum (the photospheric velocity is similar in both models from about 30 to 100 d; Fig. 11). At nebular times, model r0e2s shows narrower spectral lines, which therefore form deeper in the ejecta where the deposited power is relatively greater. A marked difference that persists at all times is the excess UV flux in model r0e2 compared to r0e2s. This arises from the excess power deposited in the outer lower-density ejecta layers, which tends to produce a greater temperature (and ionization).

The bottom panel of Fig. 12 illustrates the spectral features in model r0e2s at 427 d after explosion. The spectrum is very ionized relative to a core collapse SN. Here, the main lines are due to O  ii 7324.3 Å, O iii at 4363.2, 4958.9, 5006.8 Å, S iii at 6312.0, 9068.6, and 9530.6 Å, Ar iii 7135.8 Å, Ca ii H & K, Fe ii  at 4923.9, 5018.4, 5169.0 Å, and Fe iii at 4658.0, 4701.5, and 5011.3 Å. All these lines are also present in model r0e2 (with different relative strength), but they are more easily identified here with the reduced line overlap.

The magnetar-power deposition profile therefore impacts the SN radiation in many ways. It controls the timing for the influence felt at the photosphere as well as the different reaction of the gas for a given power (i.e., whether the power is injected at higher or lower density). These changes cause various quantitative shifts but do not change the qualitative picture. Overall, the primary impact is on the total heat content, which is controlled by the volume-integrated magnetar-power. The exact distribution is a secondary effect.

5 Influence of ejecta kinetic energy

This section presents the results for counterparts of model r0e2 for which the ejecta kinetic energy is lower (model r0e1) and higher (model r0e4; see Section 2). By adopting the same magnetar properties in all three models, one can test the influence of the ejecta kinetic energy on the SLSN Ic properties. For the magnetar power deposition profile, model r0e1 uses $V_{0}=2700$ km s*-1* and $dV=1350$ km s*-1*, model r0e2 uses $V_{0}=5700$ km s*-1* and $dV=2850$ km s*-1*, and model r0e4 uses $V_{0}=10,000$ km s*-1* and $dV=5000$ km s*-1* (see Section 2 and Fig. 3).

Figure 13 shows the bolometric light curve, the $V$-band light curve, and the $U-V$ and $V-I$ color curves for models r0e1, r0e2, and r0e4. The effect seen here is similar to employing a broader deposition profile (see previous section). With higher ejecta kinetic energy, the luminosity (or optical brightness) increases earlier (because of the lower ejecta optical depth and the broader deposition profile), peaks to a greater maximum (because it peaks earlier, when the magnetar power is greater), and is bluer early on. With higher ejecta kinetic energy, the photosphere is located further out (larger radii and velocity), although the ejecta becomes optically thin earlier (Fig. 14). At the time of bolometric maximum, model r0e4 ($E_{\rm kin}=1.23\times 10^{52}$ erg) has a photospheric velocity of $\sim$ 14000 km s*-1*, while model r0e1 ($E_{\rm kin}=1.14\times 10^{51}$ erg) has a photospheric velocity of $\sim$ 5800 km s*-1*. This reflects closely the difference in $\sqrt{2E_{\rm kin}/M_{\rm ej}}$. The photospheric temperature is greater because of the earlier influence of the magnetar, but all models have a similar $T_{\rm ph}$ around the phase of maximum.

Figure 15 shows the spectral evolution for models r0e1, r0e2, and r0e4. With higher ejecta kinetic energy, the SED becomes bluer earlier, the lines are broader and weaker, and the appearance of forbidden lines occurs earlier (e.g., O  ii 7324.3 Å and O iii at 4363.2, 4958.9, 5006.8 Å). While the ionization conditions are similar around maximum (dominance of permitted transitions from C ii  and O ii), the conditions at late times are different. Model r0e4 shows nearly exclusively lines of O ii and O iii while model r0e1 shows lines primarily from neutral and once-ionized elements (e.g., O i, Si ii, and a forest of Fe ii lines). Figure 16 illustrates this contrast for models r0e1 and r0e4 at about 240 d after explosion.

Hence, for the same magnetar properties, the nebular phase spectra can drastically differ. The origin of this feature is that in a higher energy explosion, the inner ejecta is less dense (the density profile is flatter, with more mass at large velocity) so for a given power one obtains a greater temperature and ionization at nebular times. This is exacerbated by the fact that the deposition profile is broader in model r0e4, leading to very broad O ii and O iii lines. In model r0e1, the lines are from neutral or once-ionized species like O i or Fe ii and they are narrower.

6 Influence of magnetar rotation energy

This section explores the impact of varying the magnetar initial rotational energy (or spin) for the same ejecta model. Using model r0e2 as a reference, the value $E_{\rm pm}$ is increased from 0.4 (r0e2) to 0.8 (r0e2ea), 1.2 (r0e2eb), and $5.0\times 10^{51}$ erg (r0e2ec; in this last model, $B_{\rm pm}$ is decreased from 3.5 to $2.0\times 10^{14}$ G). The change in initial spin also affects the spin-down timescale since it scales as $1/\omega_{\rm pm}^{2}$ (or $1/E_{\rm pm}$). Hence, $t_{\rm pm}$ drops from 19.1 d in model r0e2, to 9.6, 6.4, and 4.7 d (model r0e2ec, with the fastest initial spin). This causes a boost of about a factor of ten in peak bolometric luminosity (top panel of Fig. 17), while the $V$-band maximum increases by less than 2 magnitudes. This results from the shift to bluer colors in the highest energy model (bottom panel of Fig. 17). The rise time to maximum is essentially the same because all models have the same expansion time scales and diffusion times scales (here the magnetar does not affect the dynamics) — this is modulated in part by the fact that a broader deposition profile is used in higher energy models).

As shown in Fig. 18, the rise time to maximum and the luminosity at maximum obtained here are in rough in agreement with the analytic predictions of Kasen & Bildsten (2010), who use a one-zone model (compared to an extended deposition profile here) and neglect dynamics (although they account for the magnetar energy in estimating the mean ejecta velocity).

The boost to the luminosity and brightness obtained with more energetic magnetars is associated here with a boost in the photospheric temperature (Fig. 19). Indeed, the photosphere is located in the same mass shells (i.e., same radius or velocity), but the temperature rises to a maximum of $\sim$ 12000 K (at bolometric maximum) in model r0e2 and to $\sim$ 27000 K (about 10 d before bolometric maximum) in model r0e2ec. In model r0e2ec, the photospheric temperature is always greater than in model r0e2, by 2 to 10 kK. This temperature contrast is associated with a change in ionization. While O is once ionized in model r0e2 throughout the simulation, it evolves from neutral to 3-times ionized on the way to maximum in model r0e2ec, before progressively dropping to being once-ionized at late times.

Models with different magnetar rotational energy show very different spectral properties, reflecting in part the contrast in temperature and ionization discussed above. The top panel of Fig. 20 shows multi-epoch spectra (the quantity displayed is $\lambda F_{\lambda}$) for models r0e2 and r0e2ec (which differ in $E_{\rm pm}$ by a factor of 50). The higher energy model has an SED shifted to the blue, peaking in the UV at all times (the flux bias toward the UV region would look more extreme if showing $F_{\lambda}$). Because of this higher ionization, model r0e2ec shows lines of C iii (e.g., 1175.6, 1247.4, 1894.3, 2296.9, 2844.1, 4650.2, 5695.9, and 8500.3 Å) and O iii (e.g., at 4363.2 and 5006.8Å), in addition to the C ii and O ii lines that dominate the optical range in model r0e2 (see Section 3). In the UV, apart from the strong C iii lines, Fe iii causes blanketing in the range $1700-2100$ Å. The bottom panel of Fig. 20 illustrates the line contributions in model r0e2ec around the time of maximum.

7 Influence of magnetar field strength

Model r0e2ba explores the effect of a weaker magnetic field relative to model r0e2 ($10^{14}$ G compared to $3.5\times 10^{14}$ G). This yields a greater magnetar spin-down timescale (which scales as $1/B_{\rm pm}^{2}$), increasing from 19.1 in model r0e2 to 234.4 d in model r0e2ba. The magnetar power is thus smaller early on but greater after about 100 d for the weaker field case. This causes the bolometric and visual (not shown) light curves to peak later and at a fainter maximum in model r0e2ba (Fig. 21), and the color shift to the blue is also delayed (not shown). Compared to model r0e2, the photosphere in model r0e2ba heats up later and becomes more ionized later, but then the greater power at late times maintains a higher temperature and a higher ionization for longer (not shown).

Spectroscopically (no figure provided), model r0e2ba shows a similar evolution to model r0e2s early on (because the influence of the magnetar is delayed due to the smaller power initially). At late times, the enhanced power in model r0e2ba yields a highly ionized spectrum, similar to what is found in model r0e4 but with narrower lines. In model r0e4, the high ionization comes from the deposition of power in a lower density ejecta (because of the large explosion energy relative to model r0e2) while in model r0e2ba, the density is higher (the explosion energy is smaller, equal to that in model r0e2) but the magnetar power is greater.

Dessart (2018) used a similar approach to that presented here, but for a blue-supergiant star explosion model, and showed results for a SN ejecta influenced by an even weaker magnetar field ($7.0\times 10^{13}$ G). That model predicts a sustained luminosity, optical brightness, and blue color for 2 years, as observed in iPTF14hls (Arcavi et al., 2017). For these weaker magnetic fields, the longer spin down timescale implies that the dynamical effects are weaker and a smaller fraction of the magnetar energy is released at early times, making the present assumptions in cmfgen more suitable.

8 Influence of non-thermal processes

Model 5p11Bx2 is used to explore the influence of non-thermal effects on the SN ejecta and radiation. So, a counterpart is run (model 5p11Bx2th; the sequence covers from 3 to 100 d only) in which all the magnetar and decay energy is treated as heat (i.e. non-thermal effects are ignored).

Figure 22 shows multi-epoch spectra for the two model sequences at 5, 20, 40, 70, and 100 d. As can be seen, the impact of non-thermal processes throughout the photospheric phase (i.e., the high brightness phase) is negligible. A slight difference is visible at 100 d in the O ii feature at about 7324.3 Å (composed of several lines), and arises from the higher ionization in the non-thermal model (a mix of O+ and O2+ is present at the photosphere), while the thermal model is less ionized (O+ dominates at the photosphere). The photospheric properties (not shown) are essentially identical in both models, with only a slight over-ionization in the non-thermal model after bolometric maximum. At later times, the offset between the two models grows.

In model 5p11Bx2, all elements are at least once ionized throughout the ejecta at all times so the electron density is very large. This is known to quench non-thermal effects. Indeed, Kozma & Fransson (1992) find that as the ionization level increases, a growing fraction of the decay energy is channeled into heat (i.e., heating of thermal electrons), largely irrespective of the composition. In cmfgen simulations of Type II SN 1987A, Li et al. (2012) find that in the recombined layers of the H-rich ejecta ($x_{\rm e}\sim$ 0.01, where $x_{\rm e}$ is the ratio of the total ion population and the total element [atom and ions] population), 20% of the decay energy goes into heat, 40% goes into excitation, and 40% into ionization of the elements. In the partial ionization conditions of Type Ibc photosphere around bolometric maximum ($x_{\rm e}\sim$ 0.3), Dessart et al. (2012a) find that 80% of the decay energy goes into heat, and 10% goes into excitation and 10% into ionization of the elements. In the magnetar-powered model 5p11Bx2 at 30 d (close to bolometric maximum), $x_{\rm e}$ is 0.95 at the photosphere and drops to about 0.6 within a few 1000 km s*-1* above it. In this spectrum formation region, 85 to 91% of the magnetar power is channeled into heat, with the remaining $\sim$ 10% are channeled equally between non-thermal excitation and ionization. With the increase in ionization and temperature, the radiation field is also much stronger in the UV (the SED peaks in the UV in SLSNe Ic), allowing stronger photo-ionization than typically possible in cool plasmas like standard SNe Ibc photospheres. In early-time Type II (or IIb) SNe, emission from lines of elements in a relatively high ionization state have been observed [Ov and Ovi in SN 2013fs (Yaron et al., 2017); N iii and N iv in SN 2013cu (Gal-Yam et al., 2014); He ii in SN 2006bp (Quimby et al., 2007)]. Rather than non-thermal processes, these emission lines are caused by the high temperature in the spectrum formation region (Dessart et al., 2008; Gräfener & Vink, 2016; Groh, 2014; Dessart et al., 2017a). The combination of a high ionization (which favors heating of thermal electrons) and the large gas temperature (which boosts the UV flux) explain why non-thermal processes have no visible impact on the spectra. This also confirms the earlier spectral simulations of SLSNe Ic by Dessart et al. (2012b), which ignored non-thermal processes, but nonetheless reproduced the salient features of SLSNe Ic like PTF09atu (Quimby et al., 2011). Although this depends on the exact level of ionization of the gas, this seems in contradiction with Mazzali et al. (2016), who argue that non-thermal effects are important for the spectrum formation of superluminous SNe Ic, including for the formation of He i lines. The present simulations, which solve explicitly for non-thermal processes, predict that “thermal” processes dominate (this does not mean that the conditions are in LTE).

He i lines are present in the synthetic spectra computed here but generally limited to the transitions at 5875.7 and 10830.2 Å. In models 5p11Bx2 and r0e2, He i 5875.6 Å represents about 1/5 of the 5890 Å feature, whose strength is dominated instead by C ii. This (and other) C ii features are present for as long as He i 5875.6 Å, and the carbon contribution always dominates in the present models. In the near-IR, He i 10830.2 Å overlaps with O i, O ii, and Mg ii lines. The helium contribution is weak early on, but it strengthens after maximum, and at late times non-thermal effects influence that line. The photospheric temperature of about 10000 K during the high brightness phase is, however, the main reason for the presence of these lines at that time. The weakness of He i lines is probably due to the low helium abundance in the outer ejecta ($5-10$ % in model r0e2, and 30% in model 5p11Bx2, the rest being primarily carbon and oxygen). Given the large CO core in model r0e2, the helium is located too far out in the ejecta, yielding low line optical depths and consequently weak features. The progenitor CO core mass is a fundamental quantity controlling the production of a Type Ib or a Type Ic SN (Dessart et al., 2012a).

9 Influence of clumping

This section presents an exploration of the effect of ejecta clumping. The treatment of clumping is identical to that in Dessart et al. (2018). Clumping is treated through a volume-filling factor approach which does not affect the radial column density, hence it does not facilitate the escape of photons, whether of low energy (optical photons) or high energy ($\gamma$-ray photons). A value of $f_{\rm vol}$ equal to 1 corresponds to a smooth ejecta, and a value of 0.1 corresponds to 10% of the volume being filled with material (the rest being vacuum). Whatever the level of clumping, the composition is considered homogeneous at a given velocity (i.e., chemical segregation allowing for clumps of distinct composition is not considered). In Dessart et al. (2018), the main effect of clumping was to enhance the recombination efficiency of the gas at and above the photosphere, leading to a faster recession of the photosphere and a boost to the luminosity (in the context of a Type II-pec or a Type II-P SN). Here, clumping is studied for its influence on the ionization state of the gas. Indeed, Jerkstrand et al. (2017) found that a strong clumping was required to explain the nebular-phase spectra of superluminous SNe. Specifically, a volume-filling factor of 0.1% for the ONeMg rich material was required to reproduce, for example, the O i line doublet at 6316.0 Å and quench the emission from O ii and O iii lines (i.e., shift the ionization from ionized oxygen to neutral oxygen).

The next section explores the impact of various levels of clumping in model r0e1, r0e2, and r0e4 at a single nebular-phase epoch (i.e., 240 d after explosion). The subsequent section presents a time-dependent simulation for model r0e2cl, which is identical to model r0e2 but differs in that the ejecta is clumped uniformly with a 10% volume-filling factor.

9.1 Comparison of models at nebular times for different clumping factors

Figures 23 and 24 show the impact of various levels of clumping ($f_{\rm vol}$ of 1, 0.1, 0.03, and 0.01) on the ejecta temperature and oxygen ionization, as well as on the spectrum, for model r0e2. For a greater clumping (smaller value of $f_{\rm vol}$), the gas is cooler and more recombined (oxygen is a good represent of the ionization state of the ejecta material since it dominates the composition). The smooth model is more ionized and cools through emission in O ii or Fe ii lines. The clumped models are more recombined (already for $f_{\rm vol}=0.1$), with a dominance of neutral oxygen over ionized oxygen. For example, the doublet line of O i at 6316.0 Å is strong in the three clumped models, but absent in the smooth model, where the 6300 Å feature is due instead to Fe ii lines. A similar shift is seen for the near-IR Ca ii triplet, which is stronger with increasing clumping as the Ca ionization shifts to being primarily Ca+. The SED is markedly altered in clumped models, which appear redder.

In model r0e2, the effect of clumping is strong for $f_{\rm vol}=0.1$, and not much enhanced as the clumping factor if further enhanced. The value of $f_{\rm vol}$ needed to cause the ionization shift depends on the SN age, the ejecta properties (primarily density, which depends on the progenitor structure, the explosion energy), and the magnetar properties (power, deposition profile). Jerkstrand et al. (2017) performed one-zone models for various compositions and found that $f_{\rm vol}$ needed to be as low as 0.001 to yield a satisfactory correspondence to observed spectra of SLSNe Ic. This seems extreme since the O-rich material represents most of the ejecta mass (which would imply that most of the ejecta volume is essentially vacuum), but it is clear that clumping can have a profound impact on the resulting SN spectrum.

The impact of clumping depends on the ejecta properties. Indeed, in general, the ionization state of the gas for a given power depends on the density so one expects a different result for clumping in a magnetar power acting on a low or a high energy explosion. Figure 25 shows the impact of clumping in models r0e1 and r0e4. Evidently, the impact is rather modest in model r0e1, because this lower energy model is characterized by a dense inner ejecta. The ejecta ionization is rather low even in the smooth case. Introducing clumping does not alter this much, although one sees that the O i 6316.0 Å line is stronger in the clumped model. In the higher energy model r0e4, the situation is completely different. In this case, the ionization of the smooth ejecta is very high at 240 d (see also Section 5). In this case, clumping leads to a drastic change in ejecta ionization and spectral morphology. The strong and broad O ii and O iii lines present for $f_{\rm vol}=1$ (top right panel of Fig. 25) are absent in the clumped model (bottom right panel of Fig. 25) with $f_{\rm vol}=0.01$, which shows lines of O i (as well as Fe ii and Ca ii).

In model r0e4, if the ejecta is smooth, a large fraction of the flux is emitted in O ii  7323.4 Å, O iii 4363.2, and O iii 5006.8 Å. In the clumped ejecta model with $f_{\rm vol}=$ 0.01, the same amount of magnetar power is absorbed by the oxygen rich material (because the deposition profile is independent of clumping), but the oxygen line flux is much smaller. Instead, the bulk of the flux comes out in Fe ii lines. The reason for this feature is probably that the main O i forbidden line at 6316.0 Å is a poor coolant compared to O ii  7323.4 Å or O iii 5006.8 Å. The oscillator strengths are 10 to 100 times greater for these last two lines. So, if O+ and O2+ are present, the oxygen-rich material cools primarily through O ii and O iii lines, but if neutral O dominates, the cooling here is then performed in a large part by iron lines. This vividly illustrates that for the same oxygen mass and the same power absorbed by oxygen rich material, very different line strengths of oxygen lines can result, strong if O ii and O iii lines are present, weak if O i lines are present. In other words, the strength of oxygen lines is not a robust and direct measure of oxygen abundance. Allowing for chemical segregation might alter this behavior.

As discussed in Section 4, the nebular-phase properties also depend on the deposition profile, but all else being the same, this experiment demonstrates that clumping alone can completely change the spectral appearance of lower density ejecta (characteristic of higher energy explosions) in the nebular phase. In higher density ejecta, clumping is less influential because the density is higher and of low ionization even if the ejecta is smooth. Clumping has a greater impact only on smooth ejecta that are ionized. Cool recombined ejecta are less sensitive.

9.2 Evolution of the impact of clumping in a time sequence, from maximum light to late times

This section compares the evolution of two magnetar powered models that are identical, except that one is smooth (r0e2) and the other is clumped (r0e2cl). In practice, model r0e2cl starts off from model r0e2 at 29.6 d, but with a uniform volume filling factor of 10%. The time sequence for model r0e2cl is them computed as for r0e2, and until 400 d after explosion. This model shows how a fixed clumping level impacts the ejecta and the radiation as time proceeds.

Figure 26 shows the evolution of the photospheric properties for models r0e2 and r0e2cl. When they start off, the two models are identical, by design. But as the clumped model evolves, the enhanced density boosts the recombination efficiency of the gas. Initially, the effect is quite weak so the photosphere in both models have similar properties (although $T_{\rm ph}$ is markedly lower in the clumped model early on). But after 80 d, the ionization in the clumped model suddenly deviates from the smooth model, meaning that the gas recombines faster and the photosphere recedes in the ejecta (at earlier times, $T_{\rm ph}$ was dropping but not enough to cause a change in ionization). The entrance to the nebular phase occurs 30 d earlier in the clumped model as a result of the lower ionization of the gas. As shown here, oxygen, which is the dominant element in the ejecta, is recombined at 150 d but once ionized in the smooth model. Clumping has therefore a similar effect on the photosphere in this model as discussed for a Type II SN in Dessart et al. (2018). However, in contrast to the Type II SN ejecta, the bolometric light curve for models r0e2 and r0e2cl are identical to within a few percent (thus not shown). The reason is that in the Type II SN case, clumping immediately impacted the ionization and speeded up the photosphere recession (and the release of stored energy) whereas here, the choice of clumping has very little impact on the ionization up until 80 d after explosion. At that time, the ejecta has a small optical depth and the luminosity is close to the magnetar and decay power absorbed by the ejecta. A higher level of clumping would perhaps influence the ionization earlier on and impact the light curve, as obtained in Dessart et al. (2018). Quantifying the magnitude of clumping versus depth in magnetar powered SNe is left to future work.

The contrast in photospheric properties is reflected in the different spectral evolution of models r0e2 and r0e2cl (Fig. 27). The impact on the r0e2cl model spectrum is already visible at 30 d after explosion: the spectral features are essentially the same but the UV flux is reduced. As time proceeds, the SED becomes redder and eventually, the enhanced recombination is directly visible through the shift, for example, of O ii lines to O i lines. The bottom panel of Fig. 27 illustrates the line contributions for model r0e2cl at 63 d after explosion, with the preponderance of O i, Ca ii, or Fe ii rather than O ii and C ii in the smooth model r0e2 (see also Section 3). At nebular times, the results for r0e2cl are similar to those discussed in the previous section and are therefore not repeated.

Although the above results are function of the adopted clumping level and profile, they demonstrate that in a magnetar-powered SN, clumping may have a week effect up to maximum, but soon after that start to influence the ionization level and the spectrum appearance. When clumping is inferred from nebular phase spectra, it implies that clumping affects the SN radiation already soon after maximum because the ejecta optical depth is not large around maximum – the spectrum forms over the entire ejecta soon after maximum.

10 Comparison to observations

This section presents some comparison to observations. The next section first gives a preamble on the scope of these comparisons. Section 10.2 summarizes the observational data used. Section 10.3 presents a comparison of light curves and spectra to a few SLSNe Ic and in particular SN 2007bi, SN 2010gx, PTF12dam, PTF12gty, LSQ14an, SN 2015bn, and Gaia16apd.

10.1 Preamble

The radiative properties of a magnetar-powered SN are controlled by numerous processes such as the magnetar field strength and its initial spin, or the ejecta mass and its kinetic energy (Kasen & Bildsten, 2010). They also depend on how the magnetar power is deposited in the ejecta (Section 4; see also Dessart & Audit 2018), and whether the ejecta is clumped (Section 9; see also Jerkstrand et al. 2017). Even in its simplest form, there are a handful of parameters that impact the resulting SN photometric and spectroscopic evolution.

The present approach ignores dynamics and assumes a spherical ejecta and magnetar-power deposition. Given the results in Section 4, the model predictions on the rise to maximum (photometry, rise time, spectral evolution) are not robust. Similarly, little is known about the clumping, which implies that the color and spectral evolution after bolometric maximum is uncertain. At present, the treatment of clumping is 1D and neglects chemical segregation. The diversity in the properties of SLSNe Ic may come in part from asymmetry. Overall, it seems that the present simulations cannot set strong constraints on the progenitor stars (mass and composition) or how they exploded (energy). But the models give a robust qualitative description of the influence of a magnetar on an H-deficient and He-poor ejecta. Further improvements await a better knowledge of the 3D hydrodynamic properties of magnetar-powered SNe.

10.2 Dataset

There is now a vast collection of photometric and spectroscopic data for SLSNe Ic (see e.g., Quimby et al. 2011, Inserra et al. 2013, Nicholl et al. 2013, De Cia et al. 2018, Lunnan et al. 2018). The sample selected for the comparisons shown in the next section is limited to SN 2007bi, SN 2010gx, PTF12dam, PTF12gty, LSQ14an, SN2015bn, and Gaia16apd. Comparing to more SNe does not add anything substantial, since there is a lot of degeneracy in the observations and because of the limitations of the numerical approach. A summary of the relevant characteristics for the SN sample is given in Table 2. In practice, the photometric and spectroscopic data were retrieved from the SN catalog (Guillochon et al., 2017) at https://sne.space, except the maximum-light spectrum of Gaia16apd (Lin Yan, priv. comm.).

In addition, for SN 2007bi, the photometric data of Chen et al. (2015) is used to subtract the flux contamination from the host galaxy, as discussed in Jerkstrand et al. (2017).The inferred bolometric light curves of SN 2007bi and PTF12dam are from Chen et al. (2015). For LSQ14an, the nebular phase spectrum (already corrected for host contamination) is taken from WISEREP (Yaron & Gal-Yam, 2012).

10.3 Comparison of models to data

The top panel of Fig. 28 compares the inferred bolometric luminosity (with respect to the time of maximum) of a fraction of the model sample with the SLSNe Ic PTF12dam and SN 2007bi (Chen et al., 2015). This is done merely to show that the range of peak bolometric luminosities overlaps with that inferred from observations. However, in the model set, most of the light curves peak to a relatively faint maximum, exhibit a fast rise to maximum and a narrow light curve width. This results primarily from the large magnetic field used (1 to $3.5\times 10^{14}$ G) and the relatively small initial spins (most models have $P_{\rm ms}=$ 7.0 ms). Only the models r0e2ea ($P_{\rm ms}=$ 5.0 ms), r0e2eb ($P_{\rm ms}=$ 4.1 ms), and r0e2ec ($P_{\rm ms}=$ 2.0 ms) stretch to larger peak luminosities. In this paper, the goal was to understand the basic (non-dynamical) effects on SN ejecta and radiation for a given magnetar and investigate dependencies. The goal was not to match any specific SLSN Ic.

The bottom panel of Fig. 28 compares the $r$-band light curve (with respect to the time of $r$-band maximum) for a subset of models with a sample of observed SLSNe Ic. Models with faint bolometric luminosities exhibit a faint $r$-band maximum. But a faint $r$-band maximum is also obtained for the very energetic and luminous (in a bolometric sense) model r0e2ec. In this model, the bulk of the flux comes out in the UV, producing a relatively low brightness in the optical. The $r$-band or bolometric light curves can be modulated through changes in the magnetar properties. For example, PTF12dam is compatible with a faster spinning moderately magnetized magnetar, allowing a long rise time, broad light curve, and sustained brightness at late times (see e.g., Chen et al. 2015; De Cia et al. 2018; Nicholl et al. 2017). However, these works propose a factor of 5 range in magnetic field, a factor of $2-3$ in ejecta mass, as well as invoke a leakage of magnetar power at late times.

Figure 29 compares spectroscopic information for a few SLSNe Ic and a sample of our models. The top panel gives a comparison around maximum light (1 to 4 days prior to $r$-band maximum) using the observations of Gaia16apd, SN 2015bn, and PTF12dam. The UV and optical ranges are shown, since HST data is available for Gaia16apd and the UV range is always computed by cmfgen. The idea here is not to be quantitative. However, most of our model set, with the exception of the most energetic magnetars, predicts the O ii and C ii lines observed in SLSNe Ic (as revealed already in the first observations of these events; Pastorello et al. 2010; Quimby et al. 2011). Model r0e2ea is in good agreement with Gaia16apd from the UV to the optical, in the sense that the same line features appear in both spectra and with a location of the maximum absorptions at comparable wavelengths. As discussed earlier, the O ii lines are numerous and blended, while the C ii lines are only doublets and triplets and less affected by line overlap. The models tend to overestimate the observed strength of optical C ii lines, with the exception of SN 2018bsz (Anderson et al., 2018).

The bottom panel of Fig. 29 shows a spectral comparison at $1-2$ months after maximum. In models that assume a smooth ejecta, the spectrum is still very blue, with the persistence of O ii and C ii lines, as around the time of maximum. In contrast, the observations show a drastic evolution in spectral properties, with a dominance of lines from neutral and once-ionized species, in particular O i, Ca ii, and Fe ii. The only model in the present set that reproduces these features is r0e2cl, which assumes a clumped ejecta (see Fig. 27 for the full spectral evolution and the line identifications at 63 d after explosion). The disagreement with models is more severe for lower-density ejecta (i.e., higher $E_{\rm kin}/M_{\rm ej}$, as in model r0e4) or faster-spinning magnetars (i.e., higher $E_{\rm pm}$, as in model r0e2ec). In model r0e1, the $E_{\rm kin}/M_{\rm ej}$ is the lowest in our set and the ionization is lower than in model r0e2 after maximum, but not as low as in the clumped model r0e2cl. Here, clumping is not a tuning parameter for line strength. It is instead a fundamental tuning process for the ejecta ionization. Indeed, with clumping (here with a uniform volume filling factor of 10%), the ejecta ionization is low and suitable to produce most of the lines observed in the selected sample of observed SLSNe Ic (despite the diversity in their light curve properties, and hence in ejecta and magnetar properties). It thus appears that clumping is not just important at nebular times (Jerkstrand et al., 2017; Nicholl et al., 2018a). Clumping is also influential early after maximum (see Section 9.2). Physically, this makes sense since the total ejecta electron-scattering optical depth $\tau_{\rm es}$ is usually only a few tens at maximum. For example, in model r0e2, $\tau_{\rm es}$ is 65 at maximum and only 7 at about 40 d after maximum, so that the SN radiation already forms in the inner parts of the ejecta at that time. Thus, it is the same clumped inner ejecta that influences the SN radiation at nebular times and at early times after maximum. This clumping results from the dynamical influence of the magnetar, which operates most effectively early on, when the ejecta is compact, dense, and hot and the magnetar power is at its greatest. The clumpy structure of the ejecta is probably frozen within at most a few days after explosion. These results also suggest that SN 2007bi would have looked like, for example, PTF12dam, at maximum light, with a blue SED containing O ii and C ii, lines.

Figure 30 compares the observations of LSQ 14an and SN 2007bi with the models r0e2 and r0e2cl at about one year after maximum light. LSQ14an shows the presence of O iii 5006.8 Å (Inserra et al., 2017), while SN 2007bi appears similar to GRB SN 1998bw (Jerkstrand et al., 2017), with the presence of O i and Ca ii lines. Model r0e2 exhibits lines that are somewhat too broad (as a result of an adopted deposition profile that may be too broad) but all features present in LSQ14an are predicted by the model. The main contributions are from a forest of Fe ii lines up to about 5500 Å, followed by more isolated Fe ii lines at 6247.6 or 6456.4 Å, and then contributions from oxygen with O i at 7773.4 Å (the O i doublet at 6316.0 Å is a very small contribution to the Fe ii blend), O ii 7324.3 Å, and O iii 5006.8 Å. In LSQ14an, O i at 7773.4 Å is underestimated, primarily because oxygen is once ionized (see also Section 3). The O iii line is narrower and therefore forms over a smaller ejecta volume (smaller velocities) than in the model. Model r0e2 also predicts a strong Ca ii NIR triplet.

The clumped model r0e2cl predicts most of the features observed in SN 2007bi at 367 d after maximum. The contribution in the blue from Fe ii lines is somewhat overestimated, but important lines in the red part of the optical are well matched. The main contributions in the model at that time are O i 6300.0–6363.8 Å, Ca ii 7307.6 Å (with overlapping but weaker contributions from Fe ii and O ii), O i 7773.4 Å, and the Ca ii NIR triplet (at the edge of the CCD so the observed flux is probably uncertain). In the optical, magnesium lines are weaker than numerous overlapping contributions from Fe ii and Ti ii lines (for example, Mg i 4571.1 Å has the same EW as Co ii 4569.3 Å, and is about a tenth of the contribution from Fe ii 4549.5 Å  and a fifth from Ti ii 4572.0 Å). Jerkstrand et al. (2017) modeled the nebular phase spectrum of SLSNe Ic with line contributions exclusively from O and Mg, but the present model suggests that the spectrum may contain lines from other metals (this may depend on the assumption of microscopic mixing made here). In the observed spectra, there are narrow peaks corresponding to O iii lines but weaker than in LSQ14an. It is not possible to estimate from observations alone the contribution of O ii (if any) to the 7300 Å emission (since it overlaps with Ca ii), nor the contribution of Fe ii to the 6300 Å emission (since it overlaps with O i). The exact level of clumping depends on the density in the smooth ejecta, which is dependent on ejecta mass and kinetic energy, but clumping seems fundamental to reproducing the nebular-phase spectrum of SN 2007bi. This agrees with the findings of Jerkstrand et al. (2017), although the clumping of 0.1 (corresponding to a 10% volume filling factor) is sufficient here to cause an ionization shift in the present model set, while Jerkstrand et al. (2017) invokes a clumping 100 times greater. Physically, it seems extreme to fit about 90% of the ejecta mass in 0.1% of the ejecta volume, in particular if at the same time one assumes that the amount of magnetar power deposited in the ONeMg clumps is constant. Instead, as the clumps become smaller and smaller, one eventually enters a regime where little magnetar power gets absorbed by the clumps. Here, the magnetar power deposition is also assumed to be independent of clumping, but the clumping factor invoked is lower.

Magnetar power may lead to the co-existence of low-density hot ionized regions that would produce O ii and O iii lines, and dense, cool, recombined, clumps that would produce O i, Ca ii, and Fe ii lines. This seems to be the case in SN 2007bi, which exhibits weak O iii lines (and probably an O ii contribution to the 7300 Å). This co-existence may be complex. For example, this density and ionization structure may vary with radius and angle. A radial stratification could be easily discerned by a dichotomy in line width between lines from more ionized and less ionized ions. With the current treatment in cmfgen, the inter-clump medium is considered as void so that only the clumps can absorb or emit.

11 Summary and conclusions

This study presented non-LTE time-dependent radiative transfer simulations for magnetar-powered SNe. Using two solar-metallicity carbon-rich Wolf-Rayet progenitors and four ejecta in total (all endowed with ${}^{56}\rm{N}$i at the time of explosion), additional models are computed using a range of initial spins and magnetic field for the magnetar. The simulations are started at 1 d and continued until 200 to 600 d after explosion, depending on the model. The magnetar power deposition is prescribed using a simple, time-independent, analytical form function of the density and velocity structure. In each ejecta shell, magnetar power is treated the same way as radioactive decay, and the associated non-thermal effects are solved for.

Because the dynamical effects from the magnetar are ignored in cmfgen, the quantitative results are expected to be a little off, and the more so for faster-spinning magnetar models. For example, rise times or peak luminosities may be underestimated or overestimated, but the qualitative aspects and the trends should hold. The neglect of dynamics is quantified for models r0e2 and r0e2ec with the code heracles (see discussion and illustrations in Appendix A). The cmfgen and heracles results compare well for bolometric light curves (with slight offsets in rise times and light curve width) while the ejecta density structure is, as expected, most different in the case of the faster-spinning magnetar.

With magnetar fields in the range of 1.0 to $3.5\times 10^{14}$ G and initial rotational energies in the range of 0.4 to $5.0\times 10^{51}$ erg (spin periods in the range from 7.0 to 2.0 ms), the grid of magnetar-powered ejecta yield Type Ic SNe with rise times in the range of 21.6 to 51.7 d to a bolometric maximum in the range of 0.16 to $3.2\times 10^{44}$ erg s*-1*. The model counterpart without magnetar power peaks roughly at the same time, but with a luminosity ten times smaller and a narrower light curve. The main effect of the magnetar power is to raise the internal energy of the SN ejecta. This boosts the SN luminosity but also causes a drastic shift in temperature and ionization. In a standard SN Ic, oxygen is neutral or partially ionized at the photosphere at all times, while in a magnetar-powered SN Ic, oxygen is once-ionized throughout the high-brightness phase. Consequently, the ejecta optical depth is enhanced and the photosphere is located further out in radius (or velocity). The light curve is broader, the optical color is bluer, and the maximum-light optical spectra exhibit lines of O ii and C ii (instead of, for example, O i, Na i, and Fe ii in a standard SN Ic).

This temperature and ionization shift occurs in all the models presented here. In other words, even the weaker magnetar model can sizably influence the thermodynamics of the gas, alter the SN radiation, and produce the spectral appearance of a SLSN Ic. For the highest magnetar-energy model, the shift can be larger, with the O and C twice or three times ionized. This is possible in cases where the magnetar spin down timescale is shorter, allowing the temperature rise to occur while the ejecta is still relatively compact. Many different combinations of ejecta and magnetar properties are possible but the present grid of model reflects well the relative uniformity of observed properties of SLSNe Ic around maximum light.

Varying the extent of the magnetar-power deposition introduces shifts in the evolution of the SN radiation, while qualitatively, the SN model behaves in a similar fashion. With a narrower profile, the model takes longer to rise to a fainter maximum; the photospheric temperature and ionization increase later; the SED is at all times redder with a strong deficit of flux in the UV; the maximum-light optical spectra are essentially undistinguishable apart from the narrower line profiles; the late time spectra show lines from more ionized species. The differences at late times may not hold with a more realistic treatment of the magnetar power since then, the associated energy deposition may be extended in velocity space, and the more so as time progresses.

For a given ejecta mass and magnetar properties, increasing the ejecta kinetic energy (models r0e1, r0e2, and r0e4) by a factor of ten shifts the rise time from 43.0 to 25.8 d, and raises the peak luminosity from 2.53 to $4.86\times 10^{43}$ erg s*-1*. This occurs at constant time integrated luminosity since an increase in $E_{\rm kin}/M_{\rm ej}$ decreases the ejecta diffusion time (or more physically its ability to trap radiation). Around the time of maximum, the spectra are very similar, except that at higher $E_{\rm kin}$, the lines exhibit a greater velocity at maximum absorption and their emission component is weaker. However, at late times, the difference between models r0e1, r0e2, and r0e4 is large, reflecting the range in density in the inner ejecta for the different models. At higher $E_{\rm kin}$, the ejecta density is lower so that for a given magnetar power, the temperature and ionization are higher. As a result, in the nebular phase, model r0e4 exhibits lines of O ii and O iii while model r0e1 shows lines of O i, Ca ii, and Fe ii.

Varying the magnetar properties causes a change in the spin-down time scale (which goes as $P_{\rm pm}^{2}/B_{\rm pm}^{2}$), the total energy liberated (which goes as $1/P_{\rm pm}^{2}$), and the power at late times (which goes as $1/t^{2}B_{\rm pm}^{2}$). These introduce quantitative offsets but bring no new qualitative behavior. Reducing the spin period shortens the rise to maximum (from 31.7 to 21.6 from models r0e2 to r0e2eb) while lowering the magnetic field delays it. For sufficiently small spin periods, the SN can attain a very large temperature and ionization at the photosphere, making the SED bluer (with most of the flux falling in the UV) and strengthening O iii and C iii lines in the optical range. For a lower magnetic field, the power absorbed by the ejecta at late times is greater, boosting the ionization so that the ejecta cools mostly through forbidden-line emission from O ii and O iii transitions.

Non-thermal processes are explicitly treated in all simulations except one (model 5p11Bx2th) which is used to quantify the magnitude of these effects. The result shows that throughout the photospheric phase, the escaping radiation is essentially identical between model 5p11Bx2th and its non-thermal counterpart, that is, non-thermal processes have no visible impact on the plasma for this particular choice of ejecta and magnetar (in contradiction with Mazzali et al. 2016, but in agreement with expectations from previous studies, such as Kozma & Fransson 1992, Dessart et al. 2012a and Li et al. 2012). The interpretation is that the ionization is large (all the oxygen is generally once-ionized in all models), so that the free-electron density is also large, causing non-thermal electrons to lose their energy through Coulomb scattering, thereby heating the thermal bath. Since all SLSNe Ic exhibit lines of O ii during the high brightness phase, their photospheres are generally hot and ionized, and thus non-thermal effects should be small at those times. Instead, the high photospheric temperature and strong UV flux boost photoionization processes. At late times, non-thermal effects should strengthen if or when the ionization drops. In model 5p11Bx2, the He i 10830 line is indeed stronger than in the thermal model counterpart 5p11Bx2th at 100 d after explosion (early nebular phase). Overall, non-thermal effects thrive under low-ionization conditions, which are not representative of SLSNe Ic at maximum.

Following the approach presented in Dessart et al. (2018), a number of simulations are performed to test the effect of clumping on the SN radiation and gas properties. Clumping can impact the ionization level of the gas and all the above simulations for smooth ejecta (with the exception of model r0e1) exhibit a high level of ionization (with lines of O ii and O iii) until late times. There is also evidence from multi-dimensional simulations (Chen et al., 2016; Suzuki & Maeda, 2017) and nebular-phase spectra (Jerkstrand et al., 2017) that clumping is present. Here, the few simulations carried out show that the post-maximum bolometric light curve is unaffected but the UV-optical color of the clumped model is much redder. In model r0e2cl, clumping has a visible impact even at maximum light, although the greatest impact occurs a few weeks after maximum with a rapid drop in photospheric temperature and a rapid recombination of the ejecta to a partial ionization at +40 d. At this time, since the ejecta optical depth is only a few, the flux is roughly equal to the power absorbed and no longer scales with $R_{\rm ph}^{2}T_{\rm ph}^{4}$. At nebular times, the clumped model shows stronger O i and weaker O ii. Because the spectrum formation region encompasses the entire ejecta at 40 d after maximum, a clumped ejecta should exhibit distinct signatures from a smooth ejecta soon after maximum light – the effect of clumping is not confined to the nebular phase. Simulations at nebular times also show that as clumping is enhanced, the SN temperature and ionization drop as long as the conditions are ionized. Once partially ionized, a further increase in clumping has little impact. This is likely because the recombination rate drops due to the lower free-electron density as well as the strengthening of non-thermal effects at lower ionization. The effect of clumping is stronger in ejecta with a higher $E_{\rm kin}/M_{\rm ej}$.

The grid of models presented here was not designed to match any particular SLSN Ic. The choice of a relatively high magnetic field and slow spin implies that the models are located at the low-brightness and the fast-rise end of the distribution of potential SLSNe Ic. Nonetheless, the model properties reproduce the systematic presence of O ii and C ii lines in SLSNe Ic. Model r0e2 matches closely the maximum-light spectrum of Gaia16apd in the UV and optical range. The simulations suggest that clumping is essential to reproduce the colors and spectral lines observed within a month of maximum light. Without clumping, the models remain blue and persist in showing O ii and C ii lines for months. With clumping, a lower ionization is eventually enforced. For example, model r0e2cl matches well the observed spectra of the selected sample of SLSNe Ic, including SN 2007bi at +46 d or Gaia16apd at +43 d. The smooth ejecta model r0e2 reproduces some of the salient features of LSQ14an at +410 d (which exhibits a strong O iii 5006.8 Å line) and the clumped ejecta model r0e2cl reproduces well the observations of SN 2007bi at +367 d (characterized by lines primarily from neutral and once-ionized species). While clumping is supported by the present simulations, it is 100 times weaker than proposed by Jerkstrand et al. (2017). Here, the bulk of the effect of clumping is already present when adopting a 10% volume filling factor – enhancing the clumping further causes little change. In the future, this study will be extended to include a bigger set of progenitors covering from subsolar to solar metallicities, and also to improve the physical consistency with a better, less ad-hoc, treatment of magnetar-power deposition.

Acknowledgements.

LD thanks Lin Yan for providing the HST and optical maximum-light spectra of Gaia16apd. LD thanks ESO-Vitacura for their hospitality. This work utilized computing resources of the mesocentre SIGAMM, hosted by the Observatoire de la Côte d’Azur, Nice, France, as well as computing resources of the “Maison de la simulation”, CEA, Gif-sur-Yvette, France. This research was supported by the Munich Institute for Astro- and Particle Physics (MIAPP) of the DFG cluster of excellence “Origin and Structure of the Universe”.

Appendix A Comparison between cmfgen and heracles results

cmfgen is a non-LTE time-dependent radiative transfer code with no treatment of dynamics. The ejecta must be in homologous expansion and each mass shell moves at constant velocity. The density structure evolves as $1/t^{3}$. Hence, even in the presence of strong pressure gradients, as for example caused by energy deposition by a magnetar, there is no influence on the velocity or density structure. In reality, the magnetar influence on the ejecta may be strong, both on the internal energy budget (which may come out as radiation) and on the density structure, for example through a snow-plow effect (Kasen & Bildsten, 2010).

This section compares the results with cmfgen for models r0e2 and r0e2ec with the results from radiation-hydrodynamics simulations with heracles. The setup in heracles is the same as in Dessart & Audit (2018). The code starts from the cmfgen ejecta properties at 1.4 d. For simplicity, heracles uses a fixed composition with 99.9% oxygen and 0.1% iron when estimating the opacity of the material. A low floor opacity of 0.0003 cm2 g*-1* is used although it is in effect only at very low optical depth and thus has no influence on the light curve. The gray approximation is used. The code assumes an ideal gas equation of state, with $\gamma=$ 5/3, and a mean molecular weight of 10 (this is a rough estimate for an oxygen-dominated composition with partial ionization). The choice of molecular weight has no impact on the light curve (or dynamics; the plasma is radiation dominated) but it influences the gas temperature (where the albedo is high or where the optical depth is low). Radioactive decay is ignored in the heracles simulation.

Figure 31 compares the bolometric light curves computed with cmfgen and heracles for models r0e2 and r0e2ec. A time shift is applied to the heracles light curve to account for the light-travel time to the outer boundary at 1016 cm (which corresponds to 3.86 d). The light curves are in rough agreement. The main difference is in the early-time behavior and rise time, which arises from the different profile for the magnetar-power deposition. The light curve width differs also from the greater opacity returned by cmfgen  which accounts accurately for the effect of lines. At late times, the cmfgen luminosity is offset from the magnetar power absorbed for a number of reasons, including the additional contribution from 0.13 M*⊙* of ${}^{56}\rm{N}$i, the larger ejecta optical depth in cmfgen (which treats the opacity accurately), and also because the magnetar-power deposition is broader than in the heracles simulation (which caused a small boost in the observer’s frame during the photospheric phase).

Figure 32 compares the ejecta properties computed by cmfgen and heracles at 41.4 d after explosion. In the regions where the magnetar power does not affect the dynamics, the density structure is identical in both models, because it results exclusively from mass conservation. Deeper in the ejecta, the density structure differs in both models, and more strongly in the model r0e2ec characterized by a large magnetar energy ($E_{\rm pm}=5\times 10^{51}$ erg). As expected, the magnetar power has a strong dynamical effect (even with the relatively broad deposition profile used in heracles), and a snow-plow effect is visible through the presence of a dense shell at about 8000 km s*-1*. The snow-plow effect in model r0e2 is weak. The impact on the temperature (and the different results between the codes) result from the different dynamical evolution (in particular the resulting density), as well as the fact that the temperature is computed in non-LTE in cmfgen. The bottom panel of Fig. 32 shows that the velocity can deviate from homologous expansion for a large magnetar power and energy. For the model r0e2ec, the ejecta kinetic energy is roughly equal to the magnetar energy at birth ($E_{\rm pm}=5.0\times 10^{51}$ erg compared to $E_{\rm kin}=4.1\times 10^{51}$ erg).

Overall, despite the neglect of hydrodynamics in cmfgen, the resulting ejecta and radiation properties (for what can be compared) are similar to those produced by heracles. This suggests that the approach used for the main part of this paper is adequate. Furthermore, it is clear that the 1-D hydrodynamics treatment of magnetar-powered ejecta is not ideal since it overestimates the snow-plow effect and ignores multi-dimensional fluid instabilities, which are known to prevent the formation of a dense swept-up shell Chen et al. (2016); Suzuki & Maeda (2017). With the present approach, cmfgen allows for detailed non-LTE and time dependence and mimics the influence of some multi-dimensional effects. Further work is needed to improve the physical consistency of the approach.