VLBI Observations of Supernova PTF11qcj: Direct Constraints on the Size of the Radio Ejecta

We observed the field of PTF11qcj at 1.66 GHz (project code BP229A, PI: Palliyaguru) and 15.37 GHz (BP229B, PI: Palliyaguru) with the 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 the 128 MHz continuum bandwidth and were correlated at the National Radio Astronomy Observatory (NRAO) Array Operations Center 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^h13^m415100, decl. $+47^\circ 17^{\prime} 57\buildrel{\prime\prime}\over{.} 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 (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 (≈76), we used self-calibration to correct for the residual phase errors toward PTF11qcj and obtain a reliable flux density measurement. The final VLBI images are presented in Figure 1.

Zoom In Zoom Out Reset image size

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 between the target and the complex gain calibrator.

Table 1. HSA Results for PTF11qcj (Project Codes BP229A and BP229B)

	1.66 GHz	15.37 GHz
Observing date (UT)	2018-12-08	2019-04-28
Project code	BP229A	BP229B
Observing time including calibrators (hr)	8	8
Image off-source rms noise (μJy beam−1)	125	24
CLEAN restoring beam (mas)2	4.06 × 2.66	0.54 × 0.24
CLEAN beam position angle (deg)	−29	−21
Peak flux density (mJy)	4.32 ± 0.92	0.494 ± 0.054
Integrated flux density (mJy)	5.8 ± 1.4	0.681 ± 0.085
R.A. [J2000]	13h13m4147512	13h13m4147490
Decl. [J2000]	+47°17'568017	+47°17'567988
Deconvolved fitted major axis	⋯	300 ± 71 μarcsec
Deconvolved fitted minor axis	⋯	76 ± 76 μarcsec
Deconvolved fitted Pos. ang.	⋯	125° ± 19°

Note. See Section 2.1 for discussion.

Download table as:  ASCIITypeset image

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 μas (76 ± 76 μas), with a position angle of 125° ± 19° (see also Table 1). This corresponds to a diameter of (5.2 ± 1.2) × 10¹⁷ cm ((1.3 ± 1.3) × 10¹⁷ cm). Given the relatively weak radio emission from PTF11qcj, although this result gives us an estimate of the source size, we refrain from speculating about the nonsymmetrical nature of this result. The radio ejecta of PTF11qcj is not resolved at 1.66 GHz. We note an apparent positive extension in Figure 1(a), but also a significant negative peak nearby. Such an asymmetry is very likely due to residual phase errors, even after self-calibration, caused by the significant separation between target and complex gain calibrator in BP229A. We, therefore, refrain from further interpretation of this 5σ level structure.

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 3C 286 and OQ 208 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.

Table 2.  eMERLIN Results for PTF11qcj (Project Code DD8011)

	1.5 GHz	5.07 GHz
Observing date (UT)	2019-08-29	2019-08-01
Project code	DD8011	DD8011
Observing time including calibrators (hr)	6	5
Image off-source rms noise (μJy beam−1)	42	37
CLEAN beam (mas2)	59 × 19	622 × 78
CLEAN beam position angle (deg)	40	42
Peak flux density (mJy)	5.82 ± 0.87	4.15 ± 0.62
Integrated flux density (mJy)	5.78 ± 0.87	4.45 ± 0.67
R.A. [J2000]	13h13m414745	13h13m414746
Decl. [J2000]	+47°17'56795	+47°17'56800

Note. See Section 2.2 for discussion.

Download table as:  ASCIITypeset image

Similarly to what was 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

$$\begin{eqnarray}&&r={r}_{0}{\left(\displaystyle \frac{t-{t}_{e}}{{t}_{0}}\right)}^{{\alpha }_{r}},\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&B={B}_{0}{\left(\displaystyle \frac{t-{t}_{e}}{{t}_{0}}\right)}^{{\alpha }_{B}},\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&{\gamma }_{m}={\gamma }_{m,0}{\left(\displaystyle \frac{t-{t}_{e}}{{t}_{0}}\right)}^{{\alpha }_{\gamma }},\end{eqnarray} \tag{ 3 }$$

with α_r , α_B , and α_γ the temporal indices of the three quantities, respectively; t_e the explosion time, and t₀ an arbitrary reference epoch since explosion, here set to 10 days. Within the standard assumptions, one has (see Equations (9)–(10) in Soderberg et al. 2005)

$$\begin{eqnarray}&&{\alpha }_{\gamma }=2({\alpha }_{r}-1),\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&{\alpha }_{B}=\displaystyle \frac{(2-s)}{2}{\alpha }_{r}-1.\end{eqnarray} \tag{ 5 }$$

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)

$$\begin{eqnarray}&&{n}_{e}={n}_{e,0}{\left(\displaystyle \frac{t-{t}_{e}}{{t}_{0}}\right)}^{{\alpha }_{{n}_{e,0}}}\propto {r}^{-s}.\end{eqnarray} \tag{ 6 }$$

The above relations follow from assuming that the energy density of shocked particles (protons and electrons) and amplified magnetic fields are a constant fraction (_e ≈ _B ≈ 0.33 under the hypothesis of equipartition) of the postshock energy density U ∝ n_e v², where

$$\begin{eqnarray}&&v={v}_{0}{\left(\displaystyle \frac{t-{t}_{e}}{{t}_{0}}\right)}^{{\alpha }_{r}-1}\end{eqnarray} \tag{ 7 }$$

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)

$$\begin{eqnarray}&&\dot{M}=\displaystyle \frac{8\pi }{\eta }{n}_{e,0}{m}_{p}{r}_{0}^{2}{v}_{w}{\left(\displaystyle \frac{t-{t}_{e}}{{t}_{0}}\right)}^{{\alpha }_{r}(2-s)}={\dot{M}}_{0}{\left(\displaystyle \frac{t-{t}_{e}}{{t}_{0}}\right)}^{{\alpha }_{\dot{M}}},\end{eqnarray} \tag{ 8 }$$

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)

$$\begin{eqnarray}&&{\nu }_{m}={\gamma }_{m}^{2}\left(\displaystyle \frac{{eB}}{2\pi {m}_{e}c}\right)\end{eqnarray} \tag{ 9 }$$

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 ignoring synchrotron cooling effects, the SSA emission from the shocked electrons reads

$$\begin{eqnarray}&&{f}_{\nu }{\left(t)={ \mathcal F }{\left(\displaystyle \frac{t-{t}_{e}}{{t}_{0}}\right)}^{(4{\alpha }_{r}-{\alpha }_{B})/2}(1-{e}^{-{\tau }_{\nu }^{\xi }(t)}\right)}^{1/\xi }{\nu }^{5/2},\end{eqnarray} \tag{ 10 }$$

where ${ \mathcal F }$ is a normalization constant that depends on the parameters (r₀, B₀, p) with p the power-law index of the electron energy distribution (see Equations (A11) ⁶ 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)

$$\begin{eqnarray}&&\tau (t)={ \mathcal T }\,{\left(\displaystyle \frac{t-{t}_{e}}{{t}_{0}}\right)}^{(p-2){\alpha }_{\gamma }+(3+p/2){\alpha }_{B}+{\alpha }_{r}}\,{\nu }^{-(p+4)/2},\end{eqnarray} \tag{ 11 }$$

with ${ \mathcal T }$ a normalization constant that depends on the parameters (r₀, B₀, 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₀, ξ, α_r , t_e , s, B₀, 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 the strong interaction with 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, and ν_m,0 = 1 GHz, and varying r₀, ξ, α_r , s, and B₀. 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₀ 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₀ 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.

Table 3. Best-fit Parameters for the Standard SSA Model Described in Section 3.1

Parameter	Model 0	Model 1	Model 2	Model 3
r0 (cm)	1.1 × 1016	1.1 × 1016 (fixed)	1.1 × 1016 (fixed)	1.1 × 1016 (fixed)
ξ	0.23	0.19	0.19	0.34
αr	0.79	0.79 (fixed)	0.79 (fixed)	0.61
te	55842 (fixed)	55842 (fixed)	55842 (fixed)	55842 (fixed)
s	2.0 (fixed)	1.4	1.1	0.0 (fixed)
B0 (G)	6.7	5.8	3.2	1.35
p	3.0 (fixed)	3.0 (fixed)	3.0 (fixed)	3.0 (fixed)
αB	−1.0	−0.76	−0.64	−0.39
γm,0	7.3	7.8	10.5	28.1
αγ	−0.4	−0.4	−0.4	−0.8
ne,0 (cm−3)	1.5 × 105	1.0 × 105	2.4 × 104	1.6 × 103
${\alpha }_{{n}_{e}}$	−1.6	−1.1	1-0.88	0.0
${\dot{M}}_{0}\,({M}_{\odot }\,{\mathrm{yr}}^{-1})$	1.2 × 10−4	8.4 × 10−5	1.9 × 10−5	6.4 × 10−6
${\alpha }_{\dot{M}}$	0.0	0.48	0.71	1.2
χ2/dof	1793/90	1288/39	602/38	211/38

Download table as:  ASCIITypeset image

Compared to Palliyaguru et al. (2019), the reduced χ² 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 have been excluded. The improved χ² 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 the 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)

$$\begin{eqnarray}&&{\nu }_{c}=\displaystyle \frac{18\pi \,{m}_{e}\,c\,e}{{(t-{t}_{e})}^{2}{\sigma }_{{\rm{T}}}^{2}\,{B}^{3}},\end{eqnarray} \tag{ 12 }$$

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)

$$\begin{eqnarray}&&{f}_{\nu }(t)\propto {\nu }^{-p/2}\,{t}^{(6{\alpha }_{r}+(8-5p){\alpha }_{B}+2(p-2){\alpha }_{\gamma }-4p+2)/2}.\end{eqnarray} \tag{ 13 }$$

Model 1 in Table 3 predicts ν_c ≈ 50 GHz at t − t_e ≈ 7.5 yr, 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, and r₀ = 1.1 × 10¹⁶ cm, but allow α_r , B₀, 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 days (≈3.6 yr 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 dashed–triple-dotted line. Model 3 substantially improves the goodness of fit for the second radio peak compared to Models 1 and 2.

Zoom In Zoom Out Reset image size

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 nonconstant 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 (interstellar medium (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 SN 2014C, where VLBI observations revealed that the SN has substantially decelerated at late times after encountering a higher-density shell (Bietenholz et al. 2021).

Zoom In Zoom Out Reset image size

Similarly to what was 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 the PTF11qcj light curve within a constant-density environment, 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 jet 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 solution (Taylor 1950) in the late nonrelativistic regime, and a transition 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 θ₀, 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 the 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).

Table 4. Best-fit Results in the Off-axis GRB Scenario

Model	Eiso	n	θ0	θobs	p	χ2	R⊥	Angular Size	${R}_{\perp ,\exp }$	Expected Angular Size
	(erg)	(cm−3)	(rad)	(rad)			(cm)	(mas)	(cm)	(mas)
Model 1	1 × 1053	1 × 10−5	0.2	0.4	2.5	1285/41 ≈ 31	3.0 × 1019	35	4.8 × 1019	55

Note. The radius (R_⊥) and the corresponding angular diameter from modeling, and the expected radius from Equation (15) (${R}_{\perp ,\exp }$), and the corresponding angular diameter are also listed.

Download table as:  ASCIITypeset image

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 postexplosion) of r ∼ 10¹⁸ 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¹⁷ cm at 2759 days postexplosion (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¹⁹ cm in projected distance from the origin of the explosion. The projected width of the image is R_⊥,w = 1.4 × 10¹⁹ cm, and its projected height R_⊥,h = 2.5 × 10¹⁸ 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 two extreme cases of a spherically symmetric blast wave, and a nonspreading 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)

$$\begin{eqnarray}\begin{array}{rcl}R & = & {C}_{R}(k){\left(\displaystyle \frac{{E}_{j}}{{10}^{51}\mathrm{erg}}\right)}^{\tfrac{1}{5-k}}\,{\left(\displaystyle \frac{{\rho }_{\mathrm{ref}}}{{m}_{p}}\right)}^{\tfrac{-1}{5-k}}\\ & & \times {\left(\displaystyle \frac{{R}_{\mathrm{ref}}}{c\times {10}^{6}}\right)}^{\tfrac{-k}{5-k}}\,{\left(\displaystyle \frac{t}{{10}^{8}}\right)}^{\tfrac{2}{5-k}},\end{array}\end{eqnarray} \tag{ 14 }$$

where C_R (k) = 1.5 pc 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, and ρ_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¹⁷ 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)

$$\begin{eqnarray}&&{R}_{\perp ,\exp }=5\times {10}^{16}\,{\left(\displaystyle \frac{{E}_{51}}{n}\right)}^{1/6}{T}_{j}^{-1/8}{t}^{5/8},\end{eqnarray} \tag{ 15 }$$

where ${T}_{j}={({E}_{51}/n)}^{1/3}{({\theta }_{0}/0.1)}^{2}(1+z)$ is the jet break time, θ₀ is the jet half-opening angle, and n is the ISM number density (Oren et al. 2004). For the best-fit parameters of E_iso, θ₀, and n_ISM listed in Table 4, the expected radius at 2759 days postexplosion is 4.8 × 10¹⁹ 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.