The shape of SN 1993J re-analysed Free

ABSTRACT

SN 1993J is one of the best-studied Type IIb supernovae. Spectropolarimetric data analyses were published over two decades ago at a time when the field of supernova spectropolarimetry was in its infancy. Here, we present a new analysis of the spectropolarimetric data of SN 1993J and an improved estimate of its interstellar polarization (ISP) as well as a critical review of ISP removal techniques employed in the field. The polarization of SN 1993J is found to show significant alignment on the q − u plane, suggesting the presence of a dominant axis and therefore of continuum polarization. We also see strong line polarization features, including H β, He i λ5876, H α, He i λ6678, He i λ7065, and high velocity (HV) components of He i λ5876 and H α. SN 1993J is therefore the second example of a stripped-envelope supernova, alongside iPTF13bvn, with prominent HV helium polarization features, and the first to show a likely HV H α contribution. Overall, we determine that the observed features can be interpreted as the superposition of anisotropically distributed line forming regions over ellipsoidal ejecta. We cannot exclude the possibility of an off-axis energy source within the ejecta. These data demonstrate the rich structures that are inaccessible if solely considering the flux spectra but can be probed by spectropolarimetric observations. In future studies, the new ISP corrected data can be used in conjunction with 3D radiative transfer models to better map the geometry of the ejecta of SN 1993J.

1 INTRODUCTION

Core collapse supernovae (CCSNe) are the result of the death of a massive star (M_ZAMS > 8 M_⊙). They are classified as Type II or Type I according to whether their spectra show hydrogen or not (Filippenko 1997). Type I CCSNe are further subdivided into Type Ib (helium) and Type Ic (no helium – Filippenko 1997; Liu et al. 2016), and are the product of massive stars that have been stripped of their outer envelope, either through binary interaction or strong stellar winds (for a review see Smartt 2009). Type IIb SNe are a transitional class found to show spectra similar to Type II SNe within a few weeks after explosion, which then evolve to become dominated by helium as seen in Type Ib SNe.

SN 1993J is the second Type IIb SN to have been observed, after SN 1987K (Filippenko 1988). It is also one of the best-studied objects of this class owing to its proximity, as it was located near M81, ∼4 Mpc away (Bartel et al. 1993). Indeed, nearly 10 years of photometric measurements were conducted (Richmond et al. 1996; Zhang et al. 2004), as well as over 6 years of spectroscopic observations (Matheson et al. 2000). Hubble Space Telescope images also revealed the presence of a binary companion (an M-supergiant – Maund et al. 2004), and the disappearance of the progenitor of SN 1993J (a K-supergiant – Maund & Smartt 2009). In addition to classic spectrophotometric follow-up, SN 1993J was also observed with spectropolarimetry by Trammell, Hines & Wheeler (1993) and Tran et al. (1997) – subsequently T93 and T97, respectively. This makes SN 1993J the second supernova (and the first Type IIb) for which a clear, intrinsic, spectropolarimetric signal was detected.

Over the past 30 yr, spectropolarimetry has proven to be an invaluable tool to probe the three-dimensional geometry of unresolved supernova ejecta. The significant electron opacity present during the photospheric phase results in the emitted light being linearly polarized (Shapiro & Sutherland 1982). If the ejecta are spherical and homogeneous (see Fig. 1 – left-hand side illustration), the polarization components cancel out and no signal is detected. If significant asymmetry is present, however, incomplete cancellation of the normalized Stokes parameters (q and u – see Fig. 1, right-hand side) will produce an observable signal.

Asymmetry or asphericity in CCSN ejecta are usually inferred to be the result of either a global distortion of the ejecta, or of a non-isotropic distribution of the material within the ejecta (i.e. clumps). This has led to three main base cases, as illustrated in Fig. 1:

• Aspherical electron distributions, such as an ellipsoidal photosphere (e.g. Höflich 1991).

• Partial obscuration of the underlying Thomson-scattering photosphere (e.g. Kasen et al. 2003).

• Asymmetric energy input, such as heating by an off-centre radioactive decay source (e.g. Höflich 1995).

While Case (ii) produces line polarization, Cases (i) and (iii) alone result in continuum polarization (although their superposition can also create line polarization – see e.g. Stevance et al. 2019 for a discussion). All three base cases presented here result in significant (projected) axial symmetry,¹ which can be broken when multiple cases (or instances of the same case) are superposed.

Ultimately, the goal of spectropolarimetry is to infer geometrical information that is not accessible in the total flux in order to better understand and constrain the explosion process. It is the only observational method capable of directly probing the three-dimensional geometry of young supernova ejecta. On the whole, Type II SNe tend to show low polarization (∼0.2 per cent or less) at early days with a jump in the polarized signal after the plateau phase to levels around 1 per cent (e.g. Jeffery 1991; Leonard et al. 2006). Additionally, observations of stripped-envelope SNe have revealed that they exhibit high levels of continuum polarization at earlier phases compared to Type II SNe (e.g. Wang et al. 2001; Chornock et al. 2011; Mauerhan et al. 2015; Stevance et al. 2019). This has led to a widely stated conclusion that the degree of polarization increases as we probe regions closer to an intrinsically asymmetrical explosion, i.e. deeper into the ejecta (e.g. Wang & Wheeler 2008).

Type IIb SNe, although only making up ∼12 per cent of CCSNe (Eldridge et al. 2013), are remarkably well represented in the spectropolarimetry literature. The increasing number of high-quality data sets has revealed a significant diversity in their spectropolarimetric properties (for a discussion, see Stevance 2019; Stevance et al. 2019). The greater level of complexity in these data has also provided examples of SNe that do not follow the classic paradigm of ‘greater polarization at greater depths’, such as SN 2011dh and SN 2011hs (Mauerhan et al. 2015; Stevance et al. 2019). There remains a lot to learn, but these advances in our understanding of Type IIb SN spectropolarimetry offer new insights that can be used to provide a modern interpretation of the SN 1993J data.

We have so far focused on the polarization arising from the ejecta of SNe, but it is worth noting that bright light echoes – i.e. dust scattering in the circumstellar medium – can also produce observable polarization if the dust is asymmetrically distributed or inhomogeneously aligned (e.g. Wang & Wheeler 1996 and Yang et al. 2018 – see their fig. 4). In this context, the polarizing effects of dust can reveal crucial information about the environment of the progenitor star.

However, dust induced polarization usually occurs in the interstellar medium (ISM) of the host galaxy and of the Milky Way through dichroic absorption; it is then a contaminant that needs to be removed. This interstellar polarization (ISP) can be very difficult to estimate and assumptions must be made. In the case of T93 and T97, the two main assumptions were a Galactic-like ISP and an intrinsically unpolarized H α emission component. Both of these conjectures have now been found not to be necessarily justified (e.g. Patat et al. 2009, 2015; Tanaka et al. 2009; Stevance 2019; Stevance et al. 2019). Reservations about the T93 and T97 ISP estimates have already been expressed by Wang & Wheeler (2008); a few years later Chornock et al. (2011) remarked that the polarization signal of SN 1993J was unexpectedly high near 6600 Å compared to other Type IIb SNe – they suggested this could be resolved with a different ISP estimate.

Since the original spectropolarimetric studies of SN 1993J were published, ISP removal techniques have been significantly refined. In addition, spectropolarimetric data analysis methods and science focuses have evolved as more high-quality observations were obtained for CCSNe (e.g. Wang et al. 2003; Leonard et al. 2006; Maund & Smartt 2009; Tanaka et al. 2009; Chornock et al. 2011; Mauerhan et al. 2015; Reilly et al. 2016; Stevance et al. 2017). Informed by the lessons learnt in the field since the original SN 1993J spectropolarimetric studies were published, we aim to provide an updated analysis and interpretation of the data.

This paper is organized as follows: in Section 2 a summary of the observations is given. In Section 3, we review the ISP removal tools developed over the past few decades and their associated caveats, before addressing the original ISP determination techniques used in T93 and T97 and offering a new estimate. In Section 4, we analyse the resulting intrinsic degree of polarization of SN 1993J before interpreting it in Section 5 in the context of the base cases described above (see Fig. 1). We then summarize our work in Section 6.

2 OBSERVATIONS AND DATA REDUCTION

It is important to note that whenever statistics were evaluated (e.g. averages), the calculations were made using q and u before converting to a degree of polarization and debiasing. We remark that we received normalized Stokes parameter measurements, rather than the individual ordinary and extraordinary ray fluxes. The non-linear combination of fluxes performed to obtain the Stokes parameters (see Patat & Romaniello 2006) requires binning to be applied to the original flux data. As a result, we were unable to re-bin the data we obtained and are bound to the original format of T97.

Lastly, we report that we were unfortunately not able to locate the data from T93. We attempted to digitize the degree of polarization p and P.A. from the original publication, however the normalized Stokes parameters we retrieved were consistent with noise. This is expected if the digitization of p and P.A. is not sufficiently precise such that the correct P.A. and p pairs are not properly identified. This then yields normalized Stokes parameter measurements that do not reflect the original signal, but the uncertainty in the digitization process.

3 INTERSTELLAR POLARIZATION

Although in a few instances the contribution of the ISP can be negligible (e.g. Patat et al. 2011; Mauerhan et al. 2015; Stevance 2019), it is often an important contributor to the observed polarimetry signal. Because of the effects of vector addition, it can either decrease or increase p, and turn polarization peaks (troughs) into troughs (peaks). In some extreme cases, the ISP can even substantially dominate over the features intrinsic to the supernova (Stevance et al. 2019). Accurately quantifying the ISP is therefore crucial to retrieving the intrinsic polarization of the object of interest. Unfortunately, this is no easy feat, and it is often judicious to use multiple estimators and compare their results (e.g. Leonard et al. 2001, 2002a; Mauerhan et al. 2015; Inserra et al. 2016; Stevance et al. 2017, 2019).

We will first summarize the toolbox of ISP determination techniques that has been developed in the past three decades in Section 3.1. In light of these advances, we will discuss in Section 3.2 the original ISP estimates of SN 1993J and their issues, as well as compare them to new values we derive using independent methods.

3.1 The challenge of removing the ISP

3.1.1 Null intrinsic polarization of strong emission lines

The emission line flux in SN P Cygni profiles comes from resonant scattering and/or re-combinations, both of which are non-polarizing processes. Therefore, the emission line flux dilutes the flux polarized by electron scattering (assuming no electron scattering occurs in the line forming region). Furthermore, incoming polarized flux traversing the line forming region will undergo resonant scattering, resulting in depolarization, in addition to the dilution.

Based on these principles, two methods can be adopted:

• The polarization signal associated with the emission line region may be assumed to be intrinsically null and used as a direct proxy for the ISP. Typically, the strongest emission lines are considered for this purpose as the large line flux will dilute and depolarize the continuum to a greater extent. H α is most commonly employed (e.g. Leonard et al. 2001, 2006; Stevance et al. 2016), followed by the calcium triplet (e.g. Leonard et al. 2002a; Mauerhan et al. 2015; Reilly et al. 2016), although Inserra et al. (2016) also exploited the Si ii λ6355 line in the superluminous SN 2015bn.

• Alternatively, one may attempt to separate the continuum and line polarization, where in this case the latter is expected to be intrinsically unpolarized and therefore probe the ISP. Isolating the line polarization can be difficult. In the past, various methods have been used, some making assumptions in their determination of the polarization, others making assumptions when estimating the continuum flux across the line of interest. Such techniques were employed by T93 and T97 and are discussed in further detail in Section 3.2.1.

Both of these methods are susceptible to a number of caveats. First, the assumption of null intrinsic line polarization may not always be valid, either because of incomplete depolarization (e.g. Tanaka et al. 2009), or due to line blending (e.g. Stevance 2019). Additionally, Höflich et al. (1996) warn that electron scattering may also occur in the line-scattering region and therefore the continuum polarization may not be entirely determined by the deeper layers of the ejecta. Consequently, although it may be appropriate to assume null line polarization in the case of very strong emission lines where electron scattering is known to be negligible (e.g. Kawabata et al. 2002), caution is required when using these methods.

3.1.2 Null intrinsic polarization in line blanketing regions

This ISP estimate once again leverages the fact that resonant scattering is a de-polarizing process. In this case, however, we consider the blue regions of the spectrum, where a blend of a multitude of iron lines – the line blanketing regions – results in depolarization (Howell et al. 2001). If we assume that this phenomenon completely removes all intrinsic polarization, then any observed signal must be solely caused by ISP.

This method is frequently employed when a line blanketing region can be identified (Howell et al. 2001; Wang & Baade 2005; Chornock et al. 2006; Maund et al. 2009, 2013; Inserra et al. 2016; Stevance et al. 2017). Although it is expected from theory that this be valid between 4800 and 5600 Å (Howell et al. 2001), the ‘blue region of the spectrum’ used for this method is usually defined empirically and somewhat arbitrarily, in an attempt to avoid spectral regions polarized by e.g. strong absorption lines. Therefore, the wavelength range used varies from study to study – although it is typically <5000 Å – and should be considered individually when this technique is being applied to a new object.

An issue with this method is that, even though depolarization of the continuum through line blanketing is a real effect, this region of the spectrum does sometimes show significant line polarization, as can be clearly seen in SN 1993J at +30 d (Tran et al. 1997), or in SN 2008aq at +27 d (Stevance et al. 2016). Although it cannot always be used consistently to measure the ISP, it remains a useful method and a good sanity check.

3.1.3 Null intrinsic polarization at late times

Whereas the first two estimators relied on resonant scattering to depolarize or dilute the intrinsic SN signature, this technique takes advantage of the fact that intrinsic SN polarization will naturally go down to zero when the electron density has decreased sufficiently such that electron scattering no longer dominates, i.e. in the nebular phase.

Although some studies benefitted from spectropolarimetric data taken in the nebular phase (e.g. Leonard et al. 2006; Mauerhan et al. 2015), others have shown that late-time spectra obtained as the SN is transitioning from the photospheric to the nebular phase can also be used. For example, Chornock et al. (2011) and Stevance et al. (2017) focused specifically on emission line regions. Not only are they expected to be depolarized (at least to some extent – see Section 3.1.2), but they also become more prominent as the SN transitions to the nebular phase. By then, the photosphere has receded through the ejecta sufficiently that we can anticipate these emission line forming regions to be located well above the high optical depth regions where electron scattering is expected to be strongest. Consequently, the photons emitted in these lines are assumed to only carry the imprint of the ISP.

Lastly, Stevance et al. (2019) obtained a good ISP measurement by fitting over the full wavelength range of a late-time spectrum 40 d after V-band maximum, but noted that this method should be used with caution, particularly if line polarization features remain somewhat visible, indicating some level of intrinsic polarization.

3.1.4 Using the q − u plane

Another method, first described in Wang et al. (2001), takes advantage of the dominant axis (if observed). It is defined as an alignment of the data observed on the q − u plane which is expected to arise as the result of axial symmetry.

This technique assumes an intrinsic polarization with varying amplitudes but a constant direction across all wavelengths at each epoch (it can change between epochs). The ISP must then fall on one end of the dominant axis (Wang et al. 2001) – this is akin to shifting the origin of the coordinate system such that the data obey our expectation of constant P.A. The most appropriate ISP estimate can then be found by comparing several epochs (e.g. Leonard et al. 2002a), or making use of other theoretical assumptions such as the expected depolarization at the blue end of the spectrum from line blanketing (e.g. Howell et al. 2001; Wang et al. 2001). It is important to note that this is only valid if there are no large-scale inhomogeneities, which may not be the case in SNe that show prominent line polarization features.

3.1.5 Assuming Galactic behaviour

The ISP is made up of two components: the host galaxy and the Galactic ISP. The latter can be estimated by searching for local stars near the line of sight of the SN that have polarization measurements. We can assume typical stars are not intrinsically polarized and therefore provide an estimate of the Galactic ISP (Heiles 2000). Unfortunately, it can be difficult to find a polarization standard within less than 2° of the SN, and if the star is too close to Earth it will not sample the full dust column. Therefore, the two ISP components are usually not disentangled, and it is the total ISP that is considered.

Although these relationships may be valid for the Galactic component of the ISP, they may not be for the ISP originating from the host – meaning they may not be reliably used with the total ISP. There is indeed increasing evidence that the Galactic limit (equation 6) is not universally applicable to other galaxies (e.g. Leonard et al. 2002b; Stevance et al. 2019). Additionally, placing this limit necessitates an estimate of E(B − V) that often requires the use of empirical relationships between the interstellar Na i D lines and the colour excess (e.g. Barbon et al. 1990; Poznanski, Prochaska & Bloom 2012). It is also worth noting that Na i D lines are caused by discrete gas clouds and do not probe the more diffuse ISM, further adding to the uncertainty of such estimates.

It was already raised by Leonard et al. (2005) that equation (6) may not actually provide meaningful constraints to the ISP. In light of the 15 years of SN spectropolarimetric research that followed, we conclude that this practice offers little insight on extragalactic ISP. It can however serve to reveal the non-Galactic properties of interstellar dust in the host galaxies of SNe.

3.1.6 Continuum removal by means of wavelet decomposition of flux spectra

For completeness of this review, we must reference the generalized approach for continuum removal presented in Cikota et al. (2019). They derive the continuum component of the ordinary and extraordinary flux spectra using wavelet decomposition, from which they calculate the continuum component of the Stokes q and u. Finally, they subtract the continuum Stokes parameters from the total q and u to derive the line polarization spectra.

It is important to note that this technique removes both the intrinsic continuum polarization and the ISP, but it is particularly useful when studying objects where only the line polarization is of interest. In the case of SN 1993J, we did not have the opportunity to try this method as it requires the raw flux spectra, and they are no longer available.

3.2 The ISP of SN 1993J

3.2.1 Previous estimates and removal

Spectropolarimetric data of SN 1993J were published in two main studies by T93 and T97. As noted in Section 3.1.1, both teams attempted a decomposition of the continuum polarization from the H α polarization, as they assumed the latter would be a good proxy for the ISP. The isolation of the H α polarization was however done differently in both papers.

T93 assumed that the average polarization from 4900 to 6800 Å (p = 0.9 ± 0.1 per cent) in their April 20 data (+24 d) was equivalent to the continuum polarization. They removed this suggested continuum and measured the H α polarization (p_max = 1.1 per cent, θ = 150°) as proxy for the ISP. There are two main issues with this method: First, assuming that the average polarization between 4900 and 6800 Å is a proxy for the continuum may not be valid if strong line polarization is present; as seen in fig. 1(a) of T93, significant line polarization does seem to be visible. Secondly, as discussed in Section 3.1.1, incomplete depolarization of the SN continuum by the line forming region and/or blending with other spectral features can be an issue when employing this method. In the present case, we note that on April 20 (+24 d) – the date used by T93 – the H α emission profile is flat topped, and only 4 d later a significant He i λ6678 notch is visibly cutting into the hydrogen emission. The flat top at +24 d is indicative that the helium feature is already present (Swartz et al. 1993), which may invalidate the assumption of complete depolarization.

T97, on the other hand, attempted to separate the H α line polarization from the continuum in a different manner. Their method, detailed in their appendix, relies on estimating the continuum flux across the emission line profile. The uncertainties associated with this step – i.e. ‘drawing the correct continuum’ (T97) – can be problematic, as remarked by T97 and Stevance et al. (2019). This alone can lead to unreliable estimates of the ISP. Additionally, T97 specifically caution that this approach may not handle the superposition of emission and absorption lines, and as mentioned above the He i λ6678 is starting to appear in the spectral profile of H α on April 20 (+24 d). The ISP calculated by T97 at this epoch is p = 0.63 per cent and θ = 171°, which is significantly different from that of T93 (p_max = 1.1 per cent, θ = 150°).

Although their ISP estimators are different, both T93 and T97 used the same ISP removal method: A Serkowski law – see Section 3.1.5 – with p_max and θ set to their H α polarization, and λ_max = 5500 Å. The latter value was assumed based on the Galactic median value of λ_max = 5450 Å reported by Serkowski et al. (1975). However, as mentioned in Section 3.1.5, we now have numerous cases of SN host galaxies exhibiting an ISP that does not follow Galactic expectations.

3.2.2 ISP from line blanketing regions

We can employ the method described in Section 3.1.2, which assumes that the blue part of the spectrum is completely depolarized due to line blanketing (Howell et al. 2001). Wang & Wheeler (2008) performed a similar calculation to compare the T93 and T97 values, but their estimates were not explicitly reported (see their section 4.4.1).

We applied this method to the April 20, 26, 30, and May 11 data. A weighted average of the normalized Stokes parameters was performed; the errors on the weighted mean were also calculated, and propagated through to the degree of polarization and polarization angle using equations (1) and (2). The degree of polarization was also debiased following equation (3). The wavelength ranges used and the resulting ISP values are given in Table 2. Note that the average on April 20 (+24 d) is over a different range from the other epochs because no data are available between 5300 and 5500 Å.

The results obtained with this method are discussed and compared to the other ISP estimates in Section 3.2.4.

3.2.3 ISP from late-time data fitting

Another way to derive the ISP is to use the late-time data method described in Section 3.1.3. Adopting this approach, we make the assumption that at our latest epoch (1993 May 14/48 d after explosion), the polarization observed is dominated by the ISP rather than intrinsic SN polarization. This is true if the electron density has decreased sufficiently such that very little light is being polarized in the ejecta. Similarly to Stevance et al. (2019), we want to fit the late-time data with a straight line (normalized Stokes parameters versus wavelength) to estimate a wavelength-dependent ISP and we use a sigma clipping method since the data on 1993 May 14 have a low signal-to-noise ratio (SNR). Additionally, the data below 4500 Å at this epoch show a notable increase in noise (see Fig. 2). Since the spectral features of interest are located above 4500 Å, and given that two of our data sets have no data below 4600 Å (see Table 1), we consider appropriate to only fit the data with wavelength >4500 Å.

The ISP normalized Stokes parameters derived here are linear functions that do not make use of the Serkowski law (see equation 5 in Section 3.1.5). It is worth noting that the difference in gradients for q_ISP and u_ISP also results in a wavelength dependence in the P.A. of the ISP, as can be seen in the bottom panel of Fig. 3. This can be caused by the superposition in the line of sight of dust clouds with different particle sizes and orientations, or by the superposition of ISP and some latent intrinsic SN polarization (Coyne & Gehrels 1966). A similar effect was seen in Stevance et al. (2018) for the spectropolarimetric data of WR102 which contained no observable intrinsic polarization. Additionally, as we will discuss in the next section, the ISP estimates given here are supported by independent measures. On the whole, this change in P.A. over wavelength is likely to originate from the complex nature of the ISM in our line of sight.

3.2.4 Comparison of the SN 1993J ISP estimates

We now compare the various ISP estimates of SN 1993J performed in the present and previous works (see Table 2), in order to determine which ISP is most suitable.

The ISP values for the four dates obtained in Section 3.2.2 from the line blanketing method (see Section 3.1.2) significantly differ in amplitude, and an anticlockwise rotation with time of the polarization angle is observed. This is inconsistent with the properties of the ISM, which are constant on short time-scales, and shows that the results were biased by intrinsic polarization. It is interesting to note that the values obtained from May 11 (Section 3.2.3) are very close to the ISP values we derived at 5250 Å. This can be understood as the result of a decreasing amount of intrinsic polarization at later dates due to the decreasing electron density, causing the observed polarization to be dominated by ISP, independently of line blanketing.

In order to compare the three other ISP estimates (i.e. obtained from the late-time data in Section 3.2.3, T93, and T97) we plot in Fig. 4 their resulting ISP removed polarization on April 30 (+34 d). We choose April 30 for two reasons: First, it is important that the data used to evaluate the effectiveness of our estimates to remove the ISP be independent of the data used to calculate the estimates used by T93 and T97; secondly, Chornock et al. (2011) remarked that the level of polarization on April 30 near 6600 Å in the H α peak (∼0.7 per cent) reported by T97 was high (as compared to other Type IIb SNe). They suggested that a different ISP estimate for SN 1993J could lower the polarization levels of SN 1993J in this spectral region.

Consequently, we present the April 30 ISP removed data, corrected using the T93, T97, and our late-time estimate ISPs (see Fig. 4). The T97 ISP removed data were obtained directly from H. Tran (private correspondence). In order to remove the T93 ISP, we employed the same method described in their paper: A Serkowski law (see equation 5 in Section 3.1.5) with the wavelength-dependent K value from Wilking, Lebofsky & Rieke (1982) was used, and vector subtraction was performed according to equation (4) in Stevance et al. (2019). The degree of polarization was then re-calculated as prescribed in equation (1).

To compare the validity of these ISP estimates, we focus our attention on the regions devoid of strong lines in the blue, and the spectral range associated with the H α emission line. We expect these regions to be mostly depolarized (see Sections 3.1.2 and 3.1.1), although polarization features associated with e.g. He i λ6678 may be present.

First of all, it is very interesting to note that the T93 estimate actually significantly increases the overall degree of polarization compared to the observed data. The T97 ISP on the other hand does very little to reduce the overall level of polarization. By contrast, our ISP correction lowers the overall p levels of SN 1993J, as predicted by Chornock et al. (2011) – see Fig. 4. This is particularly true in the regions of line blanketing devoid of strong lines, and in the spectral region associated with the H α emission, where quasi-null levels of polarization are seen on either side of a well-defined peak correlated with the absorption component of He i λ6678. This is consistent with the expectations of intrinsic polarization, and further illustrates how the assumption of complete depolarization across the whole H α emission line can yield poor estimates of the ISP.

As another form of comparison between the ISP estimates discussed here, we can plot them on a q − u diagram (see Fig. 5) and consider them in the context of the method described in Section 3.1.4. In this case, the ISP is estimated by comparing several epochs of data on the q − u plane, knowing that the ISP must be constant over time and assuming that the ISP is located at one end of the dominant axis.

We can see in Fig. 5 that the T93 ISP falls away from the main locus of the data at all epochs. To a lesser extent this is also true of the T97 ISP, barring epoch 1 (+24 d). The line blanketing and late-time estimates are quite concentrated (although the line blanketing value found at +24 d is closer to the T97 value), and are located at one end of the locus of the data at +30 and 34 d. These estimates are, however, centrally located at +24, +45, and +48 d. At the latter two epochs this is expected if the data are, as suspected, dominated by ISP: in such a case the spectropolarimetric data should be randomly scattered around the ISP. Overall, Fig. 5 does not support the T93 and T97 ISP values, and favours our new estimates.

There is one last test of the ISP that we can perform. From theory, the absorption minima of strong lines approach the Rayleigh–Jeans limit and the scattered photons in the emission component are depolarizing. In spectral regions of strong line blending, multiple line scatterings of the emission component will cause strong depolarization. The same is true for strong lines with finite line width such as H α, as shown for SN1993J (Höflich 1995). From this model, we would expect the spectral region redward of the H α absorption component to show the lowest intrinsic polarization. However, in models the intrinsic line width depends both on the natural-line, pressure, as well as thermal and turbulent Doppler broadening, with the latter being assumed. For our last test of the ISP, we use a semi-empirical approach to minimize the variance in polarization across epochs redward of the H α absorption component (6395–6490 Å) between April 20 and all subsequent phases for a grid of q_isp, u_isp. The result is shown in Fig. 6. The empirical condition mentioned above is not sufficient to infer a specific ISP as the solutions are degenerate, but it is a useful sanity check to compare against our estimates, particularly since it is independent of theoretical assumptions.

In Fig. 6, we also plotted the ISP from the line blanketing method as well as the ISP calculated from equations (7) and (8) at 6450 Å. As we can see, the T93 and T97 values once again do not perform as well as our independent approaches. By contrast, our ISP based on the late-time fit falls right on the axis of minimum variance (light area in Fig. 6), whereas the values obtained from line blanketing are scattered around this axis.

On the whole, we rule out both the T93 and T97 ISP estimates. Additionally, ISP values derived from line blanketing are inconsistent with each other (see Table 2). The May 11 line blanketing estimate, however, is consistent with the late-time data derived ISP at similar wavelengths (5250 Å). Finally, the ISP presented in equations (7) and (8) resolved outstanding questions raised by Chornock et al. (2011) about the reported intrinsic polarization on April 30.

Consequently, we consider the ISP found from fitting the late-time data (equations 7 and 8) to be our best estimate, and it was used to correct all our epochs. Nevertheless, we note that this method (described in Section 3.1.3 and applied in Section 3.2.3) should be used with care and with full consideration of the caveats previously mentioned, particularly when fitting data that are suspected to show residual intrinsic SN polarization.

4 INTRINSIC POLARIZATION

4.1 Degree of polarization

In Fig. 7, the ISP corrected degree of polarization of SN 1993J is presented alongside the flux spectrum from +24 to +48 d (with respect to explosion date). At our first epoch the data are completely dominated by noise, but the data at +30 and +34 d show very clear line polarization features. Additionally, over these epochs, the helium features exhibit a clear increase in line strength corresponding to the transition from a Type II to a Type IIb spectrum.

At +30 d, helium and hydrogen polarization features are strong and complex, showing multiple peaks. The wavelength, velocity, degree of polarization (p), and polarization angle (P.A.) of the dominant features are summarized in Table 3. The strongest peak is that of H β; it is broad, spanning ∼100 Å and reaches a polarization level of 1.8 per cent. H α also exhibits a very broad feature, extending from ∼5960 to 6415 Å. The absorption minimum of the P Cygni profile of H α is correlated with a strong peak rising to a maximum polarization of p = 1.20 ± 0.13 per cent at 6304 Å. Three other distinct peaks are superposed on to the broad hydrogen feature at 6040, 6072, and 6152 Å, all with p ∼0.9 per cent (see Fig. 7 and Table 3). They are coincident with dips in the blue shoulder of the H α absorption line profile. We will further discuss the nature of these small features in Section 5.2.1.

Furthermore, He i λ5876 also shows a complex multipeaked feature at +30 d. The two main components are found at 5632 and 5728 Å reaching a similar degree of polarization, p ∼1 per cent (see Table 3 for more detail). The 5632 Å peak coincides with the bottom of the absorption component of He i λ5876, whereas the other peak is associated with a blue shoulder in the profile of the absorption component of the P Cygni profile. As a result, we refer to this second peak as a high velocity (HV) component of the He i λ5876 polarization. Weaker features (p ∼ 0.5 per cent) are found both on the blue and red side of those peaks, suggesting deviations from sphericity on a wide range of velocities (i.e. depth – assuming homologous expansion) in the ejecta.

Additionally, the spectral region associated with the H α emission, although it is mostly depolarized, shows a peak reaching p = 0.63 ± 0.07 per cent at 6600 Å that is likely associated with He i λ6678, which is seen to start developing in the flux spectrum at this epoch. Moreover, a polarization peak is associated with the onset of the He i λ7065 in the flux spectrum (see Table 3 and Fig. 7). Here, we have only described the main features of the polarization at +30 d, but these data are particularly rich of secondary polarization features that are difficult to interpret without a model. In the future, the spectropolarimetric data of SN 1993J, particularly at +30 d, should be revisited when adequate modelling capabilities have been developed.

The data at +34 d are noisier than those at +30 d. This is due to a combination of effects: the supernova was fading, the telescope was smaller, and the bins are smaller. None the less, very clear features similar to those observed at +30 d can be seen in the spectropolarimetric data at this epoch. For all features, excluding H α, the degree of polarization and P.A. are calculated from the weighted average of the normalized Stokes parameters in the corresponding range indicated in Table 3. The errors are propagated from the errors on the weighted mean of the normalized Stokes parameters.

At +34 d, hydrogen still exhibits a strong H β peak and a broad H α feature topped by a peak associated with the main absorption line. The H α feature actually seems to extend further in the blue than it did at +30 d: down to ∼5910 Å.

The He i λ5876 still shows two main peaks with degree of polarization similar to that observed at +30 d: one associated with the absorption minimum in the flux spectrum, and the other associated with a blue shoulder in the profile of the absorption line, interpreted as an HV component (see Table 3). Once again a He i λ6678 (∼6600 Å) peak is superposed on to the depolarized region associated with the H α emission component. Additionally, the broad polarization feature extending from 6715 to 7000 Å could be a signature of He i λ7065.

At +45 and +48 d, the main He i λ5876 and Hα peak remain visible but the noise and error bars are such that we will refrain from making quantitative estimates of the level of polarization for these lines. The decrease in polarized signal at these epochs is most likely due to the decrease in electron density in the ejecta.

4.2 Normalized Stokes q − u plane

The behaviour of spectropolarimetric data on the q − u plots can reveal large-scale and small-scale asymmetries (see Wang & Wheeler 2008 for a review). When plotting the whole data, an alignment on the q − u plane ordered with wavelength (dominant axis) can be interpreted as the result of (projected) axial symmetry; additionally, we look for ‘loops’, which we define as a gradual rotation in polarization angle across the wavelength range associated with a specific spectral line. They are a sign of departure from axial symmetry (Wang & Wheeler 2008; Tanaka et al. 2017).

4.2.1 The dominant axis

In Fig. 8, we present the ISP corrected data of SN 1993J from +24 to +48 d. Some alignment seems to be present in the data on the q − u plane at our earlier epochs, although a high level of noise and scatter is present. The data at +45 and +48 d on the other hand are centred on the origin of the plot and appear to be randomly distributed.

In order to evaluate the dominant axis in the q − u plots, we retrieved the principal components (PCs) of the data using singular value decomposition (SVD)³ on the normalized Stokes parameters and the wavelength. PCs are defined such that the first PC indicates the direction with the greatest variance. The first PC therefore corresponds to the dominant axis.

In addition to performing SVD, we also used bootstrapping in order to estimate the uncertainties on the estimated dominant axes. We processed each epoch 1000 times using a random sample containing 80 per cent of the data. The resulting PCs are overplotted on to the data in Fig. 9 and we summarize the calculated dominant axis orientations in Table 4.

The dominant axis at +30 d was calculated using data with wavelengths >4800 Å in order to exclude the H β loop since it shows a strong deviation from the rest of the data that would skew the PCs. Indeed, the dominant axis reflects the behaviour of the continuum but in Stevance et al. (2017) a single prominent line feature can bias the dominant axis fit. We further discuss the H β line in the next section.

The centroid of the data at +24 d is off-centre with respect to the origin and the data show a clear dominant axis. This suggests continuum polarization but the high level of noise prevents us from making a definitive interpretation. In subsequent epochs (+30 and +34 d), a small clockwise rotation of the dominant axis is observed. Then at +45 d and particularly at +48 d the uncertainty on the orientation of the dominant axis increases significantly, which is expected for data with enhanced levels of noise.

4.2.2 The loops

In Fig. 10, we present the H α, H β, and He i λ5876 data on the q − u plots. We only show the +30 and +34 d data since line features at subsequent epochs are noise dominated (see Fig. 7).

At both +30 and +34 d, the data corresponding to the main H α peak show a linear configuration (see Fig. 10) on the q − u plane which is in line with the +u direction (P.A. ∼45°) rather than the dominant axis at these epochs (∼66° and 61°, respectively – see Table 4). This indicates that the ejecta geometry probed by the H α line forming region is different from the configuration resulting in the observed dominant axis.

In Section 4.1, we saw that the strongest polarization feature at +30 d was that of H β. The H β data at this epoch (see Fig. 10) show one of the most complex loops observed in a CCSN, exhibiting a mixture of successive anticlockwise and clockwise rotations, with an overall orientation quasi-orthogonal to the dominant axis. At +34 d, the greater levels of noise make the case of a rotation in P.A. less obvious. However, it is clear that the data do not align with the dominant axis that is representative of the rest of the data. Therefore, at both +30 and +34 d, H β exhibits significant departures from axial geometry.

The He i λ5876 line at +30 and +34 d was associated with two main polarization peaks, the first corresponding to the absorption flux minimum and the second interpreted as an HV component, associated with a blue shoulder in the helium line profile (see Section 4.1). Both at +30 and +34 d, the main helium features show clear anticlockwise loops on the q − u plane (see Fig. 10). The loop at +34 d is more developed than 4 d prior, showing 180° rotation across the feature and a greater amplitude. Consequently, the main He i λ5876 peak also probes departures from axial geometry. The HV components at both epochs, on the other hand, show very little P.A. rotation, if any. Departures from axial symmetry are therefore less prominent, if at all present, than in the deeper layers of this line forming region.

5 DISCUSSION

5.1 Global geometry

5.1.1 Probing the continuum polarization

Previous studies of the spectropolarimetric properties of SN 1993J were published by T93 and T97. These analyses were performed at a time where intrinsic polarization had only been observed in one supernova – SN 1987A (Cropper et al. 1988). They therefore mostly focus on whether intrinsic polarization was present and on quantifying the degree of continuum polarization. In order to do this, T93 and T97 use the same method of averaging the data across a large range of wavelengths. In the former study, they use the range 4900–6800 Å and find p  = 1.6 per cent (ISP removed) on April 20 (+24 d), whereas in the latter, data are averaged over 4800–6800 Å and p ∼ 1 per cent is found on April 26 (+30 d).

Since these wavelength ranges encompass both the He i λ5876 and H α lines, this method will capture the presence of intrinsic polarization, however it is not suitable to isolate the continuum from the line polarization. In order to estimate the continuum independently, it is necessary to identify a region of the spectrum that is devoid of strong lines. A good spectral range for this purpose would be between the 7630 Å tellurics lines and the onset of the calcium infrared triplet absorption component. However, this region of the spectrum is not accessible in these data sets (see Fig. 7), and we therefore did not attempt to calculate a continuum polarization.

Determining a degree of continuum polarization is not the only way to probe the global geometry of the ejecta; it can also be inferred from the q − u plots. Specifically, the presence of a dominant axis is seen as a marker of axial geometry. In Section 4.2, we used SVD to identify the PCs of the data on the q − u plots. Although dominant axes are usually fitted in q − u space, this analysis was performed in q − u − λ space, because the wavelength dimension is informative and was necessary to obtain good fits to our data.

This wavelength dependence may be counterintuitive given that Thomson scattering is a grey-type of scattering, however, as explained by Patat et al. (2012), wavelength dependence in the continuum polarization arises from the formation of a pseudo-continuum by the blending of a large number of lines: in the blue this pseudo-continuum dominates the opacity, and it decreases to longer wavelengths. Additionally, it has been shown that polarization is dependent on optical depth (Höflich 1991): at too high optical depth (τ > 3), thermalization results in a loss of the polarization information; at too low optical depth, scattering is insufficient for significant polarization to be produced; around τ = 1.2 there is a sweet spot where polarization is maximized. From this, Patat et al. (2012) show that continuum polarization should reach a maximum around 7000 Å. This would explain why taking the wavelength dependence into account informs the fit of the dominant axis.

In this context, however, we would expect the origin of the q − u plots to coincide with the blue end of the data: we can see in Fig. 5 that this is not necessarily the case, particularly at 24 d. The cause for this discrepancy is not clear.

5.1.2 Interpreting the continuum polarization

Despite the caveat mentioned above, the clear elongation with wavelength of the data on the q − u plane at 24 d could indicate the presence of continuum polarization in the data which would not be visible in p (see Fig. 7) due to the level of noise in the data. Binning could have revealed a measurable degree of polarization, but as mentioned in Section 2, we were not able to properly re-bin the data we received.

The presence of a clear dominant axis at +24, +30, and +34 d is the marker of axial geometry. A small clockwise rotation is observed across these three epochs (see Fig. 10) as helium features start to dominate the spectrum. This is in interesting contrast with SN 2001ig that showed a drastic rotation of the continuum polarization angle by  40° between +13 and +31 d. This was interpreted as a decoupled geometry between the outer, hydrogen-dominated layers, and the inner, helium-dominated, layers of the ejecta (Maund et al. 2007b). At the time a parallel was drawn between the polarization properties of SN 2001ig and SN 1993J, however our re-analysis does not fully support that: SN 1993J shows a smooth evolution of its continuum polarization as the photosphere recedes down through the helium layers.

The axial geometry configuration probed at our first three epochs could either be the result of ellipsoidal ejecta – Case i – or of an off-centre energy source – Case iii (see Fig. 1). The spectroscopic study and modelling of SN 1993J by Swartz et al. (1993) support little mixing of heavy elements (such as radioactive nickel) to the hydrogen-rich material in the outer envelope, which would favour the ellipsoidal ejecta interpretation. However, Spyromilio (1994) found evidence of clumping of radioactive material in the emission line profiles of SN 1993J, which would support the off-centre energy source scenario. Without further modelling, spectropolarimetry cannot differentiate between these two cases.

5.2 Line specific geometry

Line polarization is discussed in T97. They remark that the helium λ5875 feature is more polarized than the Balmer lines, which they interpret as either being the result of more aspherical helium layers, or clumpy ejecta.

In this work, we saw in Section 4 (see Fig. 7 and Table 3) that the main He i feature actually shows a lower polarization peak than the main H α feature, although within uncertainties these values may be equivalent. The H β line polarization, however, does surpass the helium polarization, as well as that of H α by over 0.5 per cent. Consequently, our line polarization does not necessarily support a picture where the inner helium layers are more aspherical as suggested by T97. In the greater context of Type IIb SNe, the line polarization levels of He i and H α (p ∼ 1 per cent) recorded for SN 1993J are more prominent than in SN 2011hs (between about 0.5 to 1 per cent – Stevance et al. 2019), similar to SN 2001ig (∼1 per cent – Maund et al. 2007b) but less extreme than in SN 2008aq (>2 per cent – Stevance et al. 2016; Stevance 2019) or SN 2008ax (∼3 per cent – Chornock et al. 2011). Thus, in this regard SN 1993J appears to be a fairly average Type IIb SN.

More insight into the ejecta geometry can be obtained by looking at the behaviour of the line polarization on the q − u plane. We first consider the hydrogen features. Although the main peak of H α at +30 d behaves as expected in the case of axial geometry, H β shows a very complex loop – and so a complex departure from axial symmetry – at similar velocities, see Fig. 10. Therefore, the main feature of H α probes a different geometry from that of H β This means that global geometry effects are not the only source of polarization, and some line specific geometries must come into play. It is worth noting that the H β loop is not a frequently observed feature of Type IIb SN spectropolarimetry, and SN 2008ax is the only one that has shown a more prominent H β feature (Chornock et al. 2011).

In the case of He i λ5876, loops are seen at a depth below −9000 |$\mathrm{ km \,\, s^{-1}}\,$| at +30 and +34 d, which contrasts with the linear configuration followed by the H α data at similar depths. In addition, the helium loop strengthens between +30 and +34 d, as the helium lines in the spectrum become more prominent. Consequently, these polarization features seem to probe line specific geometries and departures from the global axial geometry. A similar loop evolution following the development of helium features in the flux spectrum has also been observed in other Type IIb SNe: SN 2008aq, SN 20011hs, and SN 2011dh (Mauerhan et al. 2015; Stevance 2019; Stevance et al. 2016, 2019); this phenomenon now stands out as a frequent feature of Type IIb SN spectropolarimetry.

Inhomogeneities in the distribution of hydrogen and helium resulting in uneven obscuration of the underling photosphere, would explain the observed departures from axial symmetry – Case ii (see Fig. 1).

5.2.1 Nature of the features in the range 6000–6300 Å

Understanding the origin of the small amplitude features seen to the blue of the absorption profile of H α (between ∼6050 and 6300 Å – see Section 4.1) is crucial to interpreting the polarization features associated with them. In T97, the most prominent peak found at 6220 Å was identified as a [O i] λ6300 feature. However, at +30 d the spectrum is still well within the photospheric phase; the collection of spectra of SN 1993J presented by Matheson et al. (2000) shows that [O i] only starts developing after 80 d.

To further investigate the nature of this feature, we performed fits using syn++ (Thomas, Nugent & Meza 2011). This code implements the elementary supernova model of Jeffery & Branch (1990). A brief description of this model and its physical assumptions is provided in Vallely et al. (2016), and more information can be found in Branch et al. (2004, 2005, 2007). We first fit the whole spectrum irrespective of the small amplitude features in the blue shoulder of H α; this is our base fit (see top panel of Fig. 11). This fit contains H i, He i, O i, Ca ii, and Fe ii. It is calculated with a photospheric velocity v_phot = 7000 |$\mathrm{ km \,\, s^{-1}}\,$| and a temperature T = 5000 K (blackbody). All values were chosen through visual comparison of the synthetic spectra to the data.

We then attempt to improve the fit in the spectral range 6000–6250 Å by individually adding HV H i, Fe ii, as well as Si ii and C ii, as done by Kumar et al. (2013) when investigating a similar profile in SN 2011fu. The results of these fits are shown in the middle and bottom panels of Fig. 11.

Both Fe ii and C ii yield unsatisfactory results, but the addition of either HV H α features or Si ii allow a much better fit to the small features. The HV hydrogen was introduced as an additional ion with v_min = −19 500 |$\mathrm{ km \,\, s^{-1}}\,$| whereas silicon required v_min = −10 000 |$\mathrm{ km \,\, s^{-1}}\,$|⁠, meaning that the silicon line forming region would be significantly detached from the photosphere. This is unexpected since Si resides in the core of the explosion, but could be the case if a plume had transported material from the core to the outer parts of the ejecta. However, it is worth noting that the P.A.s of the features at 6040, 6072, and 6152, respectively 32°, 37°, and 43°, show a smooth anticlockwise rotation that is also followed by the main H α feature, with a P.A. of 46° (see Table 3).

Consequently, we conclude that it is more likely that these features be the result of HV H α line forming regions. None the less, the velocity at peak for the polarization features is given both for the H α and Si ii case – see Table 3.

5.3 A picture of SN 1993J

Based on our interpretation of the global and line specific polarization in Sections 5.1 and 5.2, we can produce a rudimentary sketch of the ejecta of SN 1993J, which we present in Fig. 12. We created this picture by combining what we learnt from our first three epochs (+24, +30, and +34 d) as the noise in the last two epochs prevents us from coming to definitive conclusions.

First of all, we deduced in Section 5.1 that axial symmetry was present at +24, +30, and +34 d, with a small, smooth, rotation over time. This is shown as the rotation of the dominant axis in Fig. 12. The ejecta outlines presented in our sketch are aligned with the dominant axis at +24 d – we do not indicate the outlines for the following two epochs because the minimal rotation results in a significant overlap between the dashed lines that hinders visualization. The photosphere shown is located at 7000 |$\mathrm{ km \,\, s^{-1}}\,$|⁠, according to the photospheric velocity used to fit the flux at +30 d – see Section 5.2.1.

We chose to represent the ejecta through ellipses, in relation to a Case i geometry – see Fig. 1 – but note that we cannot exclude the possibility of an off-centre energy source (Case iii), see Section 5.1.

The dashed outlines representing the ejecta are ellipses with an axial ratio of ∼0.88. This number was chosen based on fig. 4 of Höflich (1991) and a visual examination of the q − u plots in Fig. 8, from which we can determine that the average continuum polarization at the first three epochs could be anywhere from ∼0.5 to ∼1 per cent.

The coloured regions in Fig. 12 represent the location of clumps resulting in the partial obscuration of the photosphere (Case ii) that is thought to be the origin of the line polarization recorded in Table 3. Note that artificial symmetry is induced by the nature of the Stokes parameters, which are pseudo-vectors (the 0° and 180° directions are the same).

It is interesting to see that the inhomogeneities in the line forming regions of H β, H α, and He i are orthogonal to the direction of elongation of the photosphere. The reason for this is not clear from our qualitative analysis. In order to better understand the geometry of the ejecta of SN 1993J, careful modelling will be required.

6 CONCLUSIONS

We presented a new analysis of five epochs of spectropolarimetric data of SN 1993J originally obtained by T97. The observations range from +26 to +48 d after explosion. We first reviewed the ISP estimation techniques used in the spectropolarimetry literature as well as the calculations made by T97 and T93. We concluded that the Serkowski assumption, as well as that of complete depolarization in the emission line profile of H α, are not suitable. We then evaluated the ISP using the line blanketing method and a late-time data fitting technique. Subsequently, all ISP estimates, including those of T93 and T97, were compared. It was concluded that the late-time data ISP was most appropriate, and it was found to reduce the seemingly enhanced levels of polarization on April 30 (+34 d) near 6600 Å seen in SN 1993J compared to other Type IIb SNe (Chornock et al. 2011).

The ISP removed data revealed a significant dominant axis on the q − u planes (particularly at early epochs), suggesting the presence of significant continuum polarization. Unfortunately, due to the limited availability of spectral regions devoid of strong lines, we were not able to directly measure the continuum polarization, which is the main probe of global asymmetry.

Significant line polarization of He i and H α was observed, including HV He i λ5876 and low amplitude polarization peaks to the blue of the main H α feature. syn++ fits showed that either HV H α or Si ii could fit the most prominent of these peaks, but the continuity of the rotation of the P.A. between these peaks and the main H α polarization features suggests they are related. SN 1993J is one of the few stripped-envelope SNe to have HV He i polarization features, alongside iPTF13bvn (Reilly et al. 2016), and it is the only one to show potential HV H α features. The reason for this is not clear from our analysis.

Overall we find that the continuum and line polarization features observed can be interpreted as a superposition of anisotropically distributed line forming regions on top of ellipsoidal ejecta geometry or an off-centre energy source within a spherical photosphere (see Fig. 12 – Section 5.3).

As emphasized in previous studies (e.g. Mauerhan et al. 2015; Stevance et al. 2019), modelling will be required to further interpret the library of spectropolarimetric observations that has been gathered for stripped-envelope CCSNe, and in particular the Type IIbs. Since we now know that toy models previously employed (Maund et al. 2007a; Reilly et al. 2016, 2017) can suffer significant degeneracies that can lead to misinterpretation or inconclusive results (Stevance 2019), we refrain from using such an approach. In the future, the new ISP corrected data presented here can be the object of further study using 3D radiative transfer models.

ACKNOWLEDGEMENTS

We are grateful to H. Tran for agreeing sharing their data. We thank Jason Spyromilio for clarifying the use of a reference and Kerry Jones for stylistic advice. We thank Luc Dessart for promptly answering a question we had regarding Dessart et al. (2018). HFS was supported through a PhD scholarship granted by the University of Sheffield during most of this project. HFS also acknowledges the support of the Marsden Fund Council managed through Royal Society Te Aparangi. PJV is supported by the National Science Foundation Graduate Research Fellowship Program Under Grant No. DGE-1343012. JCW is supported by the NSF grant AST-1813825. The following packages were used for the data reduction and analysis: matplotlib (Hunter 2007), astropy (The Astropy Collaboration 2018), numpy, scipy, and pandas (Jones et al. 2017). We also want to highlight our use of gelato.⁴

Footnotes

REFERENCES