Three-Dimensional Explosion Geometry of Stripped-Envelope Core-Collapse Supernovae. II. Modelling of Polarization

TL;DR

This study models polarization in stripped-envelope supernovae to reveal that their ejecta have a three-dimensional clumpy structure, which explains observed polarization loops and suggests large-scale convection or instabilities during explosion.

Contribution

It demonstrates that 3D clumpy models are necessary to reproduce polarization features, indicating supernova ejecta are inherently three-dimensional and large-scale.

Findings

3D clumpy structures reproduce polarization loops

Large clumps (>25% of photospheric radius) are inferred from polarization data

Supernova ejecta likely have large-scale convection or instabilities

Abstract

We present modelling of line polarization to study multi-dimensional geometry of stripped-envelope core-collapse supernovae (SNe). We demonstrate that a purely axisymmetric, two-dimensional geometry cannot reproduce a loop in the Stokes Q-U diagram, i.e., a variation of the polarization angles along the velocities associated with the absorption lines. On the contrary, three-dimensional (3D) clumpy structures naturally reproduce the loop. The fact that the loop is commonly observed in stripped-envelope SNe suggests that SN ejecta generally have a 3D structure. We study the degree of line polarization as a function of the absorption depth for various 3D clumpy models with different clump sizes and covering factors. Comparison between the calculated and observed degree of line polarization indicates that a typical size of the clump is relatively large, >~ 25 % of the photospheric radius. Such large-scale clumps are similar to those observed in the SN remnant Cassiopeia A. Given the small size of the observed sample, the covering factor of the clumps is only weakly constrained (~ 5-80 %). The presence of large-scale clumpy structure suggests that the large-scale convection or standing accretion shock instability takes place at the onset of the explosion.

Tables

Table 1. Table 1 Summary of the models

Model	${\tau }_{\mathrm{ph}}$a	${\alpha }_{\mathrm{cl}}$b	$f_{\mathrm{cl}}$c
2D-bipolar-30degd	10.0	–	0.13
2D-torus-20dege	10.0	–	0.35
3D-a0.5-f0.3	3.0,10.0,30.0,100.0	0.5	0.3
3D-a0.25-f0.3	3.0,10.0,30.0,100.0	0.25	0.3
3D-a0.125-f0.3	3.0,10.0,30.0,100.0	0.125	0.3
3D-a0.5-f0.06	30.0	0.5	0.06
3D-a0.5-f0.2	30.0	0.5	0.2
3D-a0.5-f0.5	10.0	0.5	0.4
3D-a0.5-f0.7	10.0	0.5	0.7

Equations

$$\bm{I}=\left(\begin{array}[]{c}I\\ Q\\ U\end{array}\right)=\left(\begin{array}[]{c}1\\ 0\\ 0\end{array}\right).$$

$${ϵ}_{\mathrm{out}}={ϵ}_{\mathrm{in}}\frac{1-{\mu }_{\mathrm{in}}v/c}{1-{\mu }_{\mathrm{out}}v/c},$$

$${\lambda }_{\mathrm{out}}={\lambda }_{\mathrm{in}}\frac{1-{\mu }_{\mathrm{out}}v/c}{1-{\mu }_{\mathrm{in}}v/c},$$

$$\bm{R}(\Theta)=\frac{3}{4}\left(\begin{array}[]{ccc}\cos^{2}\Theta+1&\cos^{2}\Theta-1&0\\ \cos^{2}\Theta-1&\cos^{2}\Theta+1&0\\ 0&0&2\cos\Theta\end{array}\right),$$

$$\bm{L}(\phi)=\left(\begin{array}[]{ccc}1&0&0\\ 0&\cos 2\phi&\sin 2\phi\\ 0&-\sin 2\phi&\cos 2\phi\\ \end{array}\right).$$

$${\mathbf{I}}_{\mathrm{out}}=\mathbf{L}(\pi -i_2)R(\Theta )\mathbf{L}(-i_1){\mathbf{I}}_{\mathrm{in}}.$$

$$p.d.f=\frac{1}{2}({\cos }^2\Theta +1)+\frac{1}{2}({\cos }^2\Theta -1)(\cos 2i_1{Q}_{\mathrm{in}}/I_{\mathrm{in}}-\sin 2i_1{U}_{\mathrm{in}}/I_{\mathrm{in}}).$$

Full text

1. Introduction

Core-collapse supernovae (SNe) are the explosions of massive stars. Since core-collapse SNe eject heavy elements synthesized inside of the stars, they play vital roles in chemical enrichment of galaxies. In addition, because of the large kinetic energy of the explosion ($\sim 10^{51}$ erg), SNe are also important for the galaxy formation. Despite their importance, the mechanism of the core-collapse SNe is a long-standing mystery (see Janka 2012; Kotake et al. 2012; Burrows 2013; Müller 2016 for reviews). Results of numerical simulations agree to the point that massive stars would not explode in one dimensional simulations (Rampp & Janka, 2000; Liebendörfer et al., 2001; Thompson et al., 2003; Sumiyoshi et al., 2005) except for some cases of the least massive stars (Kitaura et al., 2006; Janka et al., 2008). Therefore, multi-dimensional effects or deviation from spherical symmetry are believed to be crucial for successful explosions.

The leading scenario of core-collapse SNe is neutrino-driven explosion, where multi-dimensional effects can appear by convection (e.g., Herant et al., 1994; Burrows et al., 1995; Janka & Mueller, 1996) or standing accretion shock instability (SASI, e.g., Blondin et al., 2003; Scheck et al., 2004; Ohnishi et al., 2006; Foglizzo et al., 2007; Ott et al., 2008; Iwakami et al., 2008; Fernández, 2010; Hanke et al., 2012). In fact, some successful explosions have been reported by two-dimensional (2D) or three-dimensional (3D) simulations (e.g., Buras et al., 2006; Marek & Janka, 2009; Suwa et al., 2010; Müller et al., 2012; Takiwaki et al., 2012; Hanke et al., 2013; Bruenn et al., 2013; Takiwaki et al., 2014; Couch & O’Connor, 2014; Melson et al., 2015; Lentz et al., 2015; Roberts et al., 2016; Müller, 2016) although the obtained explosion energy is usually lower than $10^{51}$ erg. Another scenario is magneto-rotational explosion (e.g., Yamada & Sawai, 2004; Kotake et al., 2004; Takiwaki et al., 2009; Sawai et al., 2005; Burrows et al., 2007; Obergaulinger et al., 2009; Mösta et al., 2014), where the amplified magnetic fields drive the explosion. In this scenario, the bipolar explosion is generally expected.

In order to link these theoretical models with observations, it is necessary to derive multi-dimensional geometry from observed SNe. The most straightforward method is spatially-resolved observations of nearby SN remnants (e.g., Hwang et al., 2004; Isensee et al., 2010; DeLaney et al., 2010; Milisavljevic & Fesen, 2015). However, the number of accessible objects is limited. To advance our knowledge, it is therefore important to study the multi-dimensional geometry of extragalactic SNe. In fact, many efforts have been made to derive the multi-dimensional geometry from extragalactic SNe, for example by using spectral line profiles at late phases (\mathrel{\hbox to0.0pt{\lower 4.0pt\hbox{\hskip 1.0pt\sim}\hss}\raise 1.0pt\hbox{>}} 1 yr after the explosion, e.g., Spyromilio 1994; Sollerman et al. 1998; Matheson et al. 2000; Mazzali et al. 2001; Maeda et al. 2002; Elmhamdi et al. 2004; Mazzali et al. 2005; Maeda et al. 2008; Modjaz et al. 2008; Tanaka et al. 2009b; Taubenberger et al. 2009; Maurer et al. 2010; Chornock et al. 2010; Valenti et al. 2011; Shivvers et al. 2013; Roy et al. 2013; Chen et al. 2014; Milisavljevic et al. 2015; Mauerhan et al. 2016)

Polarization at early phases (\mathrel{\hbox to0.0pt{\lower 4.0pt\hbox{\hskip 1.0pt\sim}\hss}\raise 1.0pt\hbox{<}} 50 days after the explosion) is one of the most powerful methods to derive multi-dimensional geometry from extragalactic SNe (see Wang & Wheeler 2008 for a review). By observations, we can measure continuum and line polarizations. In the SN ejecta, electron scattering is the dominant source of polarization. Line scattering generally produces less polarization (Howell et al. 2001; Kasen et al. 2003, see also Jeffery 1989), and it is often assumed that line scattering works as a de-polarizer. From the spherical symmetric SN ejecta, no polarization should be detected because of complete cancellation of polarization vectors. Non-zero continuum polarization would be observed when the photosphere deviates from spherical symmetry (Shapiro & Sutherland, 1982; Höflich, 1991; Höflich et al., 1996; Dessart & Hillier, 2011; Bulla et al., 2015). In addition, even for the spherical photosphere, non-zero line polarization would be observed when the distribution of an ion producing the corresponding absorption line is not spherical symmetry (Kasen et al., 2003; Hole et al., 2010). Therefore, line polarization can be a diagnostic to multi-dimensional element distribution in the SN ejecta.

In this paper, we present modelling of line polarization in stripped-envelope SNe (SNe of Type IIb, Ib, and Ic) to obtain connections between the polarization properties and the element distribution in the SN ejecta. Compared with the cases of H-rich SNe, closer insight on the explosion mechanism can be obtained for stripped-envelope SNe, as the large hydrogen envelope is not present. In Section 2, we describe our method to compute polarization signature of the SN models. In Sections 3 and 4, we show results of 2D and 3D models, respectively. We discuss implication of our results in Section 5 and give summary in Section 6.

2. Methods

2.1. Radiation Transfer

We perform 3D radiation transfer simulations to study the properties of the line polarization. For this purpose, we use a simple 3D Monte Carlo radiation transfer code. The code takes into account the electron scattering and the line scattering. We treat only a single line at a single epoch rather than modelling time evolution of full spectra since we aim to obtain the connection between explosion geometry and properties of line polarization (see Hole et al. 2010 for a similar strategy). The code computes the polarization spectrum of the line for arbitrary 3D distribution of the line optical depth. More details of the code are given in Appendix.

We use $100\times 100\times 100$ linearly distributed Cartesian meshes. The velocity is used as a spatial coordinate thanks to the homologous expansion ($r\propto v$). The maximum velocity is 25000 km s*-1*, and thus the resolution is 500 km s*-1*, giving the resolution of $\lambda/\Delta\lambda=c/\Delta v=600$, which is comparable to a typical spectral resolution of low-resolution spectropolarimetric observations.

We start simulations by generating unpolarized photon packets from the spherical inner boundary ($v=v_{\rm in}$). The electron scattering optical depth from the inner boundary to infinity is set to $\tau_{\rm in}$. In this paper, we adopt $\tau_{\rm in}=3$ as in Kasen et al. (2003) and Hole et al. (2010). Note that the photosphere ($v=v_{\rm ph}$) is defined as the position where the electron scattering optical depth is unity, and thus, the inner boundary of the computation is located inside of the photosphere.

The photon packets are then tracked by taking into account the electron scattering and the line scattering. For the electron scattering, we use a power-law electron density profile, $n_{e}\propto r^{-n}$. The electron density is assumed to be spherically symmetric. For the power-law index, we use $n=7$ which describes the line forming region of hydrodynamic models of stripped-envelope SNe (Iwamoto et al., 2000; Mazzali et al., 2000). Although the very outermost ejecta has a steeper slope ($n\sim 10$, Matzner & McKee 1999), we use a single power-law profile since the outermost ejecta does not have a strong contribution to absorption lines. Note that polarization pattern is not affected by the slope if the slope is steep enough (n\mathrel{\hbox to0.0pt{\lower 4.0pt\hbox{\hskip 1.0pt\sim}\hss}\raise 1.0pt\hbox{>}}5, Kasen et al. 2003). We assume the photospheric velocity ($v_{\rm ph}$) and the time after the explosion ($t$), which give the photospheric radius $r_{\rm ph}=v_{\rm ph}t$. Then, with the condition that the electron scattering optical depth is unity at the photosphere, the normalization of the electron density is obtained. We adopt $v_{\rm ph}=8,000$ km s*-1* and $t=20$ days as typical values for stripped-envelope SNe around the maximum light.

For the line scattering, we use Sobolev approximation (Castor, 1970), which is a sound approximation in the SN ejecta with a large velocity gradient. For the Sobolev line optical depth, we assume a power-law radial profile above the photosphere, $\tau_{\rm line}=\tau_{\rm ph}(r/r_{\rm ph})^{-n}$. Here $\tau_{\rm ph}$ is the Sobolev optical depth at the photosphere. For simplicity, we use the same power-law index with the electron density ($n=7$).

In addition to the spherical component of the line optical depth, we assume an enhancement by a factor of $f_{\rm\tau}$ in some regions e.g., a torus or clumps. Note that this is different from the treatment by Hole et al. (2010), where the line opacity is set to be zero outside of the clumps. Such a treatment seems more suitable for Type Ia SNe (as they applied for), where a strong line is formed dominantly in a certain layer, e.g., Si lines are produced mostly in the Si-rich layers. On the other hand, for strong lines in stripped-envelope core-collapse SNe, such as those of Ca and Fe, both pre-SN and newly synthesized elements contribute to the absorption. Therefore, we assume an enhancement in addition to the spherical component. In the models presented in this paper, we adopt $f_{\rm\tau}=10.0$. Implication of this choice is discussed in Section 4. The parameters for the models are summarized in Table 1.

2.2. Comparison with Observations

We study the explosion geometry of stripped-envelope SNe by comparing results of our simulations with observations. Figure 1 shows an example of spectropolarimetric data of stripped-envelope SNe (Type Ib SN 2009jf, Tanaka et al. 2012). In this paper, we define Stokes parameters as a fraction of the total flux: $Q\equiv\hat{Q}/I$ and $U\equiv\hat{U}/I$, where $\hat{Q}$ and $\hat{U}$ are polarized fluxes, i.e., $\hat{Q}=I_{0}-I_{90}$ and $\hat{U}=I_{45}-I_{135}$, respectively ($I_{\psi}$ is the intensity measured through the ideal polarization filter with an angle $\psi$). From Stokes parameters $Q$ and $U$, the position angle of the polarization, $\theta$, is obtained by $2\theta={\rm atan}(U/Q)$.

Properties of line polarization in stripped-envelope SNe can be summarized as follows:

Non-zero line polarization is common and polarization feature shows an inverted P-Cygni profile which peaks at flux absorption minimum (e.g., Kawabata et al., 2002; Leonard et al., 2002; Wang et al., 2003b; Maund et al., 2007a, b, 2009; Tanaka et al., 2008, 2009a, 2012; Mauerhan et al., 2015; Stevance et al., 2016; Mauerhan et al., 2016). 2. 2.

When the polarization data across the line (middle panel of Figure 1) are plotted in the Stokes $Q$-$U$ diagram (right panel), the observed data commonly show a loop (e.g., Maund et al., 2007a, b, 2009; Tanaka et al., 2012; Mauerhan et al., 2015; Stevance et al., 2016; Mauerhan et al., 2016). 111Such a loop in the $Q-U$ diagram has also been observed in Type Ia SNe (e.g., Wang et al., 2003a; Kasen et al., 2003; Chornock & Filippenko, 2008; Patat et al., 2009; Tanaka et al., 2010; Porter et al., 2016; Milne et al., 2016) and Type IIn SNe (e.g., Hoffman et al., 2008) as well as in Wolf-Rayet stars (e.g., Schulte-Ladbeck et al., 1990; St-Louis et al., 2012). 3. 3.

The degree of line polarization, i.e., the maximum polarization level at the absorption lines, is generally a few percent, and tends to be higher for stronger lines (Tanaka et al., 2012).

The degree of polarization depends on the strength of the absorption. The absorption strength is mainly determined by global properties of SNe, such as ejecta mass, temperature and element abundances, and not directly by the explosion geometry. Therefore, it is important to compare features with a similar absorption strength to discuss the explosion geometry. Here we define a fractional depth (FD) of absorption at the absorption minimum, FD $=(f_{\rm cont}-f_{\rm abs})/f_{\rm cont}$, where $f_{\rm abs}$ and $f_{\rm cont}$ are the flux at the absorption minimum and at the continuum near the absorption line, respectively. Tanaka et al. (2012) showed that, in a simple configuration, the observed polarization ($P_{\rm obs}$) can be approximately described as $P_{\rm obs}\simeq P_{\rm cor}[{\rm FD}/({\rm 1-FD})]$, where a corrected polarization $P_{\rm cor}$ is defined as the polarization level if the fractional depth would be FD $=0.5$.

3. Results: 2D Models

We first study polarization properties of 2D axisymmetric models. As 2D models, we construct a bipolar model and a torus model. These models are motivated by the results of nucleosynthesis calculations for 2D bipolar (or jet-like) explosion models (e.g., Nagataki et al., 1997; Maeda et al., 2003; Nagataki et al., 2006; Tominaga, 2009). In these models, explosively synthesized elements such as Fe are preferably produced in the polar region, Our 2D bipolar model depicts such a case. In contrast, the elements produced mainly in the pre-SN stage such as O may be distributed in a torus-like geometry, which are represented by our torus model.

Polarization properties of the bipolar model are shown in the top panels of Figure 2. The half opening angle of the polar blobs is set to be 30 deg. As the opacity distribution is not spherically symmetric, non-zero polarization appears. The polarization data in the $Q-U$ diagram shows a straight line. This is always the case for every line of sight. The position in the $Q-U$ diagram represents the position angle, i.e., $\theta=(1/2){\rm atan}(U/Q)$. Therefore, the straight line in the $Q-U$ diagram means a constant position angle across the P-Cygni profile.

The observed position angle can be rotated depending on the direction of the symmetric axis of the model on the sky. However, as long as the 2D bipolar structure is kept, the polarization always shows a straight line in the $Q-U$ diagram. Also, this behavior does not depend on global parameters such as the optical depth at the photosphere ($\tau_{\rm ph}$), the enhancement factor ($f_{\rm\tau}$) since this behavior is purely caused by the geometric effect. We also test the models with different sizes of the blobs (i.e., opening angles of the bipolar structure), and confirm that, although the number of lines of sight to have a high polarization degree depends on the size of the blob, the straight line in the Q-U plane is always obtained.

Similar polarization properties are obtained for the torus model, i.e., polarization always shows a straight line in the $Q-U$ diagram. The bottom panels of Figure 2 show an example of the results for the torus model with a half opening angle of 20 deg. Note that, for a certain line of sight, a 90 deg rotation in the position angle can be observed. For example, for the line of sight of 60 deg from the pole as shown in Figure 2, the lateral part of the photospheric disk is hidden near the photospheric velocity while the bottom part of the photospheric disk is hidden at higher velocities. As a result, a positive Stokes $Q$ is obtained near the photospheric velocity while a negative Stokes $Q$ is obtained at higher velocities (middle panel). This corresponds to a 90 deg rotation in the position angle. However, only a 90 deg rotation can occur as long as the underlying model keeps axisymmetry since there is no way to produce Stokes $U$ component if the axisymmetric angle of the model is set to be north ($\theta=0$ deg) as shown in Figure 2. If the symmetric axis of the model is rotated on the sky, Stokes $U$ components can appear but the polarization data still form a straight line in the $Q-U$ diagram.

In summary, we validate the statement commonly made by previous works (e.g., Kasen et al., 2003; Wang & Wheeler, 2008), i.e., a purely axisymmetric element distribution cannot reproduce the loop in the $Q-U$ diagram. When the element distribution has a purely 2D axisymmetric structure such as bipolar blobs or a torus, the polarization shows a straight line in the $Q$-$U$ diagram.

4. Results: 3D Models

Next we study polarization properties of 3D models. Motivated by 3D simulations of neutrino-driven explosions, where various sizes of complex structure appear, we set up 3D models by randomly placing different numbers of spherical clumps with different sizes. Here we introduce two parameters to depict the model: the size parameter of the clumps $\alpha_{\rm cl}$, i.e., the radius of the clump normalized by the photospheric radius ($\alpha_{\rm cl}=v_{\rm cl}/v_{\rm ph}$) and the photospheric covering factor ($f_{\rm cl}$). Since the optical depth near the photosphere is the most important for line formation, the covering factor is evaluated by taking into account the clumps only in a shell between $v=v_{\rm ph}$ and $v_{\rm ph}+2000$ km s*-1*. Note that, as in the 2D cases, the line optical depth in our models has a spherical component and it is enhanced by a factor of $f_{\tau}=10.0$ within the clumps.

The top panels of Figure 3 show the polarization properties of the 3D model with the clump size of $\alpha_{\rm cl}=0.5$ and the covering factor of $f_{\rm cl}=0.3$. In the polarization spectrum, both Stokes $Q$ and $U$ parameters vary across the lines (middle panel), and polarization shows a loop in the $Q-U$ diagram (right panel), as also found by Hole et al. (2010).

The $Q-U$ loop in the 3D clumpy models can be understood as follows. In the 3D clumpy models, depending on the Doppler velocities, different parts of the photospheric disk are hidden by the clumps. Since the distribution of the clumps does not have a common symmetric axis, the position angle of the polarization can change depending on the Doppler velocities. In general, the change in the position angle across the line can be arbitrary large, i.e., the polarization in the $Q-U$ diagram can be scattered around. But for the relatively large size of the clumps as in the case of $\alpha_{\rm cl}=0.5$, the same clump keeps contributing to the absorption even for different Doppler velocities, and thus the change in the position angle tends to be smooth as a function of Doppler velocities. Therefore, the polarization tends to show a loop in the $Q-U$ diagram in the 3D clumpy distribution with relatively large clumps. Note that the $Q-U$ loop can also be produced by other geometries, e.g., a combination of the ellipsoidal photosphere and ellipsoidal line scattering shell whose symmetric axes are misaligned with each other (Kasen et al., 2003), But even in such a case, it is required that the axisymmetry of the system is broken.

4.1. Size of the clumps

The size of the clumps is of interest to study the origin of the 3D structure in the SNe. We show the first attempt to quantify the size of the clumps by comparing the results of the modelling and the observed polarization degrees, i.e., the maximum polarization level at the absorption line. We calculate the polarization spectra with different sizes of the clumps by keeping the covering factor of $f_{\rm cl}=0.3$ and other parameters to be the same. The middle and bottom panels in Figure 3 show the results for the 3D models with $\alpha_{\rm cl}=0.25$ and $0.125$, respectively.

As shown in the figures, for a given covering factor, models with smaller clumps show a lower polarization. In such models, the photospheric disk is hidden by many small clumps, and polarization vectors tend to be cancelled out (Figure 4). This behavior was also pointed out by Hole et al. (2010) in the context of Type Ia SNe. Since stripped-envelope SNe generally show non-zero line polarization, the typical size of the clumps should not be too small.

Since the polarization degree depends not only on the geometry but also on the absorption strength, it is important to compare the models and observations for similar absorption strengths. Therefore, in Figure 5, we compare models and observations in the plane of the polarization degree and the fractional absorption depth. The black points with error bars are observational data of the Ca ii (filled) and Fe ii (open) lines for 6 Type Ib and Ic SNe analyzed in Tanaka et al. (2012). The small dots show the polarization degree of the models for 100 lines of sight. In each panel, we show four sets of the models with the same size and distribution of the clumps but with the different line optical depth at the photosphere ($\tau_{\rm ph}$ = 3.0, 10.0, 30.0, and 100.0 from left to right).

When the clump is as small as $\alpha_{\rm cl}=0.125$ (bottom panels of Figures 3 and 5), the polarization degree cannot be $>0.5\%$ for any line of sight. For the larger sizes of the clumps, a higher polarization can be obtained. When the size of the clumps is $v_{\rm cl}=0.25$ (middle panels), the polarization degrees of these models are still short of some of the observed polarization. When the size of the clumps is relatively large, $\alpha_{\rm cl}=0.5$ (top panels), the polarization degree can be as high as $>1\%$ for the fractional depth of 0.5.

Ideally the polarization properties of the models should be compared with the statistical distribution of the observed polarization. Although the number of objects with good data is still small, Figure 6 shows a cumulative distribution of polarization properties of 6 Type Ib and Ic SNe. To define one characteristic polarization for each object, we take the average of the corrected polarization ($P_{\rm cor}$) for the Ca ii and Fe ii lines. Color lines show the cumulative distribution of the modelled polarization for 100 lines of sight. We choose models with $\tau_{\rm ph}=30.0$, which approximately give $FD\sim 0.5$ (Figure 5).

The comparison in the cumulative distribution clearly shows that the model with too small clumps ($\alpha_{\rm cl}=0.125$) is not consistent with the observations. The $p$ value for a Kolmogorov-Smirnov (KS) test is $p_{\rm KS}=0.0016$. Since the number of objects is so small, we cannot distinguish the model with the clump size of $\alpha_{\rm cl}=0.25$ ($p_{\rm KS}=0.54$) and $\alpha_{\rm cl}=0.5$ ($p_{\rm KS}=0.94$). Nevertheless, the model with $\alpha_{\rm cl}=0.25$ is already short to explain the polarization level of $>1\%$, and seems to close to the lower limit for the clump size. Here it is noted that our models adopt an enhancement factor of $f_{\rm\tau}=10.0$. For a higher enhancement factor, the polarization degree is not largely affected because models with $f_{\rm\tau}=10.0$ already give an optically thick absorption in the clumps near the photosphere. On the other hand, for a smaller enhancement factor, the polarization degree decreases for a given FD. In such cases, even larger clumps is required to reproduce a high polarization degree. Therefore, we conclude that a typical size of the 3D clumps should be \mathrel{\hbox to0.0pt{\lower 4.0pt\hbox{\hskip 1.0pt\sim}\hss}\raise 1.0pt\hbox{>}}25\% of the photospheric radius to reproduce observed polarization degrees.

4.2. Covering factor of the clumps

To obtain possible constraints on the number or the covering factor of the clumps in the ejecta, we vary the covering factors of clumps keeping their size to be $\alpha_{\rm cl}=0.5$. Figure 7 shows the model input (left) and cumulative distributions of the resultant polarization (right). For the models, we choose the line strength at the photosphere to have $FD\sim 0.5$, i.e., $\tau_{\rm ph}$ = 30 for the models with $f_{\rm cl}=0.06$, 0.2, and 0.3. $\tau_{\rm ph}$ = 10 for the models with $f_{\rm cl}=0.4$ and 0.7. The observed distribution is the same as in Figure 6.

For the model with a smaller covering factor ($f_{\rm cl}=0.06$), the probability to have a high polarization is also low. Then, by increasing the covering factor of the clumps, a higher polarization can be more frequently observed ($f_{\rm cl}=0.2-0.5$). However, if the covering factor of the clumps is too large ($f_{\rm cl}=0.7$), the distribution of resultant polarization shifts toward a lower value again since the system restores the symmetry again.

Since the observational samples are small, it is difficult to draw a firm conclusion on the covering factor of the clumps. However, the models with $f_{\rm cl}=0.06$ and $f_{\rm cl}=0.7$ are already at the edge of the distribution. By taking into account the fact that models with $f_{\rm\tau}=10.0$ tend to give an upper limit of the polarization level (see Section 4.1), it seems that current data do not support models with too small covering factors (f_{\rm cl}\mathrel{\hbox to0.0pt{\lower 4.0pt\hbox{\hskip 1.0pt\sim}\hss}\raise 1.0pt\hbox{<}}0.05) and too large covering factors (f_{\rm cl}\mathrel{\hbox to0.0pt{\lower 4.0pt\hbox{\hskip 1.0pt\sim}\hss}\raise 1.0pt\hbox{>}}0.8).

5. Discussion

We have modelled line polarization of stripped-envelope core-collapse SNe. The results of modelling are summarized as follows. (1) The observed $Q-U$ loop cannot be explained by the 2D axisymmetric models, but can be explained by the 3D clumpy models. (2) By comparing the results of the 3D clumpy models with the observed degrees of line polarization, it is found that a typical size of the clumps is relatively large, i.e., \mathrel{\hbox to0.0pt{\lower 4.0pt\hbox{\hskip 1.0pt\sim}\hss}\raise 1.0pt\hbox{>}}25\% of the radius, and a covering factor of the clumps in the ejecta is not too small and not too large ($5-80\%$).

It is intriguing that such a large-scale clumpy structure is also seen in the element distribution of Cassiopeia A (e.g., Hwang et al., 2004; Isensee et al., 2010; DeLaney et al., 2010; Grefenstette et al., 2014; Milisavljevic & Fesen, 2015), which is a supernova remnant produced by a Type IIb SN (Krause et al., 2008). The similarity suggests that the element distribution as seen in Cassiopeia A may also be able to reproduce the polarization properties observed in early phase of SNe. Here we discuss possible origins for the clumpy structure suggested by observations and modelling.

One scenario is the Rayleigh-Taylor (RT) instability, causing matter mixing in the SN ejecta. By the RT instability, many clumps are produced and metal-rich ejecta inside are delivered toward the outer layers (e.g., Hachisu et al., 1990; Fryxell et al., 1991; Herant & Woosley, 1994; Nagataki et al., 1998; Kifonidis et al., 2003; Joggerst et al., 2010; Ono et al., 2013). However, the RT instability alone usually produces small fingers in many directions. This is similar to the case of $\alpha_{\rm cl}=0.125$ in Figure 4 and not consistent with the observations.

The clumpy structure suggested by observations is more in favor of large-scale convection or SASI developed in the initial stage of the explosion. When the large-scale convection or SASI takes place, the subsequent evolution of the shock becomes asymmetric, which produces the large-scale asymmetry in the element distribution (e.g., Kifonidis et al., 2003, 2006; Hammer et al., 2010; Fujimoto et al., 2011). Also, results of long-term simulations show that the ejecta structure near the shock breakout still keeps an imprint of the large-scale asymmetry generated by neutrino-driven convection and SASI, with the small-scale structures by the RT instability added on top of it (Wongwathanarat et al., 2013, 2015). Note that such long-term simulations for neutrino-driven explosion also nicely reproduce the geometry of Cassiopeia A (Wongwathanarat et al., 2016).

It is worth noting that, although the loop in the $Q-U$ diagram does not support a purely axisymmetric element distribution (Figure 2), spectropolarimetric data do not rule out the presence of an overall bipolar structure or a dominant axis in the SN ejecta. As long as some large-scale, non-axisymmetric components exist, they can produce a large enough polarization level and the loop in the $Q-U$ diagram. In fact, analysis of the [O i] line profiles in the late-phase spectra suggest a torus-like distribution of oxygen, which is consistent with a bipolar explosion (e.g., Maeda et al., 2008; Modjaz et al., 2008; Tanaka et al., 2009b). Since polarization at the early phase and nebular line profile are sensitive to the outer and inner ejecta, respectively, the combination of early and late phase observations may indicate that global 2D structure exists more in the inner ejecta and 3D clumpy structure is added in the outer ejecta. It is noted that, even by the late phase observations, presence of clumpy structure has also been suggested by the studies of line profiles (e.g., Spyromilio, 1994; Sollerman et al., 1998; Taubenberger et al., 2009), ionization states (e.g., Mazzali et al., 2007b, a), and dust (e.g., Sugerman et al., 2006; Ercolano et al., 2007; Kotak et al., 2009; Wesson et al., 2015; Dwek & Arendt, 2015; Bevan & Barlow, 2016). Thus, the transition from the inner 2D to the outer 3D structures may be somewhat gradual. Interestingly, our studies suggest that the shape of the loop in the $Q-U$ diagram can be used as a probe of such a combined (2D $+$ 3D) geometry. As expected from the results of 2D (Figure 2) and 3D models (Figure 3), if the ejecta has an overall 2D geometry + 3D clumpy structure, it tends to produce an elongated loop in the $Q-U$ diagram. Although current observational data do not allow us to extract such information, detailed studies will be possible in future with more observational samples with high signal-to-noise ratio.

It is emphasized that our modelling includes many simplifications. For example, we parameterize the line optical depth and enhancement factor, but they must be determined by the combination of element abundance, temperature, and ionization states. Thus, our models shown in the left panels of Figures 2 and 3 are not readily connected with the element distribution. Full radiation transfer modelling using 3D hydrodynamic models is required to obtain a closer link between the explosion models and observations. Also, the comparison with observed polarization degree is done by averaging the polarization degrees of different lines. Since polarization at different absorption lines reflects the distribution of each element and ion, the direct comparison for each element is necessary when larger observed samples and full transfer calculations are available.

6. Summary

We have performed 3D radiation transfer simulations for the analysis of line polarization in stripped-envelope SNe. We demonstrate that a purely axisymmetric, 2D structure always produces a straight line in the Stokes $Q-U$ diagram, and cannot explain the commonly observed loop in the $Q-U$ diagram. On the contrary, 3D clumpy structures naturally reproduce the loop. Comparison of the results of the modelling and the observed polarization degrees enables to constrain a typical size of the clumps from polarization data for the first time. To reproduce the distribution of the observed polarization degrees (0.5-2.0 %), a typical size of the clump should be relatively large, i.e., $>25\%$ of the photospheric radius (or the radius where the clump is located). The covering factor of the clump in the ejecta is only weakly constrained i.e., to $5-80$ %. Such a large-scale clumpy structure inferred by polarization is similar to that seen in the SN remnant Cassiopeia A.

The large-scale clumpy structure is unlikely to be produced only by the RT instability as it tends to produce small fingers in many directions. Instead, the presence of the large-scale clumpy structure in the ejecta suggests that large-scale convection or SASI takes place at the onset of the explosion. Polarization properties do not necessarily exclude the presence of a dominant axis in the SN ejecta since non-axisymmetric structure on top of the 2D axisymmetric structure can also reproduce the loop in the $Q-U$ diagram. In fact, the analysis of the nebular spectra supports a bipolar geometry in the innermost layer. These observational constraints suggest that SN ejecta may have an overall 2D bipolar structure inside and 3D clumpy structure outside. We speculate that such a hybrid structure could be produced by SASI. In order to obtain further constraints on the explosion mechanism, polarization modelling using realistic SN models will be worthwhile as more and more long-term realistic simulations from core collapse to the shock breakout are becoming available.

We thank Takashi Hattori, Kentaro Aoki, Masanori Iye, Elena Pian, Toshiyuki Sasaki, and Masayuki Yamanaka for their contribution to the spectropolarimetric observations with the Subaru telescope, and the referee for valuable comments. MT thanks Thomas Janka, Takashi Moriya, and Takaya Nozawa for fruitful discussion. Numerical simulations presented in this paper were carried out with Cray XC30 at Center for Computational Astrophysics, National Astronomical Observatory of Japan. This research has been partly supported by the Grant-in-Aid for Scientific Research from JSPS (24740117, 26800100, 15H02075) and MEXT (25103515, 15H00788), and by World Premier International Research Center Initiative (WPI Initiative), MEXT, Japan.

Appendix A Three-Dimensional Radiation Transfer Code

We have developed a new 3D radiation transfer code to compute polarization spectrum of one line from arbitrary 3D distribution of the line optical depth. The code uses the Monte Carlo method, which is a common method to compute polarization by scattering processes (e.g., Daniel, 1980; Whitney & Hartmann, 1992; Hillier, 1991; Code & Whitney, 1995; Whitney, 2011). For the application to SNe, see Höflich (1991); Kasen et al. (2003, 2006); Dessart & Hillier (2011).

A.1. Spatial and Wavelength Grid

We set up the 3D Cartesian spatial mesh with the $100\times 100\times 100$ meshes. The velocity is used as a spatial coordinate because the SN ejecta are homologously expands ($r\propto v$), The outer velocity of the grid is $v_{\rm max}=$ 25000 km s*-1*, and thus, the resolution is $\Delta v=$ 500 km s*-1*. This spatial resolution gives the wavelength resolution of $\lambda/\Delta\lambda=c/\Delta v=600$, which is sufficient to make comparison with observed data.

Since the code computes only one (arbitrary) line, the wavelength range used in the computation is very small. If the rest wavelength of the line is $\lambda_{0}$, we compute the spectrum only at the wavelength range between $\lambda_{0}(1-v_{\rm max}/c)$ and $\lambda_{0}(1+v_{\rm max}/c)$. Within this wavelength range, the energy spectrum is assumed to be constant ($\lambda F_{\lambda}$ = const).

A.2. Beginning of the Simulations

Our code assumes a sharply defined inner boundary, and solves radiation transfer above the boundary by tracking the photon packets in the expanding ejecta. Every photon packet has assigned energy, wavelength, and Stokes parameters. Especially each photon packet in the simulation has a constant energy, irrespective of the wavelength of the packet. Because of this treatment, any photon is not lost during the simulation, which results in the accurate energy conservation (Lucy, 1999; Kasen et al., 2006; Kromer & Sim, 2009).

The position of the inner boundary is determined so that the electron scattering optical depth from the inner boundary to infinity is $\tau_{\rm in}$. In the simulations used in the main text of the paper, we always adopt $\tau_{\rm in}=3$ (Table 1) as in Kasen et al. (2003) and Hole et al. (2010). The radiation from the inner boundary is assumed to be thermalized, and thus, to be unpolarized;

[TABLE]

The direction of the photon is determined by $\mu=\sqrt{z}$ (Mazzali & Lucy, 1993) (hereafter we use $z$ to denote a random number, $0<z\leq 1$), where $\mu$ is cosine of the angle between the radial and photon direction. The azimuthal angle around the radial direction $\psi$ is uniformly distributed, i.e., $\psi=2\pi z$.

A.3. Scattering Events

The emitted photon packets experience the electron scattering and the line scattering, which are treated in a similar way to that by Mazzali & Lucy (1993). For the electron scattering, we assume a power-law density structure with the power-law index $n$. We also have the photospheric velocity ($v_{\rm ph}$) and the epoch from the explosion ($t_{\rm d}$) as input parameters. The photospheric radius ($r_{\rm ph}=v_{\rm ph}t_{\rm d}$) is defined to be the radius where the optical depth for the electron scattering is unity. By setting $v_{\rm ph}$ and $t_{\rm d}$, the normalization of the electron density is determined.

For the line scattering, we use the Sobolev approximation (Castor, 1970), and assume a power-law optical depth profile with the same index $n$. The parameter for the line scattering is $\tau_{\rm ph}$, the optical depth at the photosphere. In addition, we assume enhancement of the optical depth by a factor of $f_{\tau}$ in some region. The parameters used in the simulations are summarized in Table 1.

A photon packet propagating in one computational grid can have 3 possible events; (1) escaping from the grid, (2) the electron scattering, and (3) the line scattering. The event that actually occurs is judged by calculating the length to the 3 events. It is simple to compute the length to the next grid $l_{\rm grid}$ for the given position and the direction vector of the photon packet. The direction to the electron scattering event is computed by the randomly selected event optical depth $\tau_{R}=-\ln(z)$. When the optical depth reaches this value, a scattering event occurs. Thus, the distance to the electron scattering $l_{\rm elec}$ can be computed by $\tau_{R}=n_{e}(\bm{r})\sigma l_{\rm elec}$. When $l_{\rm elec}$ is shorter than $l_{\rm grid}$, the electron scattering occurs if there is no contribution of line scattering.

Since the line scattering is treated as a resonance, the distance to the line scattering event is $l_{\rm line}=ct_{\rm d}(\lambda_{0}-\lambda^{\prime})/\lambda_{0}$, where $\lambda^{\prime}$ is the comoving wavelength of the photon packet. If $l_{\rm line}$ is shortest among 3 lengths, the line scattering is taken into account. The line scattering event actually occurs when the sum of the line scattering optical depth ($\tau_{\rm line}(\bm{r})$) and the electron scattering optical depth in $s_{\rm line}$ ($\tau_{e}=n_{e}(\bm{r})\sigma l_{\rm line}$) exceeds $\tau_{R}$. If this sum does not reach $\tau_{R}$, then the electron scattering opacity is evaluated and added again, and the fate of the packet is the electron scattering or the escape from the grid. For illustration of this process, see Figure 1 of Mazzali & Lucy (1993).

When the scattering event occurs, the next direction vector of the photon packets is determined. For the electron scattering, this scattering angle depends on the polarization, which is discussed in the next Section. For the line scattering, the direction is determined by the isotropic probability function in the comoving frame.

By the scattering event, the energy and the wavelength of the packet are changed. For the energy, by the energy conservation in the rest frame,

[TABLE]

where $\epsilon_{\rm in}$ and $\epsilon_{\rm out}$ are the rest-frame energy of the incoming and outgoing packets, respectively. And $\mu_{\rm in}$ and $\mu_{\rm out}$ are the cosines of the angles between the radial direction and incoming/outgoing propagating directions, respectively. Similarly, the change in the wavelength is given by

[TABLE]

where $\lambda_{\rm in}$ and $\lambda_{\rm out}$ are the rest-frame wavelength of the incoming and outgoing packet, respectively.

A.4. Polarization Calculations

Scattering events change the polarization properties of the photon packets. For the electron scattering, the phase matrix can be written as follows (Chandrasekhar, 1960);

[TABLE]

where $\Theta$ is the scattering angle on the plane of the scattering. This matrix should be operated in the scattering frame. In general, the rotation matrix for the Stokes parameters is written as follows (Chandrasekhar, 1960);

[TABLE]

By using these matrices, the effect on the Stokes parameters is given by

[TABLE]

Here $\bm{I}_{\rm in}$ and $\bm{I}_{\rm out}$ is Stokes parameter in the rest frame before and after the scattering, respectively. The angles $i_{1}$ and $i_{2}$ are the angles on the spherical triangle defined as in (Chandrasekhar 1960, see Figure 1 of Code & Whitney 1995).

Equation A6 means that the the angle-dependence of the intensity of the scattered light depends the polarization properties of the incident radiation. From Equation A6, the probability distribution function ($p.d.f$) of the total intensity is

[TABLE]

By using this function with the rejection method as outlined in Code & Whitney (1995), we determine the scattering angle of the electron scattering.

We assume that the line scattering works as a depolarizer, i.e., the emission turns into unpolarized state by the line scattering as assumed in previous studies (see Höflich et al., 1996; Kasen et al., 2006; Hole et al., 2010).

A.5. Test Calculations

For the computation of polarization for the electron scattering, the code was tested with the analytic formulae by Brown & McLean (1977) for optically thin cases, and with numerical results by Code & Whitney (1995) for optically thick cases. For both cases, we got an excellent agreement. For the application to a SN, expanding ejecta with the steep density slope, we checked our results with those by Kasen et al. (2003). We confirmed that our code gives the consistent results on the radial profile of polarization for several power-law indexes ($n$) and the inner boundaries ($\tau_{\rm in}$).