VLBI observations of supernova PTF11qcj: Direct constraints on the size of the radio ejecta

We present High Sensitivity Array (HSA) and enhanced Multi Element Remotely Linked Interferometer Network (eMERLIN) observations of the radio-loud broad-lined type Ic supernova PTF11qcj obtained $\sim7.5$ years after the explosion. Previous observations of this supernova at 5.5 yrs since explosion showed a double-peaked radio light curve accompanied by a detection in the X-rays, but no evidence for broad H$\alpha$ spectral features. The Very Long Baseline Interferometry (VLBI) observations presented here show that the PTF11qcj GHz radio ejecta remains marginally resolved at the sub-milliarcsecond level at $\approx 7.5$ yrs after the explosion, pointing toward a non-relativistic expansion. Our VLBI observations thus favor a scenario in which the second peak of the PTF11qcj radio light curve is related to strong interaction of the supernova ejecta with a circumstellar medium of variable density, rather than to the emergence of an off-axis jet. Continued VLBI monitoring of PTF11qcj in the radio may strengthen further this conclusion.


INTRODUCTION
Supernovae (SNe) of type Ib/c are believed to mark the deaths of massive stars that are stripped of their hydrogen (type Ib), and possibly helium (type Ic), envelope before explosion (Filippenko 1997). A sub-class of Ib/c SNe dubbed broad-line (BL) Ic, estimated to constitute only ≈ 5% of the Ib/c population (Woosley & Bloom 2006;Gal-Yam 2017), is of particular interest due to its relation to long-duration gamma-ray bursts (GRBs), the most relativistic stellar explosions we know of in the universe (Piran 2004;Mészáros 2006). While all GRB-associated SNe are of type BL-Ic (Woosley & Bloom 2006;Hjorth & Bloom 2012, but see Cano et al. (2014) for the peculiar case of SN 2013ez), not all BL-Ic events make a GRB (e.g., Berger et al. 2003;Soderberg et al. 2006;Corsi et al. 2016). Thus, the question of what physical ingredients enable some stripped-envelope massive stars to launch a relativistic jet remains open (e.g., Modjaz et al. 2016).
As first demonstrated by the well-known case of SN 1998bw/GRB 980425 (Galama et al. 1998;Kulkarni et al. 1998), radio observations are particularly well suited to identify those BL-Ic SNe that may harbor GRBs (hereafter referred to as engine-drive SNe), since radio synchrotron emission traces the fastest moving ejecta (e.g., Berger et al. 2003). At the same time, because non-thermal radio photons are produced in the interaction of the SN shock with the circumstellar material (CSM), bright radio emission can also be the smok-ing gun for (non-relativistic) ejecta interacting with a high-density CSM (e.g., Chevalier 1998;Chevalier et al. 2004;Chevalier & Fransson 2006). Although strong CSM interaction is not commonly observed in BL-Ic SNe, a few cases exists such as SN 2007bg (Salas et al. 2013), PTF11qcj (Corsi et al. 2014;Palliyaguru et al. 2019), SN2018gep (Ho et al. 2019a) and, possibly, AT2018cow (Rivera Sandoval et al. 2018;Smartt et al. 2018;Ho et al. 2019b;Margutti et al. 2019).
Here, we focus on PTF11qcj, a radio-loud BL-Ic SN extensively monitored via our approved programs on the Karl G. Jansky Very Large Array (VLA; Corsi et al. 2014;Palliyaguru et al. 2019). The extraordinary radio luminosity of PTF11qcj (∼ 10 29 erg s −1 Hz −1 ) is reminiscent of the GRB-associated SN 1998bw (Kulkarni et al. 1998). As discussed in Corsi et al. (2014) andPalliyaguru et al. (2019), our radio monitoring over the first ≈ 5.5 years since explosion has revealed an unusual double-peaked radio light curve. The radio emission observed during the first light curve peak (t 200 d) can be modeled within the standard synchrotron self-absorption (SSA) model for a spherical SN shock expanding in the CSM. This model yields an estimated speed of ≈ 0.3 − 0.5 c for the fastest SN ejecta (Corsi et al. 2014), placing PTF11qcj in an intermediate class between "ordinary" BL-Ic SNe and engine-driven ones like SN 1998bw or SN 2009bb (Soderberg et al. 2010. The simple, spherically symmetric model of SN shock interaction with a smooth CSM (simple power-law density profile), however, cannot explain the second radio peak. As discussed in Palliyaguru et al. (2019), two more complex scenarios can be invoked to interpret this peculiar behavior of PTF11qcj: (i) A spherical SN shock going through a medium with extreme CSM density variations, perhaps related to eruptive progenitor mass loss; (ii) A radio-emitting SN shock (first peak) followed by radio emission from an emerging off-axis GRB jet, initially pointed away from our line of sight (second peak).
In Palliyaguru et al. (2019) we have shown that, while modeling of our VLA dataset within scenario (i) can in-deed explain the second radio peak, the presence of an off-axis jet (scenario (ii)) cannot be ruled out just based on light curve measurements. However, scenarios (i) and (ii) make rather different predictions for the angular size of the PTF11qcj ejecta at very late times. Motivated by these considerations, here we present Very Long Baseline Interferometry (VLBI) observations of PTF11qcj aimed at setting direct constraints on the size (angular diameter) of its radio ejecta. These observations ultimately provide a direct test for the presence of relativistic expansion, as expected in the case of an off-axis GRB jet.
Our paper is organized as follows. In Section 2, we present the HSA and eMERLIN observations of PTF11qcj. In Section 3, we discuss these observations within the light curve and radio ejecta size predictions of the two scenarios mentioned above. Finally, in Section 4 we summarize our results and conclude.

HSA observations
We observed the field of PTF11qcj at 1.66 GHz (project code BP229A, PI: Palliyaguru) and 15.37 GHz (BP229B, PI: Palliyaguru) with the High Sensitivity Array (HSA) on 2018 December 08.37 UTC and 2019 April 28.97 UTC. The HSA included, in both bands, the Very Long Baseline Array (VLBA, USA) and the Effelsberg 100 m antenna (Germany). At 1.66 GHz, we also used the VLA for improved sensitivity. Both observations covered 128 MHz continuum bandwidth and were correlated at the National Radio Astronomy Observatory (NRAO) Array Operations Center (AOC) in Socorro (New Mexico, USA) with averaging times of 2 s and 1 s at 1.66 GHz and 15.37 GHz, respectively. In both observations, we correlated the target data at R.A. 13 h 13 m 41.5100 s , Dec. +47 • 17 57.600 (J2000; Corsi et al. 2014). In setting up our observations, we used standard phase referencing, where scans on target are interleaved with scans on a nearby compact complex gain calibrator with known position. At 1.66 GHz and 15.37 GHz, we used J1310+4653 and J1358+4737, respectively, as our complex gain calibrators. The VLBA observations were correlated using the NRAO's implementation of the DiFX software correlator (Deller et al. 2011).
We performed all calibration and imaging procedures using the 31DEC19 release of the Astronomical Image Processing System (AIPS; Greisen 2003) and Parsel-Tongue (Kettenis et al. 2006). PTF11qcj was clearly detected at both frequencies. At 1.66 GHz, in particular, given the relatively large separation between PTF11qcj and the complex gain calibrator J1358+4737 (≈ 7.6 • ) we used self-calibration to correct for the residual phase errors towards PTF11qcj and obtain a reliable flux density measurement. The final VLBI images are presented in Figure 1.
We fitted Gaussian intensity profiles to the target images to obtain the flux density and position of PTF11qcj at each frequency. Our results are summarized in Table 1. The reported flux density uncertainties are the quadrature sum of the instrumental calibration uncertainty plus the uncertainty of the Gaussian fit. At 15.37 GHz, we adopt an instrumental uncertainty of 15%. At 1.66 GHz, we use an uncertainty of 20% to also account for residual errors due to the large separation be- tween the target and the complex gain calibrator.
To estimate the source size we fit Gaussian intensity distributions to the images. We find that the source appears marginally resolved along the major axis of our 15.37 GHz observations, while the minor axis is consistent with an unresolved source. The deconvolved major (minor) FWHM size of the fitted Gaussian model is 300 ± 71 µarcsec (76 ± 76 µarcsec), with position angle 125 ± 19 deg (see also Table 1). This corresponds to a diameter of (5.2 ± 1.2) × 10 17 cm ((1.3 ± 1.3) × 10 17 cm). Given the relatively weak radio emission from PTF11qcj, although this results gives us an estimate of the source size, we refrain from speculation about the non-symmetrical nature of this result. PTF11qcj radio ejecta is not resolved at 1.66 GHz.

eMERLIN observations
We observed the field of PTF11qcj at 5.07 GHz and 1.51 GHz with all six eMERLIN antennas on 2019 August 01.55 UTC and 2019 August 29.47 UTC via our DDT project DD8011 (PI: Perez-Torres). We used J1310+4653 as our complex gain calibrator, 1 s integration time, and 512 MHz continuum bandwidth in both bands. We used the standard eMERLIN calibrators  3C286 and OQ208 for flux density calibration and bandpass calibration, respectively. We calibrated and edited the correlated data using the eMERLIN CASA pipeline version 1.1.11 (Moldon 2018). We applied self-calibration in both bands to correct for significant residual phase errors and minor amplitude errors. We used WSClean (Offringa et al. 2014) to deconvolve the calibrated data and produce the final images. PTF11qcj is clearly detected in both bands.
We report in Table 2 the peak and total flux densities, along with the position of PTF11qcj, at each frequency, calculated by fitting Gaussian intensity profiles to the images. The quoted flux density uncertainties correspond to the sum in quadrature of the systematic, i.e., instrumental uncertainty (15%) and the image off-source RMS noise.

MODELING
The complete radio light curves of PTF11qcj are shown in Figure 2. These include the latest flux measurements (integrated fluxes from Tables 1 and 2) along with data published in Corsi et al. (2016) and Palliyaguru et al. (2019). Hereafter, we discuss the possible interpretation of these light curves within the the standard synchrotron self-absorption (SSA) scenario for radio SNe (see Soderberg et al. 2005, and references therein). We also consider an alternative interpretation within an SSA (first radio light curve peak) plus off-axis GRB (second radio light curve peak) scenario. Finally, we discuss the direct size constraints obtained via our VLBI observations in the context of both these scenarios.

Light curves within the SSA scenario
Similarly to what done in Corsi et al. (2016) and Palliyaguru et al. (2019), we can model the radio emission of PTF11qcj in the SSA scenario (Soderberg et al. 2005). As described in Soderberg et al. (2005), the temporal evolution of the shock radius, r, minimum Lorentz factor, γ m , and magnetic field, B, are parameterized as: with α r , α B and α γ the temporal indices of the three quantities respectively, t e the explosion time, and t 0 an arbitrary reference epoch since explosion, here set to 10 d. Within the standard assumptions, one has (see Equations (9)-(10) in Soderberg et al. 2005): Here, s describes the density profile of the shocked CSM where the density of the radiating electrons within the shocked CSM is given by (Chevalier 1982) n e = n e,0 t − t e t 0 αn e,0 ∝ r −s .
The above relations follow from assuming that the energy density of shocked particles (protons and electrons) and  Table 3; solid curves). The re-brightening phase is modeled in two different scenarios: (i) within the standard SSA model (dotted curves; Model 3 in Table 3) and (ii) within an off-axis afterglow model (long-dashed, blue curves; see Section 3.2). The sum of the best fit CSM model in the first peak and off-axis jet emission from the second peak is also shown (dot-dashed curve). The dash-dot-dot-dotted line in the bottom right panel shows the effects of synchrotron cooling within the SSA scenario as discussed in Section 3.1. The measurement of the shock size from HSA data rules out the off-axis jet scenario.
is the shock speed. The density of electrons in the shocked CSM is related to the progenitor mass-loss rate via the relation (see Equation (13) in Soderberg et al. 2005): where we have assumed a nucleon-to-proton ratio of 2, v w ∼ 1000 km/s is the velocity of the stellar wind, and where η (typically in the range η ∼ 2 − 10) characterizes the thickness of the radiating electron shell as r/η. In the GHz radio band, the observing frequencies ν are typically such that ν m << ν, where (see Equation (A6) in Soderberg et al. 2005): is the characteristic synchrotron frequency of electrons with Lorentz factor γ m . In the above Equation, m e is the electron mass, and c is the speed of light. In this frequency range, assuming the synchrotron cooling frequency is higher than the observing frequency, and neglecting synchrotron cooling effects, the SSA emission from the shocked electrons reads: (1−e −τ ξ ν (t) ) 1/ξ ν 5/2 , (10) where F is a normalization constant that depends on the parameters (r 0 , B 0 , p) with p the power-law index of the electron energy distribution (see Equations (A11) 7 and (A13) in Soderberg et al. 2005), and where the optical depth τ is given by (see Equations (20) and (A14) in Chevalier & Fransson 2003;Soderberg et al. 2005, respectively): with T a normalization constant that depends on the parameters (r 0 , B 0 , p, γ m,0 , η). We note that for ν m << ν, the value of γ m,0 is left largely unconstrained by the observations, and thus typically fixed so that ν m,0 ∼ 1 GHz. In addition, the thickness of the shell is typically set to a value η > 1 (Li & Chevalier 1999). Thus, for a given choice of η and ν m , the SSA model is a function of the parameters (r 0 , ξ, α r , t e , s, B 0 , p).
Within the SSA scenario, and with the data collected here, we can further test the hypothesis first presented in Palliyaguru et al. (2019) that the double-peaked radio light curve of PTF11qcj is due to strong interaction with 7 The functions defined in (A11) are not time-dependent and therefore included in the normalization constant CSM of variable density. Our results are shown in Figure 2, and are reported in Table 3. Model 0 in Table 3 is the best fit SSA model for the first peak (t 215 days since explosion) as reported in Corsi et al. (2016). The last is obtained by setting t e = 55842, p = 3 (as typically expected for Type Ib/c SNe; see Chevalier & Fransson 2006), η = 10, ν m,0 = 1 GHz, and varying r 0 , ξ, α r , s, and B 0 . Model 0 is also plotted in Figure 2 (solid line). Model 1 is a fit to the second radio peak (t 587 days since explosion) where we keep r 0 and α r fixed to their best fit values for the first peak so as to ensure a smooth radial evolution between the first and second radio light curve peaks, and allow ξ, s, and B 0 to vary. This fit is similar to the one reported in Palliyaguru et al. (2019) but updated to include the eMERLIN and HSA data presented here.
Compared to Palliyaguru et al. (2019), the reduced χ 2 for Model 1 is substantially higher, indicating a worsening of the goodness of fit. Model 2 is a fit with a model identical to Model 1 but where the 15 GHz HSA data has been excluded. The improved χ 2 value for Model 2 compared to Model 1 indicates a significant discrepancy between data and model at the highest radio frequencies, suggesting a steepening in the highest frequency light curve which may be caused by the passage of the cooling frequency in band. Indeed, the effects of synchrotron cooling may become important at the late timescales considered here. Within the SSA scenario, the synchrotron cooling frequency can be written as (see Equation (A16) in Soderberg et al. 2005): where σ T is the Thomson cross-section and e is the electron charge (Rybicki & Lightman 1986;Soderberg et al. 2005). For ν > ν c , the flux density becomes (Soderberg et al. 2005): Model 1 in Table 3 predicts ν c ≈ 50 GHz at t − t e ≈ 7.5 yrs, which is above the highest frequency of our observations, resulting in large residuals with the HSA data point. Thus, in Table 3 we also report the results of a fit where we set t e = 55842, p = 3, s = 0, η = 2, ν m,0 = 3 GHz, r 0 = 1.1 × 10 16 cm, but allow α r , B 0 , and ξ to vary (Model 3). As shown in Figure 3, with this choice the cooling frequency falls below 16 GHz at t − t e ≈ 1345 d (≈ 3.6 yrs since explosion), thus causing a steeping of the light curve at this frequency. In Figure 2 we plot Model 3 with dotted lines, and the expected steeping (f ν (t) ∝ t −2.6 according to Equation (13)) of the 16 GHz light curve due to synchrotron cooling with a dash-dot-dot-dotted line. Model 3 substantially improves the goodness of fit for the second radio peak compared to Models 1 and 2. In summary, while the SSA fits described in this Section have large limitations due to the simplifications that characterize the SSA model, overall they are indicative of the fact that (i) a non-constant CSM profile is needed to explain the PTF11qcj radio light curves, and (ii) synchrotron cooling may be playing a role at late times. We also note that in Model 3 a flat (s = 0) radial profile of the environment (ISM-like) is favored (see Chevalier  1982, for a discussion of the s = 0 case). Finally, as shown in Figure 4, the combined radial evolution of the shock as implied by Model 0 for the first peak of the PTF11qcj light curves and by Model 3 for the second peak, implies that the shock has decelerated while encountering the CSM discontinuity. A similar case is of SN2014C, where VLBI observations revealed that the SN has substantially decelerated at late times after encountering a higher density shell (Bietenholz et al. 2020).

SSA plus off-axis GRB scenario
Similarly to what done in Palliyaguru et al. (2019), we also use numerical simulations for off-axis GRB jets by van Eerten et al. (2012) to model the second peak of PTF11qcj light curve within a constant density envi-ronment, considering the new HSA and eMERLIN data points. Within this scenario, the first radio peak is the radio emission from the SN ejecta while the second radio peak is due to the delayed emission from an initially off-axis jets that enters our line of sight after spreading and decelerating. These two-dimensional hydrodynamic simulations of the GRB jets take into account the Blandford-Mckee solution (Blandford & McKee 1976) in the relativistic regime, the Sedov-von Neumann-Taylor (SNT) solution (Taylor 1950) in the late non-relativistic regime, and a transition in between regimes (Zhang & MacFadyen 2009).
In Table 4 we report the best fit results for the isotropic equivalent kinetic energy of the explosion E iso , the jet half-opening angle, θ 0 , ISM density n ISM , and the observer angle, θ obs , when the latest four data points from HSA and eMERLIN, as well as the Chandra X-ray data point reported in Palliyaguru et al. (2019), are added to the fit. These fits also assume equipartition of energy between particles and magnetic fields such that e = B = 0.33. We set the electron energy index, p = 2.5, which is typical for GRB afterglows, as expected from theoretical considerations (Kirk et al. 2000;Achterberg et al. 2001). The best fit off-axis jet light curves are shown in Figure 2 (light blue, long-dashed lines).
3.3. Size constraints: CSM-interacting vs off-axis GRB scenario Within the SSA scenario, Models 0, 1, and 2 described in Section 3.1 all imply a shock radius around the time of the 16 GHz HSA observation (2759 days post explosion) of r ∼ 10 18 cm. This in turn corresponds to an angular diameter of ∼ 1 mas at the redshift z = 0.028 of PTF11qcj, larger than the size constraints set by our HSA observations (see Table 1). On the other hand, Model 3 gives r ≈ 3 × 10 17 cm at 2759 days post explosion (see Figure 4), which corresponds to an angular diameter of ≈ 0.4 mas, compatible with the HSA observations reported in Table 1. We note however that within the SSA model, which assumes spherical symmetry, we are not able to model any potential asymmetry in the shock geometry.
Within the off-axis GRB scenario for the second radio peak, the ejecta will no longer produce a symmetric image on the sky. Instead, the resolved image will highlight the front edge of the jet heading to the observer, producing an elongated and curved shape. For model 1 in Table 4, this curved front will have traveled R ⊥,travel = 3.0 × 10 19 cm in projected distance from the origin of the explosion. The projected width of the image is R ⊥,w = 1.4 × 10 19 cm, and its projected height R ⊥,h = 2.5 × 10 18 cm. These correspond to a projected angular distance from the origin of the explosion of 35 mas, a projected angular width of 8 mas, and a projected angular height of 14 mas. Thus, if the second radio light curve peak of PTF11qcj was due to an off-axis GRB, we should have resolved a much larger image in our HSA observations.
We note that while the off-axis GRB model considered here assumes a top-hat jet, a structured jet such as the one considered in Nakar & Piran (2017) would imply sizes that can be bracketed by the the two extreme cases of a spherically symmetric blast wave, and a non-spreading relativistic cone. Indeed, for a spherically symmetric Sedov-Taylor blast wave, the analytical expression for the shock radius at a late-time t is given by (van Eerten 2018): where C R (k) = 1.5 parsec for the power law index of the CSM density profile k = 1.4 (Model 1, Table 3), E j is the total energy in the ejecta, ρ ref and R ref are the circumburst medium density and radius at a reference time (which we set to 10 days), respectively (van Eerten 2018). For the best fit parameters of Model 1 in Table 3, we find the expected radius would be 8.9×10 17 cm (1 mas), larger than our HSA constraint.
For a relativistic jet expanding into the ISM, the apparent radius at time t may be calculated using the analytical expression (Oren et al. 2004): R ⊥,exp = 5 × 10 16 E 51 n 1/6 T −1/8 j t 5/8 , where T j = (E 51 /n) 1/3 (θ 0 /0.1) 2 (1 + z) is the jet break time, θ 0 is the jet half opening angle, n is the ISM number density (Oren et al. 2004). For the best fit parameters of E iso , θ 0 , n ISM listed in Table 4, the expected radius at 2759 days post explosion is 4.8 × 10 19 cm. This corresponds to an angular size of ∼ 55 mas at the redshift z = 0.028 of PTF11qcj, much larger than our HSA constraint. Based on these results, the off-axis hypothesis for the origin of the second radio peak of the PTF11qcj light curves is disfavored.

SUMMARY AND CONCLUSION
We have presented HSA and eMERLIN observations of PTF11qcj obtained ∼7.5 yrs post explosion. The source is marginally resolved in the HSA data at 15 GHz, indicating a diameter of (300 ± 71) µas. This corresponds to an average expansion velocity of ≈ 0.036 ± 0.008 c.
Modeling of the light curves within an SSA scenario requires the interaction of the shock with a non-smooth CSM whose radial profile changes from a stellar wind profile (n e ∝ r −2 ) to a constant density medium. This variable CSM profile also affects the temporal evolution of the shock radius. We also find that synchrotron cooling may be playing a role at the highest radio frequencies. Within this SSA model with variable CSM, the derived size of the shock at the time of our HSA observations can be reconciled with our measurements. Our results are thus consistent with the CSM interaction model (considering that such model is highly simplified). The size constraints derived from our HSA observations also seem to disfavor an off-axis relativistic jet scenario for the radio re-brightening of PTF11qcj.
In conclusion, the VLBI observations reported in this paper favor the original interpretation of PTF11qcj as a strongly CSM-interacting radio SN (see Corsi et al. 2014;Palliyaguru et al. 2019). However, our conclusions are limited by the simplifications inherent in the spherically symmetric SSA model adopted here. We encourage more detailed theoretical modeling aimed at interpreting PTF11qcj complex radio light curve and available VLBI data.