On the Formation of the W-shaped O ii Lines in Spectra of Type I Superluminous Supernovae

H-poor Type I superluminous supernovae (SLSNe I) are characterized by O ii lines around 4000–4500 Å in pre-/near-maximum spectra, the so-called W-shaped O ii lines. As these lines are from relatively high excitation levels, they have been considered a sign of nonthermal processes, which may provide a hint of the power sources of SLSNe I. However, the conditions for these lines to appear have not been well understood. In this work, we systematically calculate synthetic spectra to reproduce the observed spectra of eight SLSNe I, parameterizing departure coefficients from the nebular approximation in the supernova ejecta (expressed as b neb). We find that most of the observed spectra can be well reproduced with b neb ≲ 10, which means that no strong departure is necessary for the formation of the W-shaped O ii lines. We also show that the appearance of the W-shaped O ii lines is sensitive to temperature; only spectra with temperatures in the range of T ∼ 14,000–16,000 K can produce the W-shaped O ii lines without large departures. Based on this, we constrain the nonthermal ionization rate near the photosphere. Our results suggest that spectral features of SLSNe I can give independent constraints on the power source through the nonthermal ionization rates.


INTRODUCTION
Superluminous supernovae (SLSNe) are known as a special type of supernovae (SNe) with extremely high luminosities.SNe with their peak absolute magnitudes of < ∼ −21 mag in optical bands are conventionally classified as SLSNe (e.g., Gal-Yam 2012; Moriya et al. 2018b;Gal-Yam 2019a;Nicholl 2021).Also, some SNe less luminous than ∼ −21 mag at their peak magnitudes are often classified as SLSNe by the similarity of the spectra to SLSNe (e.g.De Cia et al. 2018;Quimby et al. 2018;Gal-Yam 2019b).
The W-shaped O II lines are mainly produced by bound-bound transitions from the 3s 4 P (the spinquartet state) and 3s 2 P (the spin-doublet state) levels.The lower energy levels of these lines are ∼ 23 eV above the ground state (Figure 1).If the distribution across the excited states follows the Boltzmann distribution in local thermodynamic equilibrium (LTE), such high energy levels require high temperatures to be populated so that strong absorption lines can be produced.Therefore, non-thermal processes (in ionization and recombination and excitation) causing departure from LTE have been considered necessary to produce the W-shaped O II lines (e.g., Mazzali et al. 2016;Quimby et al. 2018).However, it is not yet clear whether non-thermal processes are required.In fact, Dessart (2019) suggested that the W-shaped O II lines can appear without the effect of non-thermal processes by performing radiative transfer simulations.
Departure from LTE is reminiscent of He I lines in spectra of He-rich SNe (SNe Ib).The prominent absorption lines of He I in SN Ib spectra arise from excited states with energies ∼ 20 eV above the ground state.Since this energy is much higher than the typical temperature of SNe Ib (T ∼ 5, 000 − 10, 000 K or ∼ 0.4 − 0.9 eV), the energy levels are not populated enough to produce the absorption lines in LTE.This suggests that they are populated by non-thermal processes.The departure from LTE can be understood well by γ−rays from 56 Ni decay (Lucy 1991;Hachinger et al. 2012).The excited states of He I are populated by high energy electrons produced by Compton scattering of γrays from 56 Ni decay.The high energy electrons ionize He I to He II, and then He II recombines to the excited states of He I (see Figure 1 of Tarumi et al. 2023).Lucy (1991) suggested that the degree of departure from LTE described by the departure coefficient b is ∼ 10 4 − 10 5 for the He I lines in spectra of SNe Ib (n = b × n * , where n and n * are the number density of atoms in the excited state in ejecta of actual SNe Ib and in LTE, respectively).
As the departure from LTE for the He I lines provides information on the power source of SNe Ib, the departure from LTE to produce the W-shaped O II lines in spectra of SLSNe-I may give us a hint of the power source of SLSNe-I.However, the detailed physical conditions where these lines appear are not understood well.In this work, we investigate the conditions for the appearance of the W-shaped O II lines by calculating synthetic spectra and comparing them with observed spectra.Then, we demonstrate that spectra of SLSNe-I can be used to give constraints on the power source.
This paper is organized as follows.Section 2 describes our samples of observed spectra of SLSNe-I.Section 3 gives an overview of the methods used for our calculations and parameters to calculate synthetic spectra.Section 4 shows results of the spectral calculations and comparison with the observed data.Section 5 discusses the behavior of the departure coefficients and implications on the power source of SLSNe-I.Finally, Section 6 summarizes the findings of this paper.

OBSERVATIONAL DATA
We use 66 spectra of eight SLSNe-I (see Table 1) taken from the WISeREP database (Yaron & Gal-Yam 2012) and corresponding photometric data from the Open Supernova Catalog (Guillochon et al. 2017, light curves are shown in Figure 2).The spectra were selected from spectra of the 28 SLSNe-I analyzed in Könyves-Tóth & Vinkó (2021) based on the following two criteria: i) availability of photometric data within three days of the spectra and ii) availability of spectroscopic or photometric data at least at one epoch covering wavelengths shorter than λ rest < ∼ 2,000 Å in the rest frame.The former condition enables accurate flux calibration of the spectra while the latter one secures an accurate determination of temperatures of SLSNe-I (typically T ∼ 10, 000 − 15, 000 K). Accurate temperature estimates are necessary to evaluate the departure coefficient for the appearance of the W-shaped O II lines since the LTE population strongly depends on temperature in the Boltzmann factor.For flux calibration, we scaled the fluxes of the observed spectra to the corresponding photometry.To obtain a photometric flux at each epoch when a spectrum was taken, we linearly interpolated two photometric data points in the magnitude space.Wavelength dependence in the flux scaling factor was considered as a linear function of wavelength if a spectrum has two or more corresponding photometric data.If a spectrum has only one corresponding photometric data, the whole spectrum was multiplied by a constant scaling factor.
Milky Way extinction was corrected for based on the dust map by Schlegel et al. (1998).Note that extinction within the host galaxy was not corrected for since its amount is often uncertain.We will discuss possible inaccuracies in Section 5.2.When the observed spectra are compared with synthetic spectra (described in Section 3), the distance modulus of each SLSN-I is applied to the synthetic spectra.As SNe in our sample are at moderate redshifts (Table 1), distance moduli are obtained from the redshifts with the cosmological parameters H 0 = 70 km s −1 Mpc −1 , Ω m = 0.3, and Ω λ = 0.7.

Code
For calculations of synthetic spectra, we utilize a onedimensional Monte Carlo radiative transfer code (Mazzali & Lucy 1993;Mazzali 2000).Here, we briefly describe the code by focusing on the assumptions relevant to this work.The code calculates the plasma condi- tions under the modified nebular approximation.The modified nebular approximation takes into account the fact that, in the SN atmosphere, the densities are low and radiative processes dominate, giving the largest influence on the level population (Abbott & Lucy 1985;Mazzali & Lucy 1993).Thus, the code does not assume LTE both for ionization and excitation.The modified nebular approximation is known to work well for spectra of normal SNe (e.g., Pauldrach et al. 1996;Sauer et al. 2006;Tanaka et al. 2008;Hachinger et al. 2012;Teffs et al. 2020).For more details of the code, we refer readers to Mazzali & Lucy (1993) and Mazzali (2000).
The code assumes a sharp photosphere inside the ejecta.At the photosphere, blackbody radiation is emitted.Photon packets replicating flux are then propagated through the SN ejecta outside the photosphere.For the photon propagation in the ejecta, electron scatterings and bound-bound absorptions are taken into account.The degree of scattering and absorption is determined by the level populations and ionization in the ejecta, which depend on density and temperature.A temperature structure in the ejecta outside the photosphere is established by tracing the photon packets, differentiating between a local radiation temperature (T R ) and an electron temperature (T e ).The local electron temperature is crudely assumed to be T e = 0.9 T R .The temperatures are estimated by an iterative process as the temperatures give level populations and ionization degrees, which in turn affect the energy flux, and thus the temperatures.
The level population (number density) of the j-th ionized element with atomic number i at the l-th excited level n i,j,l is (at some given radius): where g i,j,l is the statistical weight, E i,j,l is the excitation energy from the ground state, k is the Boltzmann constant, and T R is a radiation temperature.W in Equation ( 1) is the dilution factor calculated as where J is the frequency-integrated mean intensity from the simulation and B(T R ) is the frequency-integrated Planck function with the radiation temperature T R at some given radius.
Ionization is estimated by the modified nebular approximation (Mazzali & Lucy 1993): (2πm e kT R ) 3/2 h 3 e −Ii,j /kTR , (3) where n e is the number density of electrons, m e is the mass of the electron, h is the Planck constant, and I i,j is the ionization potential of the j-th ionized element with atomic number i.In Equation (3) η is defined as where δ is a correction factor for the optically thick region at wavelengths shorter than λ 0 = 1, 050 Å (the Ca II edge), and ζ is the fraction of recombinations going directly to the ground state.Typically, η ∼ 0.4 for O just outside the photosphere (W ∼ 0.5).Line scattering by bound-bound transitions from a lower level l to an upper level u in homologously expanding ejecta is treated in the Sobolev approximation (Sobolev 1957).A particularly important quality in this context is the Sobolev optical depth: where B lu and B ul are Einstein B-coefficients, f lu is the oscillator strength of the transition, t expl is the time since explosion, and λ lu is the wavelength of the transition.The effect of line fluorescence is also taken into account (Mazzali 2000).

Set up of calculations
The density structure of the ejecta outside the photosphere follows a power law (ρ ∝ r −n , where ρ is density and r is radius) with index n = 7.The radius can be expressed by a velocity v in homologously expanding ejecta (r = vt expl ).In our fiducial model (iPTF13ajg, see Section 4), the ejecta mass above v = 12, 000 km s −1 is M ej (v ≥ 12, 000 km s −1 ) = 3.7 M ⊙ .This ejecta mass is scaled by f ρ parameter as described below.Note that the total ejecta mass cannot be estimated from our modeling as the radiative transfer is solved only outside of the photosphere.Thus, in the following sections, we give the mass outside of a typical velocity, for which we adopt v = 12, 000 km s −1 , as a representative value.Abundances in the ejecta are assumed to be homogeneous, and set to be the same as those adopted by Mazzali et al. (2016) for the first spectrum of iPTF13ajg for all the 66 spectra analyzed in this work: X(He) = 0.10, X(C) = 0.40, X(O) = 0.475, X(Ne) = 0.02, Parameters of the spectral calculations are as follows: • f ρ : density scaling factor • v ph : velocity at the photosphere [km s −1 ] • t expl : time since explosion in the rest frame [days] • b neb : departure coefficient for the excited states of O II.
Note that the purpose of this work is to understand the condition for the W-shaped O II lines as compared with those of normal SNe.Thus, we define the departure coefficient b neb as a departure from the population expected from the modified nebular approximation (n OII = b neb × n OII,neb , where n OII and n OII,neb are the number density of O II in the excited state in ejecta of SLSNe-I and that in the modified nebular approximation, respectively).These parameters in our calculations to model the observed spectra can almost independently be determined from the following physical quantities.We constrain the density scaling parameter f ρ from the dilution factor so that the dilution factor W at v ph estimated by Equation ( 2) is converged to ∼ 1/2 (i.e., close to the value of the geometric dilution factor) for a given combination of the other parameters.f ρ is kept the same for the spectral series of each SN.L bol is constrained from the observed flux (with an accuracy of ∼ 10%).v ph is estimated from wavelengths of blueshifted absorption lines.t expl is then constrained using , where σ is the Stephan-Boltzmann constant and T eff is the effective temperature.We estimate t expl for each spectrum so that the entire spectral series for each SN are reproduced by the same explosion date.The estimated explosion date gives a rise time to the peak of the light curve (Table 1).The light curves in Figure 2 are shown as a function of the time since the estimated explosion date, and our estimated rise time is consistent with (or longer than) that estimated from an extrapolation of the light curve.Finally, b neb is estimated from depth of the W-shaped O II lines.iPTF13ajg: First, we show results of modeling for iPTF13ajg, which we use as our fiducial model as this object was also modeled in Mazzali et al. (2016).This object was reported in Vreeswijk et al. (2014).The ab- An example of modeling for the first spectrum of LSQ14mo (MJD = 56688).The black line shows the observed spectrum.The light blue line shows a synthetic spectrum with the parameters: L bol = 10 44.2 erg s −1 , v ph = 10, 500 km s −1 , t expl = 33 days, b neb = 1 (no departure), and fρ = 0.2.The estimated radiation temperature just outside the photosphere is TR ∼ 15, 000 K. The red line is a synthetic spectrum with the same parameters as those of the light blue line but with departure coefficient b neb = 10.The square points show photometric fluxes.Only the photometric point near 5,000 Å (the r-band) was used for flux calibration of the observed spectrum.The photometric points near 3,000 Å were not used because the coverage of the observed spectrum was not wide enough to be compared with the photometry.The crosses are integrated fluxes of the synthetic spectrum (b neb = 1) with the filters of UV where there is no coverage of the observed spectrum.
For the modeling, we adopt the same parameters as in Mazzali et al. (2016).By definition, the density scaling factor is f ρ = 1.0, which corresponds to the ejecta mass outside above v = 12, 000 km s −1 of M ej (v ≥ 12, 000 km s −1 ) = 3.7 M ⊙ .This object shows the W-shaped O II lines until t expl = 54 days.The required departure coefficients to reproduce the depths of the W-shaped O II lines are b neb ∼ 50.The depths of the W-shaped O II lines of this object require one of the  largest departure coefficients among all the objects in our sample.
PTF09atu: This object was reported in Quimby et al. (2011b).The absolute AB magnitude at u-band peak is ∼ −21.5 mag (Quimby et al. 2011b).To reproduce the time series of the spectra, the density scaling factor for this object is found to be f ρ = 0.5, which corresponds to the ejecta mass outside v = 12, 000 km s −1 of M ej (v ≥ 12, 000 km s −1 ) = 1.8 M ⊙ .This object shows the W-shaped O II lines at all the epochs at which spectra were taken (t expl = 23 − 46 days).The required departure coefficients to reproduce the depths of the Wshaped O II lines are b neb ∼ 1 − 50.
PTF09cnd: This object was also reported in Quimby et al. (2011b).The absolute AB magnitude at u-band peak is ∼ −22 mag (Quimby et al. 2011b).The spectra taken at MJD 55055 (2009 August 12) and MJD 55068 (2009 August 25) were also modeled in Mazzali et al. (2016).We updated flux calibration of the spectra using the newly published photometric data (De Cia et al. 2018); the fluxes are found to be higher than those in Mazzali et al. (2016) by a factor of ∼ 2.5.Thus, the bolometric luminosity of the best model in this work is also higher than Mazzali et al. (2016) by a factor of ∼ 2.5.With the increase of the bolometric luminosity, the time since explosion is longer than that in Mazzali et al. (2016) to have a larger photospheric radius.This updated time is also consistent with the light curve.The density scaling factor for this object is found to be f ρ = 0.7, which corresponds to the ejecta mass outside v = 12, 000 km s −1 of M ej (v ≥ 12, 000 km s −1 ) = 2.6 M ⊙ .Despited the differences in the luminosity and time, the overall properties of the model are consistent with those derived by Mazzali et al. (2016), e.g., the inferred density scaling factor (or mass outside of a certain velocity) for PTF09cnd is smaller than that for iPTF13ajg.
This object shows the W-shaped O II lines until t expl = 65 days.The required departure coefficients to reproduce the depths of the W-shaped O II lines are b neb ∼ 1 − 10.Note that the absorption line around 2,700 Å cannot be reproduced by our models because of the fixed homogeneous abundance.Mazzali et al. (2016) reproduces those lines by Mg II with the abundance X(Mg) ∼ 0.1, but we fixed the abundance to X(Mg) = 5 × 10  (Quimby et al. 2011b).For this object, the density scaling parameter is found to be f ρ = 3.0, which corresponds to the ejecta mass out-side v = 12, 000 km s −1 of M ej (v ≥ 12, 000 km s −1 ) = 11 M ⊙ .This object shows the W-shaped O II lines until t expl = 25 days.The required departure coefficients to reproduce the depths of the W-shaped O II lines are b neb ∼ 1 − 3.After t expl = 25 days, the W-shaped O II lines disappear.The observed spectra taken after MJD 55323 (t expl = 63 days) could not be reproduced maybe because of the fixed abundance and/or the single power density law and the mass assumption.Thus, we only show models for the observed spectra before MJD 55323.
SN 2011kg: This object was reported in Quimby et al. (2011a).The absolute AB magnitude at g-band peak is ∼ −20.8 mag (Inserra et al. 2013).For this object, the density scaling parameter is found to be f ρ = 0.1, which corresponds to the ejecta mass outside v = 12, 000 km s −1 of M ej (v ≥ 12, 000 km s −1 ) = 0.4 M ⊙ .This object shows the W-shaped O II lines until t expl = 37 days.The required departure coefficients to reproduce the depths of the W-shaped O II lines are b neb ∼ 1 − 10 (small departure).After t expl = 37 days, the W-shaped O II lines disappear.The absorption lines around ∼ 3, 000 Å and ∼ 3, 500 Å may be due to Fe II and/or Si II, which are not reproduced with our fixed abundance.However, this does not affect the departure coefficient required for the O II lines.
PTF12dam: This object was reported in Quimby et al. (2012).The absolute AB magnitude at g-band peak is ∼ −22 mag (Chen et al. 2017).For this object, the density scaling parameter is found to be f ρ = 0.4, which corresponds to the ejecta mass outside v = 12, 000 km s −1 of M ej (v ≥ 12, 000 km s −1 ) = 1.5 M ⊙ .This object shows the W-shaped O II lines until t expl = 59 days.The required departure coefficients to reproduce the depths of the W-shaped O II lines are up to b neb ∼ 50.After t expl = 59 days, the W-shaped O II lines disappear.
LSQ14mo: This object was reported in Nicholl et al. (2015).The absolute AB magnitude at g-band peak is ∼ −21.2 mag (Leloudas et al. 2015).The spectra of this object were also modeled in Chen et al. (2017).The first spectrum (at −7 days) is modeled with t expl = 18 days.However, the model spectrum shows a somewhat higher velocity of the W-shaped O II lines as compared with the observed spectrum.To have a better match in the velocity, we adopt the photospheric velocity smaller than that estimated in Chen et al. (2017) by ∼ 5, 000 km s −1 .Accordingly, the time since explosion in our model for this spectrum is estimated to be t expl = 33 days.This revised time is compatible with the light curve, and reproduces the later spectra as well.With this choice, the density scaling parameter for this object is found to be f ρ = 0.2, which corresponds to the ejecta mass outside v = 12, 000 km s −1 of M ej (v ≥ 12, 000 km s −1 ) = 0.7 M ⊙ .
This object shows the W-shaped O II lines until t expl = 54 days.The required departure coefficients to reproduce the depths of the W-shaped O II lines are b neb ∼ 1 (no departure).After t expl = 54 days, the W-shaped O II lines disappear.
SN 2015bn: This object was reported in Nicholl et al. (2016).The absolute AB magnitude at g-band peak is ∼ −22 mag (Nicholl et al. 2016).For this object, the density scaling parameter is found to be f ρ = 0.4, which corresponds to the ejecta mass outside v = 12, 000 km s −1 of M ej (v ≥ 12, 000 km s −1 ) = 1.5 M ⊙ .This object lacks the W-shaped O II lines throughout the epochs of available spectra (t expl = 63 − 135 days).This is the only object that does not show the W-shaped O II lines at any epochs in our sample.

Temperature and departure coefficient
In this section, we show the behavior of the temperature and the departure coefficient as a function of time.The left panel of Figure 5 shows radiation temperatures just outside the photospheres in the models with the best parameters as a function of time since explosion.The temperatures decrease with time.It is clear that the only spectra with temperatures higher than T R ∼ 12, 000 K show the W-shaped O II lines (see also Könyves-Tóth 2022).Dependence of the W-shaped O II lines on temperature is discussed in more detail in Section 5.1.The W-shaped O II lines can appear only up to t expl ∼ 60 days.In fact, more luminous objects tend to show the W-shaped O II lines somewhat longer (see the points with the black edges in Figure 2).This is because, for more luminous objects, higher temperatures can be achieved for a longer time.
The right panel of Figure 5 shows relations between time since the explosion and the departure coefficients.The W-shaped O II lines in many of the observed spectra in the early phases (t expl < ∼ 50 days) are reproduced well with the departure coefficients b neb ∼ 1.This means that the W-shaped O II lines are formed with little departure from the population in the nebular approximation.Some of the spectra with the W-shaped O II lines around the phases t expl ∼ 50 − 70 days require somewhat larger departure coefficients b neb ∼ 10 − 50.For the observed spectra without the W-shaped O II lines, we obtained upper limits of the departure coefficients.When the departure coefficients exceed the upper limits, the synthetic spectra would produce the W-shaped O II lines.The upper limits of the departure coefficients become larger at later times.
Temperatures and departure coefficients are plotted one against another in Figure 6.There is a tendency that spectra with higher (lower) temperatures require smaller (larger) departure coefficients.The spectra with the temperatures T R > ∼ 14, 000 K show the W-shaped O II lines with the departure coefficient b neb = 1 (no departure).As the temperatures decrease, larger departure coefficients are requisite for the spectra with the W-shaped O II lines.The departure coefficient for the W-shaped O II lines in the spectra of the SLSNe-I in our sample is b neb ∼ 50 at most.This is in contrast to the departure coefficient for the He I lines in spectra of SNe Ib (b ∼ 10 4 − 10 5 ; Lucy 1991).

Dependence of line strength on temperature
Here, we discuss the dependence of the strength of the W-shaped O II lines on temperature.Figure 7 shows a series of synthetic spectra with various temperatures.For the calculations of the synthetic spectra, the bolometric luminosity L bol was parameterized to change temperature.The model with the radiation temperature just outside the photosphere T R ∼ 15, 000 K shows the strongest W-shaped O II lines.Neither the models with temperatures T R < ∼ 13, 000 K nor with T R > ∼ 18, 000 K show the W-shaped O II lines.
To understand this behavior, we evaluate the strength of the absorption lines from the Sobolev optical depth as in Equation ( 5) under the nebular approximation (see Hatano et al. (1999) for a similar analysis for normal SNe under LTE).Absorption lines can appear when the Sobolev optical depth is τ lu > ∼ 1 above the photosphere over a range of velocities.A typical velocity range of the line forming region is v ∼ v ph − 1.5 v ph to reproduce the width of the absorption line seen in the O II lines.
The left panel of Figure 8 shows the fraction of O II at typical photospheric parameters evaluated with the same abundance adopted for the spectral synthesis calculations in Section 3.2.Near the photospheres, density is typically ρ = 3 × 10 −14 g cm −3 , which is almost constant over time.The right panel of Figure 8 shows the Sobolev optical depth τ lu of the W-shaped O II lines at the photosphere obtained from the ionization fraction and Equation (1) and ( 5).The Sobolev optical depth shown here is that of one of the most prominent O II lines, λ lu = 4, 649 Å.For the line, we apply f lu = 0.34 (Wiese 1996) and t expl = 40 days in Equation ( 5).
The Sobolev optical depth peaks at T R ∼ 15, 000 K with a value τ lu ∼ 4. As the temperature decreases from T R ∼ 15, 000 K, the Sobolev optical depth steeply decreases since the number density of the excited state decreases with the Boltzmann factor.On the other  hand, as the temperature increases from T R ∼ 15, 000 K, the Sobolev optical depth decreases because O II is ionized to O III.Therefore, a temperature around T R ∼ 15, 000 K makes the Sobolev optical depth peak.Only at T R ∼ 13, 000 − 18, 000 K (T R ∼ 14, 000 − 16, 000 K), the Sobolev optical depth exceeds τ lu ∼ 1 (τ lu ∼ 3).
This dependence of line strength on temperature can also be seen in Figure 6.The synthetic spectra with higher temperatures than T R ∼ 14, 000 K show the Wshaped O II lines with a departure coefficient b neb = 1.
On the other hand, the synthetic spectra with temperatures lower than T R ∼ 14, 000 K require departure (b neb > 1) to show the W-shaped O II lines.Most of the observed spectra with temperatures T R < ∼ 12, 000 K do not show the W-shaped O II lines.
Although SLSNe-I are sometimes classified into two types based on the absence or presence of the W-shaped O II lines in their spectra (e.g., Könyves-Tóth & Vinkó 2021), our results indicate that the difference between the two types would be only temperature.It is natural that all the spectra of SN 2015bn do not show the Wshaped O II lines because of its low temperatures.

Effects of host extinction on departure coefficients
We now examine the effect of host extinction on the departure coefficients via temperature.Although the temperatures estimated in the model are affected by host extinction, it was not corrected for as mentioned in Section 2. Here, we investigate a possible host extinction of iPTF13ajg, which requires one of the largest departure coefficients (b neb < ∼ 50) among our sample.The upper limit of the host extinction of iPTF13ajg is estimated to be E(B − V ) < 0.12 from the absence of Na I D lines (Vreeswijk et al. 2014).Thus, we applied E(B − V ) = 0.05 to the second spectrum of iPTF13ajg taken at MJD 56391 (one of the best spectra).
The host-extinction-corrected spectra and the models with the best parameters are shown in Figure 9.For the model of the spectrum with E(B − V ) = 0.0 (before the host extinction correction), the best parameters are f ρ = 1.0,L bol = 10 44.5 erg s −1 , v ph = 12, 250 km s −1 , t expl = 50 days, and b neb = 30.These parameters give a radiation temperature just outside the photosphere T R ∼ 13, 000 K.When we assume E(B − V ) = 0.05, the intrinsic spectrum becomes brighter and bluer.Then, the best parameters are f ρ = 1.0,L bol = 10 44.6 erg s −1 , v ph = 12, 250 km s −1 , t expl = 47 days, and b neb = 1.These parameters give a radiation temperature just outside the photosphere T R ∼ 14, 500 K.By this somewhat higher temperature, the model spectrum with b neb = 1 produces a deeper W-shaped O II absorption feature.
These results demonstrate that the required departure coefficients is quite sensitive to the host extinction, which is always difficult to estimate.Even a small host extinction correction largely affects the UV fluxes, increasing the estimated temperature.Thus, the required departure coefficient tends to be smaller when the host extinction is applied.

Constraints on ionization rate and implication to power sources
We here examine non-thermal processes in ejecta of SLSNe despite our finding that population of the excited states of O II does not largely deviate from population in the nebular approximation.This may yield constraints on the power source of SLSNe-I because population would be influenced by non-thermal excitation/ionization to some extent whenever there are γ-rays from either 56 Ni or a magnetar in ejecta.We test this for temperatures of T = 8, 000 − 12, 000 K, where the nebular approximation does not lead to formation of the W-shaped O II lines.
In analogy to He I lines in spectra of SNe Ib, the ionization rate can be constrained by the population of the excited states of O II.For the appearance of the He I lines, the excited states of He I are populated by high energy electrons that are produced by Compton scattering of γ-rays from 56 Ni decay (Lucy 1991;Hachinger et al. 2012).The high energy electrons ionize He I to He II with a rate determined by the energy deposition rate of γ-rays from 56 Ni decay.Then, He II recombines to the excited states of He I (see Figure 1 of Tarumi et al. 2023).Therefore, the population of the excited states of He I is related to the non-thermal ionization rate.While the presence of He I lines in spectra of SNe Ib gives the non-thermal ionization rate, the absence of the W-shaped O II lines in spectra of SLSNe-I can give an upper limit to the non-thermal ionization rate through population of the excited states of O II.
The number density of the excited states of O II should not exceed a certain value that produces the Wshaped O II lines.The typical value of the Sobolev optical depth producing absorption lines is of the order of τ lu ∼ 1.This corresponds to densities n i,j,l > ∼ 1 cm −3 for the W-shaped O II lines obtained from Equation 5by adopting λ lu = 4, 649 Å, f lu = 0.34 (Wiese 1996), and t expl = 40 days.This number density is discussed further below and shown to give an upper limit of the ionization rate.The ionization rate can then be converted to an energy deposition rate through the fraction of the deposition energy spent for ionization (among ionization, excitation, and heating of thermal electrons, see e.g., Fransson & Chevalier (1989)), which is called work  1) and ( 5), where f lu = 0.34 (Wiese 1996), t expl = 40 days, and λ lu = 4, 649 Å.The horizontal gray line in the right panel marks the Sobolev optical depth τ lu = 1, above which the W-shaped O II lines can appear.
per ion pair.Finally, this energy deposition rate can be related to power sources through the γ-ray spectra of power sources (and γ-ray transport).
In order to estimate the population of the excited states of O II under certain ionization and excitation conditions, we solve a simplified rate equation as illustrated in Figure 10.We only consider three levels: the ground state of O II, the ground state of O III (35.1 eV higher than the ground state of O II), and the excited O II spin-quartet state (23.0 eV higher than the ground state of O II).We do not consider the case of the excited state of O II in the spin-doublet state, which is expected to be similar to that in the spin-quartet state.This is because both the excited states would be balanced to be comparably occupied via O III.This is also seen in the excited state of He I in the spin-singlet state and that in the spin-triplet state achieved by the balance via He II (Lucy 1991).Also, we separately solve two equations below: an equation for the balance between the ground state of O II and the ground state of O III, and an equation for inflow to and outflow from the excited O II spin-quartet state.
We first consider balance between the ground state of O II and that of O III by non-thermal processes.The balance between the ionization from O II to O III and the recombination from O III to O II gives the number density of O III n OIII : where β is the non-thermal ionization rate (s −1 ), n OII is the number density of the ground state of O II, n e is the number density of electrons, and α OIII→OII (T ) is an effective recombination rate from the ground state of O III to the ground state of O II.Since we consider only non-thermal ionization here, the ionization rate β is expressed as where ė is the energy deposition rate per particle given by power sources (erg s −1 ), and w is the energy (erg) required to ionize an ion (work per ion pair).Here, w is usually obtained from the Spencer-Fano equation (Spencer & Fano 1954;Kozma & Fransson 1992).Treatment of the recombination processes is not trivial.The direct recombination to the ground state is often largely suppressed because it is immediately followed by absorption.In fact, in a typical density we consider, the optical depth of the recombination photons is quite high (τ > 10 6 ).In a realistic ejecta including various elements other than O, the recombination photons can also be absorbed by other ions, and thus, some direct recombination can still occur.Another path is recombination through the excited states.However, since we do not solve a full rate equation including the excited states of O II (see below), it is difficult to accurately evaluate the effective recombination rate through many excited states.Under these circumstances, we approximately adopt the direct recombination rate to the ground state from Nahar (1999) as an effective recombination rate α OIII→OII (T ) in Equation ( 6).This corresponds to the limit of the most efficient recombination to O II.In reality, the recombination rate would be lower than what is adopted here, and the number density of O III ions would be enhanced (implications of this assumption are discussed below).
Next, we solve the equation for inflow to and outflow from the excited O II spin-quartet state.This gives the

O III (Ground)
Ionization where σ PI (ν) is a photoionization cross section.The transition rate by electron collisions q if is computed as (Eissner et al. 1969), where g i is the statistical weight of an initial state, γ if is an effective collision strength from an initial state to a final state, and E i and E f are an energy level of an initial state and a final state, respectively.The escape probability β esc is computed as (Castor 1970) with the Sobolev optical depth τ lu in Equation ( 5).
The procedure above gives the number density of the excited stated of O II n OII,ex as a function of the ionization rate (Figure 11).Here, we assumed that all the ejecta consist of O, and that O II is dominant as shown in the left panel of Figure 8.This provides a condition of n e = n OII + 2n OIII .Also, we adopted the following parameters: ρ = 3 × 10 −14 g cm −3 as a typical density near the photosphere, α OIII→OII,ex (T ) from Nahar (1999), σ PI = 10 −17 cm 2 (Nahar 1999), g i = 4, γ if = 1, E i = 23.0 eV, E f = 0.0 eV for the collisional de-excitation, E f = 35.1 eV for the collisional ionization, A ex→gs = 9.8 × 10 8 s −1 (Nahar 2010), and W = 1/2; for τ lu , f lu = 4.3 × 10 −2 (Wiese 1996), n i,j,l = n OII , λ lu = 539 Å, and t expl = 40 days.Note that the photoionization cross section is assumed to be constant because the ionization cross section has many resonances just above the ionization energy.As the peak of the blackbody radiation with the temperature T = 8, 000 − 12, 000 K is located at much lower energy than the ionization threshold, the assumption of the constant cross section is effectively applied just above the ionization threshold.
Figure 11 shows the number densities for the excited O II spin-quartet near the photosphere as a function of the ionization rate.The number density of the excited O II spin-quartet exceeds the critical density n OII,ex ∼ 1 cm −3 (the Sobolev optical depth τ lu ∼ 1) when the ionization rate is β ∼ 3 × 10 −6 s −1 at T = 10, 000 K. At T = 8, 000 K and T = 12, 000 K, the number density for the excited state of O II exceeds the critical density when the ionization rate is β ∼ 10 −6 s −1 and β ∼ 3 × 10 −5 s −1 , respectively.
We note that the estimated number densities for the excited state of O II are lower-limit estimates because we assume blackbody radiation even at short wavelengths.In actual SNe, photons at short wavelengths are suppressed because of line blanketing.Smaller photoionization rates should increase the estimated number densities of the excited stated of O II.To show this effect, we compare the number densities for the excited O II spinquartet with photoionization rate to those without photoionization since the photoionization rate is uncertain because of line blanketing.The results at all the temperatures T = 8, 000 K, 10, 000 K, and 12, 000 K without photoionization are almost the same as the result at T = 8, 000 K with photoionization.At the temperatures T = 10, 000 K and 12, 000 K, the effect of the photoionization is the strongest to de-populate the excited state of O II in comparison to all the other processes (collisional de-excitation, collisional ionization, and spontaneous emission with self-absorption).In light of line blanketing, the result with radiation temperature T = 8, 000 K may best represent the actual condition in the ejecta of SLSNe.
The ionization rate can be translated to the energy deposition rate if the work per ion pair in Equation ( 7) is given.The work per ion pair becomes larger when the electron fraction is larger, since a larger fraction of the deposition energy is spent for heating of thermal electrons, not for ionization (Kozma & Fransson 1992).In electron-rich environments where ionization of metals produces many electrons, a typical value of the work per ion pair (calculated for Type Ia SNe) is w ∼ 30 I, where I is an ionization potential for an ion (Axelrod 1980).Hereafter, we assume a work per ion pair of w = 30 I for ionization of O II to O III in SLSNe-I.The deposition rates then corresponding to an ionization rate β that makes n OII,ex > ∼ 1 cm −3 are ė ∼ 2 × 10 −15 erg s −1 , ė ∼ 5 × 10 −15 erg s −1 , and ė ∼ 5 × 10 −14 erg s −1 at T = 8, 000 K, T = 10, 000 K, and T = 12, 000 K, respectively.
Finally, we demonstrate that the energy deposition rate can be related to power sources if the γ-ray spectrum emitted by the power sources and subsequent γ-ray transport are known.Here, to give crude constraints on the mass of 56 Ni in SLSNe-I, we hypothetically assume that 56 Ni is the power source for SLSNe-I.For simplicity, 56 Ni is assumed to uniformly supply the deposition energy to the ejecta (full mixing).The optical depth of γ-rays is τ γ > ∼ 1 in photospheric phases.Thus, for fully mixed 56 Ni distribution, the energy deposition rate given by 56 Ni decay can be approximately written as where L decay is the decay luminosity of 56 Ni (erg s −1 ) and N is the number of atoms in the ejecta.
where M Ni is the mass of 56 Ni.Here, we apply t expl = 40 days.The number of atoms is N A = M ej /⟨A⟩m p , where M ej is an ejecta mass, ⟨A⟩ is the average mass number in the ejecta (here, ⟨A⟩ = 16 for pure O), and m p is the mass of a proton.The ejecta mass is fixed to a typical value for SLSNe-I, M ej ∼ 6 M ⊙ (Blanchard et al. 2020).Then, the deposition rate needed for appropriate ionization rates can be converted to a mass ratio of 56 Ni vs the ejecta M Ni /M ej .This ratio is shown as to the upper x-axis in Figure 11.The mass ratio corresponding to the ionization rate that makes n OII,ex > ∼ 1 cm −3 is M Ni /M ej ∼ 0.05 at T = 10, 000 K. This is translated to 56 Ni mass of M Ni ∼ 0.3 M ⊙ for M ej = 6 M ⊙ .If SLSNe-I are entirely powered by 56 Ni, the required 56 Ni mass is > ∼ 3 M ⊙ , but it would lead to too strong a non-thermal ionization.
We emphasize that our estimates above involve a number of assumptions and uncertainties, and the upper limit of the 56 Ni mass is still quite uncertain.For example, the recombination rate adopted in the ionization balance (i.e., the direct recombination rate to the ground state) is expected to be lower because of immediate reabsorption.A lower recombination rate would increase the number density of O III ions and enhance the population of the excited state of O II, giving stronger W-shaped O II lines.Then, the upper limit of 56 Ni would become smaller (or tighter).On the other hand, there may also be the effects making the upper limit of the 56 Ni mass higher.For example, we only evaluate the Sobolev optical depth just outside the photosphere.But, in reality, the line is formed over a velocity range of up to v ∼ 1.5 v ph (Section 5.1).To have τ lu > 1 for a wider velocity region, more 56 Ni would probably be required, making the upper limit of the 56 Ni mass higher.Also, we assumed full mixing to translate ionization rate to the 56 Ni mass.Under less mixing, energy of non-thermal electrons is spent inside the photosphere, which decreases the non-thermal ionization rate outside the photosphere.Thus, with less mixing, the upper limit of the 56 Ni mass would also become higher (i.e., less stringent).
There are also other assumptions/uncertainties: the simplified rate equation, the uncertain photoionization rate because of the line blanketing, and the uncertain work per ion pair in ejecta of SLSNe-I.Nevertheless, our work demonstrates that spectroscopic properties can be in principle used to give constraints on the mass of 56 Ni in SLSNe-I, which roughly corresponds to the mass ratio M Ni /M ej of the order of 0.1.In addition, in the cases where SLSNe-I are powered by magnetars, spectroscopic properties can also be used to give constraints on spectra of γ-rays from magnetars via ionization rates if combined with detailed transfer calculations (Vurm & Metzger 2021;Murase et al. 2021).

SUMMARY
We have performed systematic spectral calculations to model the observed spectra of eight SLSNe-I to quantify the conditions for the formation of the W-shaped O II lines.We find that many of the pre-/near-maximum spectra with the W-shaped O II lines can be reproduced well with the departure coefficient b neb ∼ 1 (i.e., without departure) at the temperatures T R ∼ 14, 000 − 16, 000 K near the photosphere.This suggests that departure from nebular-approximation conditions is not necessarily large for the formation of the W-shaped O II lines in spectra of SLSNe-I.The appearance of the W-shaped O II lines is very sensitive to temperature.Thus, to understand the physical conditions for the line formation, it is important to estimate accurately the temperature of the ejecta from spectral modeling (rather than a simple fitting of the spectral energy distribution).We also highlight the importance of the extinction correction in the host galaxy; even a small extinction correction (E(B − V ) ∼ 0.05) can increase the intrinsic UV fluxes, which tends to increase the estimated temperatures by ∼ 2, 000 − 3, 000 K.
Finally, we have shown that the absence of the the W-shaped O II lines in spectra with a lower temperature (T < ∼ 12, 000 K) can be exploited to constrain the non-thermal ionization rate in the ejecta.Solving the simplified rate equation gives an upper limit to the nonthermal ionization rates.Under the several assumptions, this upper limit is roughly translated to an upper limit on the mass ratio M Ni /M ej of order 0.1.Similar methods can also be applied to give constraints on γray spectra of magnetars if detailed γ-ray transport is considered.Although our estimate involves a number of assumptions for simplification, our work demonstrates that spectroscopic properties can be used to give independent constraints on the power sources of SLSNe-I.

A. RESULTS OF SPECTRAL MODELING
The best parameters of the models for all the available observed spectra are summarized in Table A.

Quartet
Figure 1.A schematic energy diagram of O II mainly contributing the W-shaped O II lines.

Figure 2 .
Figure 2.Light curves of the eight SLSNe-I in our sample.We show light curves in the bands with the best temporal coverage.The difference in the colors correspond to the different SLSNe-I as shown in the legend.The points show the photometric observations.The points with the black edges show epochs when spectra with the W-shaped O II lines have been taken.The points with the gray edges show epochs when spectra without the W-shaped O II lines have been taken.

Figure 3
Figure 3 shows an example of comparison between synthetic spectra and an observed spectrum (the first spectrum of LSQ14mo).The observed spectrum is reproduced reasonably well with a departure coefficient b neb = 1 (no departure from the nebular approximation) at an estimated radiation temperature just outside the photosphere T R ∼ 15, 000 K. The model with b neb = 10 produces too deep W-shaped O II lines.The wavelengths of the W-shaped O II lines are matched well with a photospheric velocity v ph = 10, 500 km s −1 .The UV fluxes are strongly affected by metal absorption, which makes it difficult to estimate temperatures when observed spectral energy distributions are simply fitted by blackbody functions.Models with the best parameters for all the spectra in our sample are shown in Figure 4 and the best parameters are summarized in Table A (Appendix A).The properties of each SN are descried below.iPTF13ajg:First, we show results of modeling for iPTF13ajg, which we use as our fiducial model as this object was also modeled inMazzali et al. (2016).This object was reported inVreeswijk et al. (2014).The ab- Figure 3.An example of modeling for the first spectrum of LSQ14mo (MJD = 56688).The black line shows the observed spectrum.The light blue line shows a synthetic spectrum with the parameters: L bol = 10 44.2 erg s −1 , v ph = 10, 500 km s −1 , t expl = 33 days, b neb = 1 (no departure), and fρ = 0.2.The estimated radiation temperature just outside the photosphere is TR ∼ 15, 000 K. The red line is a synthetic spectrum with the same parameters as those of the light blue line but with departure coefficient b neb = 10.The square points show photometric fluxes.Only the photometric point near 5,000 Å (the r-band) was used for flux calibration of the observed spectrum.The photometric points near 3,000 Å were not used because the coverage of the observed spectrum was not wide enough to be compared with the photometry.The crosses are integrated fluxes of the synthetic spectrum (b neb = 1) with the filters of UV where there is no coverage of the observed spectrum.

Figure 4 .
Figure 4.All the results of modeling.Note that some spectra with other spectra taken at close epochs were excluded for visibility.For the observed spectra with the W-shaped O II lines reproduced well with b neb = 1, we only show one model with b neb = 1 (the light blue lines) for each observed spectrum.For the observed spectra with the W-shaped O II lines reproduced well with b neb > 1, we show two models: a model with b neb = 1 (the light blue lines) and a model with b neb > 1 (the red lines).For the observed spectra without the W-shaped O II lines, we also show two models: a model with b neb = 1 (the light blue lines) and a model with b neb > 1 (the orange lines) as an upper limit of b neb .The square points show UV photometric fluxes.The crosses show integrated fluxes of the synthetic spectra (b neb = 1) with the filters of UV.Vertical offsets for each spectrum are shown by the horizontal lines on the right y-axis.

Figure 5 .
Figure 5. Left panel: Relation between time since explosion and temperature of each spectrum.The different combinations of the colors and the shapes of the points relate to different SNe as shown in the legend.The filled and open points represent spectra with and without the W-shaped O II lines, respectively.Right panel: same as the left panel, but departure coefficient is shown as a function of time since explosion.Note that the open points show upper limits of the departure coefficients since the respective spectra lack the W-shaped O II lines.

Figure 6 .
Figure 6.Temperature vs departure coefficient.The data and the symbols are the same as in Figure 5.

Figure 7 .
Figure 7. Series of synthetic spectra with various bolometric luminosities L bol for a given combination of the other parameters to show the gradual influence of temperatures.Here, the bolometric luminosity L bol is parameterized with a step of log L bol = 0.1 erg s −1 from log L bol = 44.2erg s −1 (the red line) to log L bol = 45.0 erg s −1 (the dark blue line).Other parameters are fixed to v ph = 12,750 km s −1 , t expl = 40 days, and b neb = 1.Accordingly, the density scale factor fρ is varied from 1.3 to 0.6 so that the dilution factor converges to 1/2.The redder and bluer colors show models with lower temperatures (lower bolometric luminosities) and higher temperatures (higher bolometric luminosities).Radiation temperatures just outside photospheres are shown in the legend.The black line shows a synthetic spectrum with the deepest W-shaped O II lines (TR ∼ 15, 000 K).The gray shaded area highlights the W-shaped O II lines.

Figure 8 .
Figure 8. Fraction of O II among and the Sobolev optical depth at typical photospheres (W = 1/2 and ρ = 3 × 10 −14 g cm −3 ) with the same abundance as that adopted to the calculations of the synthetic spectra in Section 3.2.Left panel: dependence of the fraction of O II on temperature, estimated by Equation (3).Right panel: dependence of the Sobolev optical depth of one of the most prominent lines among the W-shaped O II lines estimated by Equation (1) and (5), where f lu = 0.34(Wiese 1996), t expl = 40 days, and λ lu = 4, 649 Å.The horizontal gray line in the right panel marks the Sobolev optical depth τ lu = 1, above which the W-shaped O II lines can appear.

Figure 9 .
Figure 9. Modeling of observed spectra of iPTF13ajg taken at MJD 56391 with the host extinction E(B −V ) = 0.0 (the lower spectrum) and E(B − V ) = 0.05 (the upper spectrum) corrected.The black lines show the observed spectra.The red line shows the model for the spectrum with the host extinction E(B − V ) = 0.0 with the best parameters: fρ = 1.0,L bol = 10 44.5 erg s −1 , v ph = 12, 250 km s −1 , t expl = 50 days, and b neb = 30.The model with b neb = 1 and the other parameters same as above is shown with the lower light blue line.The upper light blue line shows the model for the spectrum with the host extinction E(B − V ) = 0.05 with the best parameters: L bol = 10 44.6 erg s −1 t expl = 47 days, b neb = 1, and the other parameters are the same as above.The gray shaded area highlights the W-shaped O II lines.Vertical offsets for each spectrum are shown by the horizontal lines on the right y-axis.
Figure 10.A schematic diagram showing the populating mechanism of the excited O II spin-quartet state.Here, we consider 3s 4 P, which absorbs photons with 4,349 Å and 4,649 Å (see also Figure 1).

Figure 11 .
Figure 11.The number densities for the excited O II spin-quartet state nOII,ex as a function of the ionization rate (the lower x-axis).The upper x-axis shows the MNi/Mej ratio corresponding to the ionization rate.Different colors show different temperatures, as shown in the legend.Within this range of temperatures, O II is secured to be dominant as shown in the left panel of Figure 8 with typical density ρ ∼ 3 × 10 −14 g cm −3 at a photosphere.The number densities for the excited O II spin-quartet state, nOII,ex, in the nebular approximation at the temperature shown are too small to produce the W-shaped O II lines.The dashed lines and solid lines show the number densities with and without photoionization from the excited O II spin-quartet state.All solid lines (the cases without photoionization) overlap near the dashed line for T = 8, 000 K. The gray line represents nOII,ex = 1 cm −3 , beyond which the W-shaped OII lines can appear.

Table 1 .
Sample of SLSNe-I 1The rise time in the bands shown in Figure2.
number density of O II in that state n OII,ex :n OIII n e α OIII→OII,ex (T ) = n OII,ex (R PI +n e q OII,ex→OII + n e q OII,ex→OIII + β esc A ex→gs ), OIII→OII,ex (T ) is the recombination rate from the ground state of O III to the excited O II spin-quartet state, R PI is a photoionization rate, q OII,ex→OII is the collisional transition rate from the excited O II spinquartet state to the ground state of O II, q OII,ex→OIII is the collisional transition rate from the excited O II spin-quartet state to the ground state of O III, β esc is the Sobolev escape probability, and A ex→gs is the Einstein coefficient for the transition from the excited O II spin-quartet state to the ground state of O II.Note that the last term on the right side of Equation 8 includes both self-absorption and spontaneous emission.The photoionization rate R PI is computed as

Table A .
Summary of the best-fit parameters for each spectrum.NameMJD W-shaped O II lines f rho log L bol (erg s −1 ) v ph (km s −1 ) t expl (days)