Spectral sequences of Type Ia supernovae. II. Carbon as a diagnostic tool for explosion mechanisms

TL;DR

This study uses tomographic techniques to measure unburned carbon in a Type Ia supernova, providing constraints on explosion models and highlighting discrepancies with previous spectral identifications.

Contribution

It introduces a tomographic method to quantify carbon in supernova ejecta, offering a new approach to test explosion mechanisms against observational data.

Findings

Measured carbon mass fraction between 0.001 and 0.05 in specific velocity range.

Found an upper limit of 0.005 for carbon inside that region.

Identified discrepancies between observed features and model predictions.

Abstract

How an otherwise inert carbon-oxygen white dwarf can be made to explode as a Type Ia supernova remains unknown. A promising test of theoretical models is to constrain the distribution of material that is left unburned, in particular of carbon. So far, most investigations used line identification codes to detect carbon in the ejecta, a method that cannot be readily compared against model predictions because it requires assumed opacities and temperatures. Here, we instead use tomographic techniques to investigate the amount of carbon in the inner layers of SN~2011fe, starting from the previously published tomographic analysis of Mazzali et al. From the presence of the carbon feature in the optical at early epochs and its disappearance later on, we derive an average carbon mass fraction between 0.001 and 0.05 for velocities in the range $13500 \lesssim v \lesssim 16000\ \rm{km\ s^{-1}}$, and an upper limit of 0.005 inside that region. Based on our models and the assumed density profile, only small amounts of carbon should be in the neutral state, too little to be responsible for features seen in near-infrared spectra that were previously identified as due to neutral carbon; We discuss possible reasons for this discrepancy. We compare our results against a suite of explosion models, although uncertainties in both the models and our simulations make it difficult to draw definitive conclusions.

Figures

Tables

Table 1. Table 1: Summary of models.

Model	Scenario	References
W7	Fast deflagration	Nomoto et al. (1984)
N100	Delayed det.	Seitenzahl et al. (2013)
	transition	Röpke et al. (2012)
		Sim et al. (2013)
GCD200	Gravitationally	Seitenzahl et al. (2016)
	confined det.
(1.1, 0.9)$M_{\odot }\dots$	Violent Merger	Pakmor et al. (2012)
$1.0M_{\odot }$	Double det.	Fink et al. (2010)
$1.0M_{\odot }$	Spontaneous det.	Shen et al. (2018)

Table 2. Table 2: Relevant quantities.

	EpochbbRelative to time of explosion, with maximum light at 19.1 d. (d)
QuantityaaAll quantities for carbon are for a fiducial model where $(X_C^i,X_C^m)=(0.002,0.01)$.	3.7	5.9	9	12.1	16.1	19.1
Global properties
${\log }_{10}L/L_{\odot }$	7.903	8.505	9.041	9.362	9.505	9.544
$v_{\mathrm{inner}}\left(\mathrm{km}{s}^{-1}\right)$	13300	12400	11300	10700	9000	7850
Inner part ()
$M_{\mathrm{tot}}\left(M_{\odot }\right)$	$0.0205$	$0.0855$	$0.2127$	$0.2984$	$0.5132$	$0.6508$
$M_C\left(M_{\odot }\right)$	$0.0001$	$0.0003$	$0.0005$	$0.0007$	$0.0011$	$0.0014$
$M_{CI}/M_C$	$0.000006$	$0.000002$	$0.000000$	$0.000000$	$0.000000$	$0.000000$
$M_{C\mathrm{II}}/M_C$	$0.9780$	$0.8432$	$0.3692$	$0.2652$	$0.2153$	$0.2359$
max $\tau \left(\text{C I}\lambda 10693\right)$…	$0.29421$	$0.03836$	$0.00426$	$0.00111$	$0.00041$	$0.00022$
max $\tau \left(\text{C II}\lambda 6580\right)$	$16.25423$	$7.64200$	$1.15749$	$0.41434$	$0.18804$	$0.17041$
Middle part ()
$M_{\mathrm{tot}}\left(M_{\odot }\right)$	$0.0683$	$0.0683$	$0.0683$	$0.0682$	$0.0682$	$0.0682$
$M_C\left(M_{\odot }\right)$	$0.0007$	$0.0007$	$0.0007$	$0.0007$	$0.0007$	$0.0007$
$M_{CI}/M_C$	$0.000015$	$0.000005$	$0.000002$	$0.000001$	$0.000001$	$0.000001$
$M_{C\mathrm{II}}/M_C$	$0.9953$	$0.9929$	$0.9778$	$0.9751$	$0.9932$	$0.9966$
max $\tau \left(\text{C I}\lambda 10693\right)$…	$0.23596$	$0.02896$	$0.00318$	$0.00082$	$0.00032$	$0.00017$
max $\tau \left(\text{C II}\lambda 6580\right)$	$10.07514$	$1.99005$	$0.83938$	$0.30879$	$0.05106$	$0.01895$
Outer part ()
$M_{\mathrm{tot}}\left(M_{\odot }\right)$	$0.0316$	$0.0316$	$0.0316$	$0.0316$	$0.0316$	$0.0316$
$M_C\left(M_{\odot }\right)$	$0.0033$	$0.0033$	$0.0033$	$0.0033$	$0.0033$	$0.0033$
$M_{CI}/M_C$	$0.001179$	$0.000062$	$0.000009$	$0.000004$	$0.000003$	$0.000004$
$M_{C\mathrm{II}}/M_C$	$0.9988$	$0.9999$	$0.9997$	$0.9996$	$0.9999$	$0.9999$
max $\tau \left(\text{C I}\lambda 10693\right)$…	$0.22876$	$0.01900$	$0.00216$	$0.00056$	$0.00018$	$0.00009$
max $\tau \left(\text{C II}\lambda 6580\right)$	$0.31105$	$0.10694$	$0.04343$	$0.01573$	$0.00264$	$0.00101$
Outskirts ($v>19000\mathrm{km}{s}^{-1}$)
$M_{\mathrm{tot}}\left(M_{\odot }\right)$	$0.0095$	$0.0095$	$0.0095$	$0.0095$	$0.0095$	$0.0095$
$M_C\left(M_{\odot }\right)$	$0.0093$	$0.0093$	$0.0093$	$0.0093$	$0.0093$	$0.0093$
$M_{CI}/M_C$	$0.004206$	$0.000176$	$0.000012$	$0.000004$	$0.000003$	$0.000004$
$M_{C\mathrm{II}}/M_C$	$0.9958$	$0.9998$	$1.0000$	$1.0000$	$1.0000$	$1.0000$
max $\tau \left(\text{C I}\lambda 10693\right)$…	$2.56294$	$0.12246$	$0.01079$	$0.00248$	$0.00073$	$0.00037$
max $\tau \left(\text{C II}\lambda 6580\right)$	$0.00463$	$0.01429$	$0.01587$	$0.00880$	$0.00243$	$0.00088$

Full text

supernovae: general — supernovae: individual (SN 2011fe)

††software: Astropy (Astropy Collaboration et al., 2013), TARDIS (Kerzendorf & Sim, 2014; Kerzendorf et al., 2018)

1 Introduction

There is strong evidence that type Ia supernovae (SNe Ia) originate from explosions of white dwarfs (WDs) in binary systems. However, both the explosion mechanism and the stellar evolutionary path remain unknown (for a review see Livio & Mazzali 2018).

At least in part, the confusion surrounding the explosion mechanism arises from the fact that in a detonation the burning products are determined primarily by the local density of the fuel being consumed. For instance, carbon will remain mostly unburned at densities $\lesssim 10^{5}\ \rm{g\ cm^{-3}}$ (Shigeyama et al., 1992; Fink et al., 2010). In principle, models with different progenitor masses can still produce similarly plausible spectra provided that in the evolution up to the actual detonation the density profile of the star becomes similar. For instance, in the delayed detonation transition model (DDT, e.g. Khokhlov 1991; Blondin et al. 2013), a near-Chandrasekhar WD expands during the deflagration phase; the detonation that follows then propagates across lower densities, producing a distribution of elements similar to that of a pure detonation of a sub-Chandrasekhar WD (Sim et al., 2010).

In this series of papers we use the fast Monte Carlo radiative transfer code tardis (Kerzendorf & Sim, 2014) to try to determine whether the spectral features of goups of SNe Ia can be explained as variations of a single physical parameter. In paper I (Heringer et al., 2017), we showed that both normal and subluminous events can be reproduced using the same ejecta structure as long as the temperature was properly adjusted, suggesting that these events share a common explosion mechanism.

Here, we turn our attention to the distribution of carbon in the ejecta of SN 2011fe, a proto-typical normal event that was discovered possibly within a day of explosion, exhibited carbon and oxygen features in its pre-maximum spectra, and was extensively observed and analyzed (e.g., Nugent et al. 2011; Parrent et al. 2012; Vinkó et al. 2012). SN 2011fe reached its rest-frame $B$-band peak brightness on 2011 September 10 (MJD 55814.30 $\pm$ 0.06), corresponding to a time after explosion of $\sim$19 days (Mazzali et al., 2014, M14 hereafter). The exact rise time, however, will depend on how the 56Ni produced by the explosion is distributed in the ejecta–if there is little radioactive material near the surface, a “dark phase” may follow before any light can escape (e.g. Piro & Nakar 2013). M14 estimated this phase to last $\sim$1 day for SN 2011fe.

In principle, the abundance distribution of any given element could serve as a diagnostic tool for explosion mechanisms. However, carbon should be less susceptible to uncertainties pertaining to the precise nucleosynthesis, because the presence of unburned material is determined by a single burning process, compared to the multiple reactions required to form heavier elements. For instance, different groups have computed how much 56Ni is produced by the detonation of a WD of a given mass, finding discrepancies as large as a factor of 2, possibly due to differences in the nuclear reaction network adopted (see Shen et al. 2018 and references therein). Similarly, oxygen is less suitable because it is also a burning product of carbon (though in the nebular phase, which probes the deep interior, is very valuable: for instance, for SN 2010lp, an atypical SN Ia, it was used to infer an asymmetric explosion; Taubenberger et al. 2013).

The usefulness of carbon as a diagnostic tool arises because theoretical models predict the presence of unburned material at widely different regions of the ejecta. For instance, the violent merger model of Pakmor et al. (2012) leaves a large mass fraction ($\sim\!0.1$) of carbon mixed throughout the whole ejecta, while 3D simulations of the DDT scenario predict much more moderate amounts ($\sim\!0.01$) of carbon at velocities of $\sim 10000\ \rm{km\ s^{-1}}$ (Seitenzahl et al., 2013). In contrast, carbon might be concentrated in an intermediate layer in the double detonation scenario, as a result of the combustion of a helium shell (Fink et al., 2010; Sim et al., 2012). In models in which the supernova is approximated as the spontaneous detonation of a sub-Chandrasekhar WD with negligible helium mass, unburned material is expected only in the very outermost regions of the ejecta, at high velocities (Shen et al., 2018).

Wheeler et al. (1998) were among the first to investigate carbon in SNe Ia. They used a full radiative transport code (Hoeflich et al., 1995) to study SN 1986G and SN 1994D, finding the presence of Mg and O at similar layers in these events. Because Mg is a product of explosive burning of C, they were able to infer, based on delayed detonation models, that the transition between C and O burning takes place at $\sim 15000-16000\ \rm{km\ s^{-1}}$. Höflich et al. (2002) used the same code to study the subluminous SN 1999by, leading to the identification of a neutral carbon feature in the near-infrared (NIR) spectrum near maximum.

Subsequent research has used the parametrized spectral synthesis code SYNOW (Branch et al., 2002) to identify carbon signatures in the optical spectra of individual SNe Ia (SN 1998aq, Branch et al. 2003; SN 1999ac, Garavini et al. 2005 and SN 2006D, Thomas et al. 2007). More recently, the SYNAPPS package (an automated version of SYNOW, Thomas et al. 2011a) has been used to suggest the presence of unburned material in SN 2011fe at $\sim 10000-16000\ \rm{km\ s^{-1}}$ (Parrent et al., 2012).

Results from these semi-empirical studies include that the fraction of normal SNe Ia that exhibit a carbon feature if observed early enough is $\sim\!30\%$, and possibly higher if in some cases the carbon feature is blended with other features due to high Doppler shifts (Thomas et al., 2011b; Parrent et al., 2011; Folatelli et al., 2012; Silverman & Filippenko, 2012). Compared to normal SNe Ia, there appears to be a paucity of subluminous and overluminous events that show carbon (Parrent et al. 2011; see Chomiuk et al. 2016 for a short summary of distinct subtypes of SNe Ia). On the other hand, events that belong to the so called super-Chandrasekhar mass Ia subclass show a conspicuous carbon feature, significantly stronger than in normal events (e.g. Hachinger et al. 2012).

While parametrized packages, such as SYNOW, suffice to provide identifications and velocity estimates, they require temperature scales and reference opacities passed by hand and thus cannot be used to infer abundance ratios between elements, let alone quantitative abundance profiles across the ejecta. Detailed radiative transport codes (i.e. full NLTE or multi-dimensional) can be used to check abundance profiles, but their computational cost is too high to explore a large parameter space. In order to run the large number of simulations required to obtain abundance distributions of elements, one has to resort simplifications, and, e.g., use 1D codes with simplified excitation/ionization approximations.

Such simplified codes have been used for tomographic analyses, where the distribution of elements above the photosphere111Defined here as an inner boundary of the simulation domain where the radiation field is assumed to be that of a blackbody. is estimated by fitting the spectra at different epochs. As the ejecta expand and the photosphere recedes, deeper layers of the ejecta become optically thin, allowing the inference of the distribution of various elements. Here, we follow this approach to estimate the mass fraction of carbon relatively deep in the ejecta of SN 2011fe ($7850\lesssim v\lesssim 16000{\rm\,km\,s^{-1}}$) and compare it against model predictions. The study of the outermost layers could also potentially discriminate between models; However, for reasons we discuss in appendix A, we find that such analysis would require a more realistic tomographic analysis of the earliest spectrum, at $t_{\rm exp}$=3.7 d. Redoing the tomography of SN 2011fe is beyond the scope of this work and therefore we reserve the analysis of the high velocity layers to an appendix and limit ourselves to place only conservative constraints on the carbon mass fraction in that zone.

Following paper I (Heringer et al., 2017), as a starting point we use the tomographic analysis of SN 2011fe by M14, which led to density and abundance profiles with which the observed spectra during the photospheric phase (up to one week after maximum) can be reproduced, at least qualitatively (so far, no code can reproduce observed spectra well enough for meaningful quantitative assessments).

From their tomographic analysis of SN 2011fe, M14 suggest that carbon is present (mass fraction $\sim\!0.01$) down to $\sim\!8000\ \mathrm{km\ s^{-1}}$. This velocity is lower than predicted by many explosions models (see Section 3) and also lower than typically inferred for other SN Ia. For instance, Folatelli et al. (2012) and Silverman & Filippenko (2012) derived typical Doppler velocities (measured from the absorption minimum)222Note that the Doppler velocity of the red wing of the carbon feature in the optical may extend to $\sim 8000\ \mathrm{km\ s^{-1}}$, but this does not imply that unburned material is present at such velocities, as further discussed in §3.1. of $\sim\!12000\ \mathrm{km\ s^{-1}}$, while Parrent et al. (2011) used SYNOW to estimate velocities, obtaining similar values for the core-normal (Branch et al., 2006) SNe Ia in their sample. Inspired by this hint of carbon at depth and the scarcity of reliable measurements of the carbon profile in SNe Ia, we decided to re-assess the distribution of unburned material in SN 2011fe.

Focusing on a single element allows us to not just derive a best mass fraction, but a range of mass fractions that can describe the data, thus giving better hope of placing realistic constraints on the location and distribution of carbon that are useful for comparison against explosion models. It also helps to study more in depth the physical conditions that are relevant for the formation of the carbon feature (at least under the approximations we adopt), both in the optical and in the NIR, which may help to shed light on why the feature is more often seen in some subtypes than others.

This paper is divided as follows: in §2 we describe our methodology and in §3 we constrain the carbon mass fraction in the ejecta of SN 2011fe, which is then compared against model predictions in §4. In §5 we discuss the physical conditions that are relevant for the formation of the singly-ionized carbon feature in the optical, while in §6 we focus on neutral carbon and its possible impact on the NIR part of the spectrum. In §7 we discuss possible sources of uncertainties and the limitations of our analysis. Our conclusions are presented in §8.

2 Methods

To constrain the distribution of carbon in SN 2011fe, we use the spectral synthesis code TARDIS v1.5dev3181 (Kerzendorf & Sim, 2014). This code uses Monte Carlo radiative transfer through spherically symmetric ejecta. It treats the densest parts of the ejecta as an optically thick region that produces a blackbody distribution of photons. The photons are represented by packets which can interact with the layers above either via line interaction (in the Sobolev approximation) or electron scattering. TARDIS follows the approach and is implemented similarly to the Lucy (1999) code used by M14, except for possible differences in the atomic physics (see §3 of Kerzendorf & Sim 2014 for a detailed description), and generally produces consistent predicted spectra.

We chose TARDIS as a compromise between physical accuracy and computational cost in order to explore the parameter space of carbon mass fractions: opacities and temperatures across the ejecta are computed self-consistently, which is a clear improvement over line identification codes, but, as in the code used by M14, several approximations are made in order to reduce to computation time. The approximations that might affect our conclusions most are: (i) The ejecta are assumed to have no energy deposition and thermalization processes above the photosphere. (ii) The ejecta are treated as spherically symmetric, i.e., we simulate just the radial dimension333We usually take velocity as the radial coordinate, but this can readily be mapped to a radius at each given epoch since the ejecta are in homologous expansion at the timescales relevant for this work. and cannot capture any possible line of sight effects if the explosion was asymmetric. (iii) The distribution of ions and the population of atomic levels are computed using only a simple NTLE approximation, where these quantities are corrected by a “dilution” factor444We adopt the nebular and dilute-lte ionization and excitation modes in TARDIS, respectively.. We do not solve the full set of statistical equilibrium equations and thus ignore non-thermal processes (see Mazzali & Lucy 1993 and Kerzendorf & Sim 2014 for detailed descriptions.) (iv) Parameters other than the carbon abundance are fixed at the values found by M14, i.e., we do not attempt to see whether a change in the carbon feature due to a change in carbon abundance can be compensated by changes in the other parameters. Therefore, our results and conclusions are valid only within the scope of these approximations and comparison to models inherit these limitations.

2.1 Input parameters

Like in Paper I (Heringer et al., 2017), we start from the tomography of SN 2011fe by M14.555Since Paper I, we have further improved our model inputs by interpolating the abundances provided in Fig. 10 of M14, leading to slightly different synthetic spectra. This, however, does not affect the conclusions presented in paper I. We use the model parameters reported in their Table 6 and Fig. 10 and we generally find good agreement with the data (see Fig. 1 below). Small discrepancies between the synthetic spectra are likely due to minor differences in the codes (and secondary input parameters). The only significant difference we find is for the 3.7 days spectrum, near $\lambda=4000\,\text{\AA}$; this discrepancy is likely due to small abundance mismatches of Cr, Ti and/or Fe above the photosphere, located at $v_{\rm inner}=13300\ \mathrm{km\ s^{-1}}$, which can strongly influence the opacity at the blue end, even if the other parameters are similar (e.g., we infer a temperature of $10740$ K at the photosphere, which agrees with the value of $10800$ K reported by M14). Given their location, however, these abundance differences are unlikely to affect any of the conclusions we draw regarding carbon.

For the density profile, we use M14’s “$\rho$-11” profile. At the velocities relevant for our work ($v\lesssim 19000{\rm\,km\,s^{-1}}$), this is similar to other profiles that were tested, viz., W7 of Nomoto et al. (1984) and WDD1 of Iwamoto et al. (1999), but differs at higher velocities, where it is intermediate between W7 and WDD1, leading to a better match to the UV part of early spectra. We later show in §4 that, in the velocities of interest, this density profile is similar to those of the explosion models we compare against.

The times after explosion ($t_{\rm{exp}}$) for the spectra used in this work are listed in Table 2 below, together with the corresponding velocities of the inner boundary in our simulations ($v_{\rm{inner}}$, also referred to as “photospheric” velocity; these follow directly from M14). For comparison with our simulations, both the observed and simulated spectra were normalized by the mean flux computed in the region defined by $4000\,\text{\AA}\leq\lambda\leq 9000\,\text{\AA}$.

3 The Distribution of Carbon in SN 2011fe

In order to constrain the distribution of carbon in SN 2011fe, we first determine a velocity scale $v_{\rm{lim}}$ – and, thus, equivalently an enclosed mass – below which carbon is not required to explain the data. Next, we abandon the M14 carbon distribution and instead approximate it assuming constant $X_{\rm C}$ in two regions of interest: immediately below and above $v_{\rm{lim}}$, aiming to determine how much carbon could be “hidden” below the velocity threshold and what is the minimum amount of carbon required above it. To this purpose, we inspect the most prominent carbon feature in the optical spectra, due to the C II $\lambda$6580 line (e.g. Folatelli et al. 2012).

We note that in all cases, variations in the carbon mass fraction are made at the expense of the most abundant element in each layer so that the relative change of that element and the impact on the overall spectra are minimized.

3.1 The depth down to which carbon has to be present in SN 2011fe

We start our analysis by revisiting the distribution of carbon in the ejecta of SN 2011fe inferred by M14 (see the bottom panel of their Fig. 10). Their adopted distribution includes carbon at velocities as low as $\sim\!8000\ \mathrm{km\ s^{-1}}$, comparable to the velocity of the inner boundary at maximum light, $v_{\rm inner}=7850\ \mathrm{km\ s^{-1}}$, suggesting the presence of unburned material relatively deep in the ejecta.

Since the carbon feature appears stronger in the synthetic spectra of M14 than seen in the data (see their Fig. 8 and 9, near $\lambda\approx 6300$  Å), we first aim to determine quantitatively the velocity $v_{\rm{lim}}$ below which carbon is not required to explain the data. For this purpose, we perform a “depth scan,” where we adopt the carbon profile of M14 but remove all carbon (i.e., set $X_{\rm C}=0$) below a given velocity cut, $v_{\rm{cut}}$, which we vary in steps of $500\ \mathrm{km\ s^{-1}}$. $v_{\rm{lim}}$ is then determined as the largest $v_{\rm{cut}}$ for which the simulated spectra are consistent with the data.

Our simulations are shown in Fig. 1, where the top panels (a–f) correspond to different epochs and the bottom panel (g) shows the corresponding pseudo equivalent widths (pEW666The routine we employ here to compute pEW is made publicly available alongside TARDIS., e.g. Garavini et al. 2007) of the carbon feature for each choice of $v_{\rm{cut}}$. The models generally provide a good fit and one can see that a choice of $X_{\rm C}=0$ at $v\leq 13500\ \rm{km\,s^{-1}}$ suffices to explain the data. In fact, with the original carbon distribution below $13500{\rm\,km\,s^{-1}}$, for which $0.007\lesssim X_{\rm C}\lesssim 0.02$, the carbon features in the simulated spectra are stronger than those observed (while not providing a better fit anywhere else in the spectra). Based on the pEWs shown in panel g, we adopt $v_{\rm lim}=13500{\rm\,km\,s^{-1}}$.

We note that it is not clear why M14 included carbon at much lower velocities. They refer to Parrent et al. (2012) and possibly they were motivated by the fact that the spectra presented in Fig. 4 of Parrent et al. show the red wing of the C II $\lambda$6580 feature extending down to $\sim\!8000\ \mathrm{km\ s^{-1}}$ in Doppler velocity. The insets in Fig. 1 help to illustrate why the red-most extent cannot be used directly to determine where carbon-rich layers are located in physical space. For instance, the models for which $v_{\rm{cut}}=10500{\rm\,km\,s^{-1}}$ show absorption down to Doppler shifts as low as $-5000{\rm\,km\,s^{-1}}$, even though no carbon is actually present at those depths. This is because the line-of-sight component of the velocity of parts of the ejecta that scatter a line can be smaller than the true velocity, i.e., the absorption seen at low velocity can be due to parts of the ejecta that are at higher velocity but are not directly moving towards us.

3.2 Constraining the carbon profile

In order to constrain the range of carbon mass fractions that can reproduce the data, we next aim to determine for what $X_{\rm C}$ in a given velocity range the trough has a depth roughly consistent with the observations. Here, since the simulations are not yet at a level where they can reproduce spectra in detail, we refrain from employing any particular quantitative metric of consistency. Instead, we aim to determine upper and lower limits to the carbon mass fractions beyond which the simulated features are in obvious disagreement with the observed ones.

Since we compare with a feature that has only few parameters, it is clear that only limited information on the carbon profile can be gleaned. Indeed, we find that a simple two-zone model suffices, where we divide the velocity range into inner and middle regions, which are separated at $v_{\rm lim}=13500{\rm\,km\,s^{-1}}$, and which each have a constant mass fraction: $X_{\rm C}^{\rm i}$ at $7850\lesssim v\lesssim 13500{\rm\,km\,s^{-1}}$ and $X_{\rm C}^{\rm m}$ at $13500\lesssim v\lesssim 16000{\rm\,km\,s^{-1}}$. Here, the lower velocity limit of the inner region is set by the photospheric velocity of our latest spectrum, while the upper limit of the middle region is somewhat arbitrary. We explore the effects of moving the boundary between the inner and middle regions in §7, and discuss a trial in which we include an outer zone in appendix A.

In Fig. 2 we compare simulated spectra for which the carbon mass fraction is varied in the inner and middle regions with observed carbon profiles. For each epoch, we run simulations for combinations of $X_{\rm C}^{\rm i}$, $X_{\rm C}^{\rm m}$ in the range $0\leq X_{\rm C}^{\rm i}\leq 0.05$ and $0\leq X_{\rm C}^{\rm m}\leq 0.05$. This choice of range covers both extremes, resulting in spectra that either contain too little carbon for the feature to manifest or too much carbon, producing a flux depression much stronger than seen in the data. Note that for all but the last epoch, the photosphere is outside the lower velocity limit of our inner region. Hence, with increasing time, one constrains the mass fraction at increasing depth, thus justifying the choice of our inner region being coarser ($\Delta v\sim 5500{\rm\,km\,s^{-1}}$) than the middle region ($\Delta v\sim 2500{\rm\,km\,s^{-1}}$).

Inspecting Fig. 2, we infer for the middle region an upper limit of $X_{\rm C}^{\rm m}<0.05$: for larger values the carbon feature is stronger than observed for any choice of $X_{\rm C}^{\rm i}$ at $t_{\rm exp}=9\,$d. For the inner region, the carbon fraction is constrained more effectively by the later epochs, since more of the region is exposed. Hence, for $t_{\rm exp}=9$ days, we find $X_{\rm C}^{\rm i}<0.01$, while for later epochs the carbon feature is overpredicted unless $X_{\rm C}^{\rm i}<0.005$.

Setting a lower limit to $X_{\rm C}^{\rm m}$ is more complicated because a deficit of carbon in the middle region may be compensated by larger fractions of carbon in an outer region. From Fig. 2 alone, one would be tempted to set a lower limit of $X_{\rm C}^{\rm m}>0.005$ based on the strength of the carbon feature at $t_{\rm exp}=5.9\ \rm{d}$. However, we show in Fig. 8 in appendix A that if we allow the carbon fraction outside this region to be larger, somewhat lower values may work. Based on the spectra at $t_{\rm exp}=5.9\,$d, we infer a more conservative lower limit of $X_{\rm C}^{\rm m}\geq 0.001$.

Here, we note that inferences are hindered by the red wing of the Si feature not being predicted correctly, i.e. being offset by $\sim 50$  Å (the same can be seen in M14). As a consequence, the emission portion of this Si feature (i.e. a “limb brightening” effect [Hoeflich 1990]) will also be somewhat shifted, thus affecting the precise “pseudo-continuum” across the wavelength region relevant for the C II$\lambda$6580 line. Looking at the choice of $X_{\rm C}^{\rm i}\ =X_{\rm C}^{\rm m}\ =0$ in Fig 2, we expect a shift in the pseudo continuum to have a minor effect for a qualitatively analysis like ours, especially at $t_{\rm exp}\geq 9.0\,$d. Thus, the poor fit of the Si feature should not affect our analysis of the carbon fraction at greater depth, but will be important for any future assessment of carbon at high velocities.

In Fig. 3, we indicate schematically the parts of the carbon mass fraction parameter space that we believe are ruled out, labeling each region with the epoch based on which we excluded it. These are meant to be conservative limits. In particular, not all combinations of $(X_{\rm C}^{\rm i},X_{\rm C}^{\rm m})$ between their respective limits will produce spectra that match the data. Rather, the limits for one region indicate that there exists a choice of carbon fraction for the other region for which the data can be understood.

The above abundances are strictly for the adopted density profile. Since the optical depth scales, to first order, with density, the corresponding mass of carbon $M_{\rm C}$ is a less model-dependent quantity. For the middle region ($13500\lesssim v\lesssim 16000{\rm\,km\,s^{-1}}$), we find that $10^{-4}<M_{\rm C}<5\times 10^{-3}\,M_{\odot}$ is required to form the carbon feature. In the inner region ($7850\lesssim v\lesssim 13500{\rm\,km\,s^{-1}}$), we find an upper limit of $M_{\rm C}<2.5\times 10^{-3}\,M_{\odot}$.

Finally, we note that carbon is predominantly singly ionized in the middle region, so uncertainties in the degree of ionization should be small and our inferred range of $X_{\rm C}^{\rm m}$ reasonably reliable. For the inner region, a significant fraction of carbon can be doubly ionized, although the singly-ionized fraction never becomes small. In Table 2, we list inferred masses, ionization fractions and optical depths for reference values of $X_{\rm C}^{\rm i}=0.002$ and $X_{\rm C}^{\rm m}=0.01$.

4 Comparison to Explosion Models

We compare our constraints on the carbon mass fraction distribution with the predictions from different explosion models777Density and carbon abundance distributions for the explosion models were retrieved from HESMA (Kromer et al. 2017; https://hesma.h-its.org)., as summarized in Table 1. We show the comparison graphically in Fig. 4 and discuss the individual models below, but note here that, naturally, they predict not only distinct abundance profiles but also different density profiles, which are equally important in determining whether a feature manifests or not. Fortunately, as can be seen in Fig. 5, the density profiles predicted by different models are remarkably homogeneous: they agree within a factor of 5 for $v\lesssim 19000{\rm\,km\,s^{-1}}$ and are even closer in our regions of interest (i.e. $7850\lesssim v\lesssim 16000{\rm\,km\,s^{-1}}$). Indeed, we find that our conclusions below are unchanged if we compare predictions in terms of carbon mass (which is arguably more logical for line strengths, but we opted to show mass fraction instead as it more closely reflects the properties of the exploding white dwarf).

4.1 Fast Deflagration (W7)

The W7 model (Nomoto et al., 1984) was part of a suite of spherically symmetric pure deflagration models. For W7, the predictions are a good match to observations, but it required a flame propagation speed much higher than expected theoretically or found in more recent pure deflagration models (which also produce significantly different outcomes; e.g., Fink et al. 2014). Nonetheless, observables derived from the W7 model remain a standard comparison baseline, and thus we include it for reference.

From Figure 4, one sees that W7 predicts a sharp transition between regions where carbon is fully burned and not burned at all. This is consistent with no carbon in our inner region, but seems to overpredict the carbon in the middle region. Indeed, we find a carbon feature that is too broad and blue if we run a simulation with a W7-like carbon distribution.

4.2 Delayed Detonation (N100)

In delayed detonation models, run-away carbon fusion in a near-Chandrasekhar mass WD at first leads to a deflagration, which heats up and expands the white dwarf before transitioning to a detonation. These models are among the best studied and provide fairly good matches for observations.

For our comparison, we used data from the “N100” model (Seitenzahl et al., 2013), which was computed as part of a sequence of 3D simulations of delayed detonations, where the explosion strength is parametrized by the number of ignition kernels. The deflagration to detonation transition occurs when a pre-defined condition based on the local turbulent velocity is met (for a thorough description, see Ciaraldi-Schoolmann et al. 2013).

In this model, the input parameters are such that the predicted peak brightness and spectra resemble those of normal SNe Ia. The synthetic light curve for this model was computed by Sim et al. (2013) and compared to the normal SN 2005cf (Pastorello et al., 2007), exhibiting a reasonable agreement and only a modest dependency with the viewing angle, which is particularly small at the wavelengths relevant for the carbon feature. Some challenges to this model include the peak $B-V$ color being too red and the post maximum decline rate of the light curve being too fast when compared to typical values of normal SNe Ia.

In general, the DDT models of Seitenzahl et al. (2013) predicted optical colors that were systematically redder than expected at peak and did not reproduce the width-luminosity relation (Phillips, 1993) in its entirety. One significant difference of this model compared to 1D DDT models (where the flame transition to detonation is typically assumed to occur at a density threshold; e.g., Blondin et al. 2013) is that, because of the buoyancy of the hot ignition kernels, the deflagration burning products are not necessary centrally located at low velocities and might be mixed through zones rich in intermediate mass elements, such as Si.

We find that the carbon mass fraction predicted, on average, by this model is formally within the allowed region we derive here. We note, however, that this model has relatively little carbon in both the middle and outer regions. In Fig. 8 we see that $X_{\rm C}^{\rm m}\ \sim 0.001$ requires $X_{\rm C}^{\rm o}\ \gtrsim 0.2$ for a clear carbon trough to be formed at $t_{\rm exp}=$ 5.9 and 9 d, but N100 exhibits $X_{\rm C}^{\rm o}\ \lesssim 0.1$. We thus consider N100 to be a borderline case.

4.3 Gravitationally Confined Detonation (GCD200)

A gravitationally confined detonation (GCD) might occur if, once a hot deflagration bubble breaks out of the surface, some or most of the material were to remain bound. This could then wrap around the star, and trigger an off-center detonation upon convergence at the opposite side (Plewa et al., 2004). A general problem with these models is that they overpredict the abundances of iron-group elements at high velocities (Seitenzahl et al., 2016).

Nevertheless comparing our results to the simulations of Seitenzahl et al. (2016), we find that the carbon abundance distribution is roughly consistent, if perhaps a little too high, especially around $v\simeq 11000{\rm\,km\,s^{-1}}$. Given the approximations in our modeling, we also consider this to be a borderline case which we cannot rule out with confidence. We caution, though, that Seitenzahl et al. note that their model is too bright to be representative of normal SNe Ia.

4.4 Violent Merger (1.1M⊙+0.9M⊙)

Some SNe Ia might arise from the “violent” merger of two WDs, where the explosion is triggered at hot spots during the merging process. Existing studies for this scenario are not currently considered likely to be responsible for normal SNe Ia, as it predicts, e.g., far stronger polarization that is observed (Maund et al., 2013; Bulla et al., 2016).

We find that it also is inconsistent with our constraints on the carbon abundance profiles, producing, on average, too much unburned material mixed throughout the ejecta.

4.5 Double Detonation

In the double detonation scenario, a SN Ia arises when a detonation in a surface helium shell triggers a secondary detonation in the WD core (Livne, 1990). One of the outstanding questions for this class of models is whether the elements produced in the helium detonation would produce signatures, e.g., of iron-group elements, that are not observed (Kromer et al., 2010).

Fink et al. (2010) calculated element distributions for a range of progenitor masses, under the assumption that each had a helium layer with the minimum mass required to trigger a detonation. An interesting feature in these simulations is that unburned material is mostly confined to an intermediate layer – the carbon in the outer layers is partially burned by the initial shell detonation. The velocity of this intermediate layer increases systematically with the strength of the explosion, because for the more massive white dwarfs that lead to brighter explosions, a smaller minimum helium mass is required, which leads to unburned carbon further out in the ejecta.

For our comparison, we used model 3 of Fink et al. (2010) – which corresponds to the double detonation of a 1.0 M⊙ WD surrounded by a 0.055 M⊙ He shell – as it has the predicted 56Ni masses closest to the mass of $\sim\!0.45\,M_{\odot}$ estimated for SN 2011fe (Nugent et al., 2011). We find that the predicted distribution of carbon agrees reasonably well with the allowed region we infer in this work, especially considering that the predicted position of this layer may vary at the $1000{\rm\,km\,s^{-1}}$ level due to viewing angle effects.

4.6 Detonation of a Cold White Dwarf

Simulations of the detonation of cold WDs, i.e., with a detonation started in the centre without assigning a specific cause, have been shown to provide remarkably good matches to observations of SN Ia (Shigeyama et al., 1992; Sim et al., 2010). These models could be seen as examples of the double detonation model in the limit where the helium shell can be disregarded, perhaps representative of, e.g., models in which the helium is accreted very rapidly, at the onset of a merger, leading to helium masses of $\lesssim\!10^{-2}\,M_{\odot}$ (Shen & Moore, 2014).

For this class of models, the density gradient on the outside is steep and one thus expects carbon only at very high velocity. Indeed, in the recent simulations of detonations of cold WDs by Shen et al. (2018), carbon is virtually absent at the velocities we cover: $X_{\rm C}$ $<10^{-4}$ for $v<21000\ \rm{km\ s^{-1}}$. Shen et al. note, however, that a proper simulation of a double detonation would have to take into account the inward shock launched by the helium shell detonation, which might tamp the explosion and reduce the velocities of the outer layers.

4.7 Summary of the Model Comparisons

The carbon distributions predicted by the violent merger model of Pakmor et al. (2012) and by the detonation of cold WDs model of Shen et al. (2018) both seem inconsistent with the constraints on the carbon mass fractions (and masses) we derived for SN 2011fe, with the former predicting, on average, too much and the latter too little unburned material at intermediate velocities (see Fig. 4). We suggest that these differences are large enough for these models to be ruled out, despite the approximations involved in computing our synthetic spectra.

The carbon distributions predicted by the gravitationally confined detonation model of Seitenzahl et al. (2016) and by the N100 delayed detonation model of Seitenzahl et al. (2013) are roughly consistent with our results. This might provide a clue to what physical processes are needed to explain carbon distribution in the ejecta: in both mechanisms, the initial deflagration leads to relatively low density material at the locations that end up at the intermediate velocities where we infer carbon is present.

If the carbon distribution is indeed a consequence of the deflagration phase, then one might expect it to be somewhat stochastic and dependent on viewing geometry, thus providing a possible explanation for the variation in strengths of the carbon features in SN Ia spectra (Parrent et al., 2011). This hypothesis, however, needs to be further investigated. For instance, Sim et al. (2013) showed that the light curves and spectra derived for the N100 model had only a mild dependency on viewing angle, despite significant anisotropies in the flame propagation. A similar conclusion was reached by Fink et al. (2018), who computed their observables in the context of DDTs in differentially rotating WDs.

Finally, the predicted carbon distribution for the double detonation model is also consistent with our results and thus cannot be ruled out.

5 The Physical Conditions Relevant for the Carbon Feature

In this section, we attempt to answer the following questions: Where in the ejecta is the C feature formed?; Why does the C feature usually disappear near maximum light? and How much carbon can be hidden in normal SNe Ia that do not exhibit a C feature? We will do this in the context of our models; uncertainties in those, including the NLTE approximations made in TARDIS, are addressed in Section 7. Unless otherwise stated, we will compare with a fiducial model that reproduces the observed features of SN 2011fe well, with $X_{\rm C}^{\rm i}=0.002$, $X_{\rm C}^{\rm m}=0.01$ (and the carbon profile of M14 outside).

The carbon $\lambda 6580$ feature arises from the $\rm{2s^{2}3s}~{}^{2}{\rm S}\rightarrow\rm{2s^{2}3p}~{}^{2}{\rm P}$ transition in C II (Kramida et al., 2018). Because the ejecta are expanding, one expects a P-cygni profile. This will be superposed on the red wing of the main Si $\lambda 6355$ feature, which means that, depending on the latter’s strength, one may note the presence of carbon either through an actual trough or just through flattening of the silicon feature (as can be seen in Fig. 2 for the different choices of $X_{\rm C}^{\rm m}$ at, e.g., $t_{\rm exp}=5.9\,$d and $X_{\rm C}^{\rm i}=0$).

5.1 Where in the Ejecta is the Carbon Feature Formed?

To trace where a given feature is formed, we inspect the ejecta’s relevant line opacities. In TARDIS, these are calculated in the Sobolev approximation, in which interaction is considered only when a photon packet is Doppler shifted into the line (as appropriate for SN Ia ejecta with their large velocity gradients [Sim 2017]), and the code allows one to retrieve the Sobolev optical depth (i.e., the opacity integrated over the velocity width).

In the top panel of Fig. 6 we show the Sobolev opacities of the C II $\lambda$6580 line for our fiducial model. One sees that the opacity peaks in a relatively small velocity range, between just above the pseudo-photosphere and $\sim\!14000{\rm\,km\,s^{-1}}$ (consistent with the findings of Parrent et al. 2012).

5.2 Why does the C feature usually disappear near maximum light?

The line opacity depends on the column density of the relevant ions and will thus decrease with time as $t^{-2}$ due to the expansion of the ejecta if the level distributions do not change. Of course, these do also change and furthermore the amount of material involved in generating the line increases as the pseudo-photosphere recedes.

In the bottom panel of Fig. 6, we show how the fraction of ions at the relevant level for forming the carbon feature depends on the local density and temperature (under the NLTE approximations adopted here), with the temperature-density profiles of SN 2011fe at different epochs overdrawn. One sees that for typical conditions, the hotter and denser regions are favored for forming the carbon feature; while in the hottest layers a significant fraction of carbon might be more than singly ionized (see Table 2), this is compensated by increased population of the relatively high excitation state from which the carbon feature arises.

One sees that near maximum, the inner regions of the ejecta are still sufficiently hot and dense that a significant fraction of carbon is in the relevant excitation state for the formation of the feature, but nevertheless the opacities have become quite small (see top panel of Fig. 6 and Table 2). This indicates that the primary factor why the carbon feature weakens is the dilution of the column density due to the expansion of the ejecta. Nevertheless, at $t_{\rm{exp}}=16.1$ and 19.1 d, the peak opacity is still of the order of 0.2, indicating that even though no carbon trough is formed, carbon still contributes to the flattening of the main silicon feature.

An important clue is provided by the fact that a carbon feature is almost never observed near maximum in normal SNe Ia, even though the density and temperature near the photosphere are typically ideal for forming the carbon signature. This implies that the inner regions of the ejecta must generally be carbon poor, with $X_{\rm C}<5\times 10^{-3}$ inferred from Fig. 2 (where $X_{\rm C}^{\rm i}=5\times 10^{-3}$ does lead to the formation of a feature at $t_{\rm{exp}}=16.1$ and $19.1\,$d).

In this respect, an apparent exception to the rule may be informative: SN 2002fk did exhibit a carbon feature even one week past maximum (Cartier et al., 2014). This SN Ia was relatively normal in its lightcurve, but stood out in that it exhibited a small expansion velocity ($v\lesssim 10000{\rm\,km\,s^{-1}}$, as measured near and post maximum from the Si II $\lambda$6355 line). Using near-maximum velocities for a crude estimate, SN 2002fk’s ejecta was slower by a factor of $\sim 0.94$ and thus, if conditions were otherwise similar to those of SN 2011fe, the column density in the former would be increased, and the time at which the carbon feature would disappear lengthened, by the inverse of that factor squared. This correction alone cannot account for carbon being observed at such late epochs in SN 2002fk, which suggests either even lower pre-maximum velocities or different conditions relatively to SN 2011fe.

5.3 How much carbon can be hidden in normal SNe Ia that do not exhibit a C feature?

For objects similar to SN 2011fe, Fig. 2 indicates that towards maximum light, for inner carbon mass fractions $X_{\rm C}^{\rm i}\lesssim 10^{-3}$ there is little impact on the spectra, while for $10^{-3}\lesssim X_{\rm C}^{\rm i}\lesssim 2\ \times 10^{-3}$ the red wing of the main silicon feature is flattened, and for $X_{\rm C}^{\rm i}\gtrsim 5\ \times 10^{-3}$ there should be a detectable carbon feature. As expected, these ranges correspond to a transition from negligible to substantial optical depth ($\tau\lesssim 0.09$, $0.09\lesssim\tau\lesssim 0.18$, and $\tau\gtrsim 0.46$, resp.).

5.4 Overluminous versus “super-Chandrasekhar mass” events

For the typical conditions in SN Ia ejecta, one generally expects that hotter and denser layers are more favorable for the formation of the carbon feature (see Fig. 6, bottom panel). Given this, what can one learn for other sub-types of SN Ia? The carbon feature is rarely observed in overluminous (91T-like) events (Parrent et al., 2011). These are thought to have somewhat higher temperatures than normal SN Ia like SN 2011fe, but otherwise similar density and velocity structures, which suggests the carbon feature should be, if anything, easier to detect: even if carbon were to become doubly ionized in the innermost layers, the conditions would still be favorable for forming the carbon feature in the layers above, which are generally more carbon rich. The lack of observed carbon signatures would then suggest that in these events the amount of unburned material is truly lower than in their normal counterparts.

On the other hand, for super-Chandrasekhar mass candidates, which are also associated with relatively hot ejecta, the carbon feature is one of the defining features. This would be unexpected if these objects were part of the same family: since they arise from even more energetic explosion than the 91T-like SN Ia, one would naively expect them to contain even less unburned material. It is less clear whether the presence of carbon indicates a much larger abundance of unburned material, because the velocities inferred from the carbon feature in super-Chandrasekhar mass candidates are relatively low ($7500\lesssim v\lesssim 10000{\rm\,km\,s^{-1}}$; Taubenberger 2017), and thus the carbon column density will be diluted relatively slowly (as discussed for SN 2002fk in Sect. 5.2).

Independent of the precise reasons, the above suggests a distinct explosion mechanism. It might be particularly fruitful to study overluminous and super-Chandrasekhar events together with possible “transitional” events like SN iPTF13asv (SN 2013cv, Cao et al. 2016).

6 Neutral Carbon

So far, we have only considered singly-ionized carbon, and in particular its $\lambda 6580$ transition. However, Hsiao et al. (2013) presented evidence also for the presence of neutral carbon in SN 2011fe, from a flux deficit in the NIR attributed to C I $\lambda 10693$ (which has the strongest oscillator strength in the near-infrared region of the spectra; Marion et al. 2006), a flux deficit that seems to become more relevant at epochs near maximum (see their Fig. 4), and appears to move to lower velocity.

This detection is puzzling, since in our models carbon is mostly singly ionized throughout the ejecta (consistent with early work; e.g., Tanaka et al. 2008). Indeed, we find that neutral carbon does not contribute significantly to our synthetic spectra, except possibly at very early times, $t_{\rm{exp}}=3.7$ and 5.9 d, where the optical depth in the C I $\lambda 10693$ feature reaches $2.6$ and $0.12$, respectively (see Table 2). But this happens only in the outermost layers, above $19000\ \rm{km\ s^{-1}}$, which, with $X_{\rm C}=0.97$, consist of almost pure carbon888M14 note that a C/O abundance ratio closer to 50/50 would worsen the fit to the oxygen features.. The line opacity decreases drastically at later epochs; e.g., $\tau<0.02$ for $t_{\rm{exp}}\geq 9\,$d and, consistently, we do not see clear signatures of C I $\lambda 10693$ in these simulated spectra.

We note that the fits are not particularly good in the NIR region, e.g., not reproducing the neighbouring Mg feature near $1.05\,\mu$m (see Fig. 1). This is possibly because the approximation of a sharp photosphere is less accurate to simulate the NIR region, implying that long-wavelength radiation from below the assumed photosphere could also interact with the ejecta above. Importantly, in our models, we find neutral carbon to have higher opacities in the cooler layers, which are located above the photosphere; Therefore, the approximation of a sharp photosphere is unlikely to affect the opacities we derive for neutral carbon. Also, we note that changing the abundance of Mg is unlikely to change the fact that carbon is predominantly ionized, rather than neutral, in our simulations. To verify this, we ran an additional model at $t_{\rm{exp}}=9.0$  d where the Mg abundance is increased by a factor of 10 in the carbon rich region; in this model, the opacity of the C I $\lambda 10693$ line remained small.

The above raises the question of whether neutral carbon is relevant or not in SN 2011fe and, more generally, in normal SNe Ia. Could the identification by Hsiao et al. (2013) be wrong? Or is the fraction of neutral carbon underestimated in TARDIS, perhaps because it does not fully take into account all NLTE effects?

In support of their identification, Hsiao et al. (2013) use a physical argument first advanced by Marion et al. (2006), that the fractions of both neutral carbon and neutral oxygen can be quite high because high opacities by Fe line transitions in the ultraviolet deplete the photons that can ionize these elements. Here, Marion et al. (2006) noted that carbon and oxygen have similar ionization and excitation levels, and that the oscillator strength of the C I $\lambda 10693$ transition is about 50 times larger than that of the O I $\lambda 7773$ one. They thus suggest that the detection of the O I feature in SNe Ia implies that, if any carbon is present in the same layers as oxygen, it would be partially neutral and produce a strong feature.

Since Marion et al. (2006) found no clear evidence of carbon in the objects they studied, they concluded that those events had low carbon-to-oxygen fractions. But using the same argument that neutral carbon might be expected, Hsiao et al. (2013), Marion et al. (2015) and Hsiao et al. (2015) used line identification codes to infer that neutral carbon was present in the ejecta of SN 2011fe, SN 2014J and iPTF 13ebh, respectively, because it improved spectral fits in the NIR region.

As noted, our simulations do not show significant neutral carbon, but they do show an appreciable fraction of neutral oxygen. For instance, at $t_{\rm{exp}}=5.9$ d, in the outer layers with $v\gtrsim 19000{\rm\,km\,s^{-1}}$ which have the most favorable conditions, about 10% of oxygen but only 0.02% of carbon is predicted to be neutral. Such differences are expected in thermal equilibrium despite carbon and oxygen having similar ionization potentials ($11.2$ and $13.6{\rm\,eV}$ respectively). For instance, for a typical temperature in the outermost layers of $T=7000{\rm\,K}$, the Saha equation estimates a neutral oxygen fraction about 100 times larger than the neutral carbon fraction.

To determine whether the identifications are correct, it may help to compare with sub-luminous SNe Ia, for which, given their cooler ejecta, the fraction of neutral carbon should be higher than in normal SNe Ia. Indeed, the subluminous SN 1999by was one of the first SNe Ia for which neutral carbon was identified, by Höflich et al. (2002), using a self-consistent radiation code that included a detailed NLTE treatment.

Interestingly, Hsiao et al. (2015) detected neutral carbon in iPTF 13ebh, which is considered a transition object: it is on the fainter/cooler end of the Phillips (1993) relation, but its NIR maxima peaked before the $B$-band maximum and it did not exhibit a strong Ti trough. The C I $\lambda 10693$ feature is strong in the earliest spectra (i.e. prior to $-10$  d with respect to the $B$-band maximum), but quickly fades away, showing little velocity evolution. This behaviour is closer to what we would expect from cooler SN 2011fe models, where the temperature across the ejecta is artificially lowered so that the spectra resemble that of fainter SNe Ia (as described in Paper I).

Given the above, it is worth considering whether something else than C I could cause the observed flux deficit near 1.03$\rm{\mu m}$ found in SN 2011fe and SN 2014J. One alternative might be high velocity Mg. This was discussed by Marion et al. (2015) but discarded because of the lack of other high-velocity features at the same epoch (although the authors note that high-velocity Si might be present, but they cannot be sure because of possible blending with Na I).

Alternatively, the flux deficit in the NIR might be due to helium. Studies of the normal SN 1994D by Meikle et al. (1996) and Mazzali & Lucy (1998) already hinted that He I (and/or Mg II) could impact the NIR spectra under NLTE conditions, although it was noted that the absence of this feature in post-maximum spectra set quite stringent constraints. More recently, Boyle et al. (2017) investigated the possible signatures of He I in NIR spectra in the context of double detonation models. Their synthetic spectra exhibited a stronger (deeper) He feature than what is seen in SN 2011fe (compare Fig. 4 of Marion et al. 2015 and Fig. 11 of Boyle et al. 2017), but perhaps this would be different for other progenitor models. Unfortunately, we cannot test this directly in our models, since for helium, with its unusual excitation levels, a more sophisticated NLTE treatment is necessary (Boyle et al., 2017).

From the discussion above, we conclude that in cooler SNe Ia the C I $\lambda 10693$ signature has likely been detected; This is also corroborated by the identification of the other C I features in SN 1999by and iPTF 13ebh. However, for normal SNe Ia a more detailed investigation is needed, in which the relevant parameter space is explored with a code that can produce self-consistent spectra and adequately take into account NLTE effects.

7 Model uncertainties

Our simulated spectra have systematic uncertainties related to the approximations used in TARDIS. One is related to the line treatment. For the simulations shown here, we adopted the downbranch treatment, because this is the default approximation used in other Monte Carlo codes in the literature, thus allowing better comparisons. In principle, however, the “macroatom” treatment is more physical, as it allows for de-excitation cascades and multiple-photon excitation. As already shown in §4 and Fig. 4 of paper I, the choice of line treatment does not have a large impact on spectra formed at the temperatures found in the ejecta of SN 2011fe. We find that also for this study the results do not depend on the treatment adopted. To illustrate this, we show in Fig. 7 the spectra of our fiducial model for SN 2011fe with both line treatments. One sees that the main difference is in the flux level at the redder wavelengths999The different continuum shape means that with macroatom the fits to the spectra are no longer correct, i.e., for a proper analysis the tomography would have to be redone.; the shape of most spectral features are very similar.

Fig. 7 also illustrates the effect of changing the placement of the boundary between the inner and middle regions from the $v=13500{\rm\,km\,s^{-1}}$ used in our analysis. This placement was somewhat arbitrarily taken from our “depth scan” (see Sect. 3.1). When we move this boundary to $v=15000{\rm\,km\,s^{-1}}$ in our default model, where $X_{\rm C}^{\rm i}=0.002$ and $X_{\rm C}^{\rm m}=0.01$, the carbon feature weakens earlier and is shallower than observed at $t_{\rm{exp}}=$ 9 and 12.1 d. Conversely, when we place it at $v=12000{\rm\,km\,s^{-1}}$, the carbon feature becomes stronger than observed for $t_{\rm{exp}}\geq$ 9 d. Both changes are as expected, and show that in interpreting our results, one has to keep in mind that our abundances are averages over fairly large regions.

Another possible source of uncertainty in our synthetic spectra arises from the incomplete treatment of NLTE effects. It is known that deviations from LTE can be important: for instance, Thomas et al. (2007) studied the CO-rich zone of their W7 model and found significant departures when computing the population of atomic levels of C II that are relevant for producing lines in the optical. While Thomas et al. do not mention the explicit velocities and epoch used to compute their departure coefficients, the CO-rich zone of W7 is at velocities $v\gtrsim 14000{\rm\,km\,s^{-1}}$ in the W7 model (and $v\gtrsim 19000{\rm\,km\,s^{-1}}$ in our simulations), where the densities and temperatures are the lowest and thus more likely to deviate from LTE.

Despite the CO-rich zone being composed of nearly $100\%$ carbon in our calculations, it still does not contribute significantly to the carbon feature. This is consistent with the observations: for instance, using Fig. 4 of Parrent et al. (2012), one can see that Doppler velocities of the carbon trough do not exceed $15000{\rm\,km\,s^{-1}}$. Thus, NLTE effects are unlikely to greatly strengthen the carbon feature. We cannot exclude, however, the possibility that they greatly weaken it, in which case all our carbon mass fractions would be underestimates. Fortunately, for the deeper layers in which the feature is more likely to form, NLTE effects are likely to be less important.

Finally, as discussed in §6, one specific NLTE effect is that line blanketing by Fe depletes the pool of UV photons that could otherwise ionize C and O, possibly causing our calculations to underestimate the amount of neutral carbon everywhere in the ejecta. TARDIS only uses a simple approximation to treat this problem, which is implemented (following Mazzali & Lucy 1993) as a correction factor $\delta$ when estimating the balance of ions under the nebular ionization mode. We tested the influence of this approximation by comparing our fiducial spectra to models that enforce $\delta=1$ (i.e. no correction). The spectra did not differ in any significant way and thus this does help to determine neither whether depletion of UV photons is accounted for properly, nor whether the fraction of neutral carbon is estimated correctly in TARDIS.

Despite these uncertainties, we stress that compared to work that relies on line identification routines (e.g. Parrent et al. 2011, 2012; Hsiao et al. 2013; Marion et al. 2015), our approach has the advantage of computing temperatures and ionization balances in a more consistent way. Compared to full tomographic analysis (e.g. Stehle et al. 2005; Mazzali et al. 2008; Tanaka et al. 2011; Sasdelli et al. 2014; Ashall et al. 2016), our focus on one particular element helps to derive not just approximate values but also get a sense of the range of acceptable values. Nevertheless, in all these approaches, for a more complete assessment of their results, it will be important to compare with detailed NTLE treatments.

8 Conclusions and Ramifications

In this work we have used a previously published tomographic analysis (Mazzali et al., 2014) to investigate in detail the distribution of carbon at depth in the ejecta of SN 2011fe.

Based on the local temperatures and densities in the ejecta of SN 2011fe, we find that the carbon feature should be strong if significant amounts of carbon are present near the photosphere for epochs up to maximum. We thus interpret the lack of observed events with a clear feature near maximum as evidence for the lack of carbon deep in the ejecta.

We also find that for the relevant conditions, carbon is predominantly singly ionized and unlikely to produce the near-infared C I $\lambda 10693$ line; This agrees with the findings of Tanaka et al. (2008) and leaves as a bit of a mystery the flux deficit near $1.03{\rm\,\mu m}$ observed in some normal SN Ia. It remains an open question whether it is due to neutral carbon (if full NLTE effects are relevant), neutral He (Boyle et al., 2017) or high velocity Mg (Marion et al., 2015).

We compare our results with predictions for various SN Ia explosion scenarios and find that the violent merger model of Pakmor et al. (2012) predicts, on average, too much carbon mixed across the whole ejecta. In contrast, the spontaneous detonation of cold white dwarfs with no He shell, which might be seen as a first-order approximation to, e.g., the dynamically-driven double-detonation model of Shen et al. (2018), produce too little carbon at the velocities where we detect it. The discrepancies for these models are so large (factor of $\sim\!10$) that we believe these are excluded even taking into account the approximations used in our approach.

The delayed detonation model of Seitenzahl et al. (2013) and the gravitationally confined detonation model of Seitenzahl et al. (2016) are roughly consistent with our results, but might well be excluded by determining the distribution of unburned material in sub- and over-luminous events, or by investigating in detail the outermost layers of the ejecta. For instance, the double detonation models predict that the carbon abundance should peak at intermediate velocities, with the velocity of maximum abundance likely depending on mass of the exploding white dwarf (and thus peak luminosity; Fink et al. 2010), whereas in DDT models, carbon abundance is predicted to be higher for fainter SNe Ia (e.g. Wheeler et al. 1998; Höflich et al. 2002).

We thank Stephen Ro for insightful discussion regarding the physics of exploding WDs. This work made use of the Heidelberg Supernova Model Archive (HESMA), https://hesma.h-its.org. W. E. Kerzendorf was supported by an ESO Fellowship and the Excellence Cluster Universe, Technische Universität München, Boltzmannstrasse 2, D-85748 Garching, Germany. SAS acknowledges support from STFC via grant ST/P000312/1. .

Appendix A Outer zone

As argued in §5, carbon in the outer parts of the ejecta experience a lower density and temperature, meaning that under the NLTE approximations in TARDIS, the excitation levels that contribute to forming the optical feature are less likely to be populated. However, this effect may be compensated by greatly increasing the mass fraction of carbon in these layers. Here, we investigate the possible contribution from carbon at relatively high velocities by adding another zone to our analysis, where the mass fraction of carbon ($X_{\rm C}^{\rm o}$) is varied in a region defined by $16000\lesssim v\lesssim 19000{\rm\,km\,s^{-1}}$. For reference, M14 adopts $X_{\rm C}\gtrsim 0.02$ at comparable velocities.

As in §3.2, we perform our analysis pair-wise, i.e., we investigate the acceptable mass fractions in the outer region, while allowing the carbon mass fraction in the middle region to vary. These zones are best constrained by the early spectra, for which the photosphere is the closest to these zones and at a time of minimal dilation. Note that the inner region is basically below the photosphere at $t_{\rm exp}$=3.7 d and thus of less importance here.

In Fig. 8 we show the synthetic spectra for which $0.01\leq X_{\rm C}^{\rm o}\leq 0.5$ (color coded) and $0.0002\leq X_{\rm C}^{\rm m}\leq 0.05$. This range covers both the extremes of a carbon feature not forming at all in the early epochs, or the feature being either too strong or too blue. Here one has to be careful, because similarly to the spectra in M14, the red end of the main Si $\lambda$6355 is inconsistent with the data and will make a precise assessment of the blue end of the carbon feature less reliable. Correcting this discrepancy is beyond the scope of this work. Unfortunately, it also means that we can only place a conservative upper limit on $X_{\rm C}^{\rm o}$ and lower limit on $X_{\rm C}^{\rm m}$.

We conclude that $X_{\rm C}^{\rm o}\ \geq 0.4$ leads to a carbon trough that is inconsistent with the data at $t_{\rm exp}$=3.7, 5.9 and 12.1 d for any choice of $X_{\rm C}^{\rm m}$. Assessing a lower limit here is also difficult, given that none of the simulations seem to reproduce particularly well the spectra $t_{\rm exp}$=3.7 d.