# Radio evolution of a Type IIb supernova SN 2016gkg

Nayana A. J. Indian Institute of Astrophysics, II Block, Koramangala,
Bangalore 560034, India.  Poonam Chandra National Radio Astronomy Observatory,
520 Edgemont Road, Charlottesville VA 22903, USA  National Centre for Radio
Astrophysics, Tata Institute of Fundamental Research, PO Box 3, Pune, 411007,
India  Anoop Krishna  N.S.S College, Ottapalam, Palakkad 679103, India  G.C.
Anupama  Indian Institute of Astrophysics, II Block, Koramangala, Bangalore
560034, India.

(Received April 8, 2022; Revised June 16, 2022)

###### Abstract

We present extensive radio monitoring of a type IIb supernova (SN IIb), SN
2016gkg during $t\sim$ 8$-$1429 days post explosion at frequencies $\nu\sim$
0.33$-$25 GHz. The detailed radio light curves and spectra are broadly
consistent with self-absorbed synchrotron emission due to the interaction of
the SN shock with the circumstellar medium. The model underpredicts the flux
densities at $t\sim$ 299 days post-explosion by a factor of 2, possibly
indicating a density enhancement in the CSM due to a non-uniform mass-loss
from the progenitor. Assuming a wind velocity $v_{\rm w}\sim$ 200 km s-1, we
estimate the mass-loss rate to be $\dot{M}\sim$ (2.2, 3.6, 3.8, 12.6, 3.7, and
5.0) $\times$ 10-6 $M_{\odot}$ yr-1 during $\sim$ 8, 15, 25, 48, 87, and 115
years, respectively before the explosion. The shock wave from SN 2016gkg is
expanding from $R\sim$ 0.5 $\times$ 1016 to 7 $\times$ 1016 cm during $t\sim$
24$-$492 days post-explosion indicating a shock deceleration index, $m$ $\sim$
0.8 ($R\propto t^{m}$), and mean shock velocity $v\sim$ 0.1c. The radio data
being inconsistent with free-free absorption model and higher shock velocities
are in support of a relatively compact progenitor for SN 2016gkg.

Supernovae: general — supernovae: SN 2016gkg — radiation mechanisms: non-
thermal — circumstellar matter — radio continuum: general

††journal: ApJ††facilities: Giant Metrewave Radio Telescope, Karl J. Jansky
Very Large Array††software: CASA (McMullin et al., 2007), AIPS (van Moorsel et
al., 1996), emcee (Foreman-Mackey et al., 2013)

##  1 Introduction

Type IIb supernovae (SNe IIb) are a sub-class of core-collapse supernovae
(CCSNe) characterized by the presence of broad HI absorption features in the
early optical spectra. At later times, these HI lines disappear and He I
feature becomes dominant in the spectra (Filippenko, 1997), placing SNe IIb in
between hydrogen-rich type II SNe and hydrogen poor type Ibc SNe. The
progenitors of SNe IIb are understood to be stars that have lost most of their
hydrogen envelope but not all. The hydrogen envelope could be lost either via
radiatively driven winds (Smith & Conti, 2008) or via mass transfer by a
binary companion (Yoon et al., 2010).

Progenitor candidates have been identified for a few SNe IIb from pre-
explosion images. A K-type supergiant in a binary system for SN 1993J
(Aldering et al., 1994; Maund et al., 2004), a yellow supergiant for SN 2011dh
(Arcavi et al., 2011; Maund et al., 2011; Van Dyk et al., 2013; Sahu et al.,
2013), and a massive star of mass $M\sim$ 20$-$25 $M_{\odot}$ for SN 2008ax
(Crockett et al., 2008; Taubenberger et al., 2011). Besides these direct
detection efforts, the luminosity evolution of early light curve that maps the
cooling envelope phase after the shock breakout can also put constraints on
mass and radius of the progenitor star (Baron et al., 1993; Swartz et al.,
1993).

Independent constraints on the progenitor properties can be obtained by
studying the non-thermal radio emission from SNe IIb that arises as a result
of the dynamical interaction between the supernova (SN) shock and the
circumstellar medium (CSM; Chevalier, 1982, 1998). Radio observations uniquely
probe the density structure of the CSM and thereby a longer period of mass-
loss of the progenitor star (Weiler et al., 1986; Chevalier, 1982, 1998).
Several SNe IIb exhibit luminous radio emission; a few examples are SN 1993J
(Weiler et al., 2007), SN 2001gd (Stockdale et al., 2003, 2007), SN 2001ig
(Ryder et al., 2004), SN 2003bg (Soderberg et al., 2006), SN 2011dh (Krauss et
al., 2012; Soderberg et al., 2012), SN 2011hs (Bufano et al., 2014), and SN
2013df (Kamble et al., 2016).

Chevalier & Soderberg (2010) compiled a sample of radio bright SNe IIb and
divided them based on their radio properties. The authors proposed two
populations of SNe IIb; one with compact progenitors (SNe cIIb) and the other
with extended progenitors (SNe eIIb). The SNe cIIb group shows faster shock
velocities, less dense CSM, and compact progenitors in comparison with that of
SNe eIIb. However, there exist a few examples (e.g., SN 2011dh Bersten et al.,
2012; Horesh et al., 2013; Maeda et al., 2014) that suggest that the radio
properties may or may not be a good indicator of the progenitor size. The
progenitors of SNe IIb could be a continuum of objects between compact and
extended stars instead of a sharp split like SNe eIIb and SNe cIIb.

This paper presents the radio follow-up observations of a type IIb supernova
(SN IIb) SN 2016gkg from $t\sim$ 8 to 1429 days over a frequency range of
0.3$-$25 GHz. The data include our observations taken with the Giant Metrewave
Radio Telescope (GMRT) and the archival data from the Jansky Very Large Array
(JVLA). We model the radio emission to investigate the mass-loss history of
the progenitor system, the evolution of SN shock radius & magnetic field, and
irregularities in the CSM density.

The paper is organized as follows. In §2, we present the compilation of
various results on SN 2016gkg from the literature. The details of observations
and data reduction are presented in §3. We discuss the radio emission model
and derive various parameters of the progenitor and environments in §4 and §5.
The results are discussed in §6 and conclusions are drawn in §7.

##  2 SN 2016gkg

SN 2016gkg was discovered by Buso & Otero on 2016 Sep 20.18
(UT)111http://ooruri.kusastro.kyoto-u.ac.jp/mailarchive/vsnet-alert/20188 at a
position $\alpha_{\rm J2000}$ = 01h34m14.40s, $\delta_{\rm J2000}$ =
$-$29∘26′24.20′′. The SN is located at a distance of 26.4 $\pm$ 5.3 Mpc in the
galaxy NGC 613 (Nasonova et al., 2011). The SN was classified as a SN IIb
based on the optical spectroscopic observations (Tartaglia et al., 2017).
Kilpatrick et al. (2017) model the early time optical light curve of SN
2016gkg and derive the date of explosion to be $t_{\rm 0}=$ 2016 September
20.15${}^{+0.08}_{-0.10}$ UT. We adopt this as the date of explosion through
out this paper and all epochs $t$ mentioned are with respect to $t_{\rm 0}$.

SN 2016gkg was extensively followed in the optical bands soon after the
discovery which provided excellent coverage of its early evolution. The SN
showed double peak structure in the optical light curve, the early peak due to
the shock cooling of the hydrogen envelope of the progenitor star, and the
later peak powered by radioactive decay (Bersten et al., 2018; Kilpatrick et
al., 2017; Tartaglia et al., 2017). Pre-explosion images of the field
containing SN 2016gkg taken in 2001 with the Hubble Space Telescope (HST) Wide
Field Planetary-Camera 2 (WFPC2) is available in the archive (Bersten et al.,
2018).

Various groups (Kilpatrick et al., 2017; Tartaglia et al., 2017; Piro et al.,
2017; Arcavi et al., 2017; Bersten et al., 2018) attempted to constrain the
properties of the progenitor of SN 2016gkg by modeling the early peak of the
light curve by shock cooling models. The progenitor radius estimates from
these studies span a wide range $R_{*}\sim$ 40 $-$ 646 $R_{\odot}$ depending
on the model and assumed structure of the hydrogen envelope of the progenitor.
Piro et al. (2017) investigated the early peak of the optical light curve by
numerically exploding a large number of extended envelope models and constrain
the radius to be $R_{*}\sim$ 180–260 $R_{\odot}$. Arcavi et al. (2017) fit the
observed light curve with analytical shock cooling models (Nakar & Piro, 2014;
Piro, 2015; Sapir & Waxman, 2017) and estimate the progenitor radius to be
$R_{*}\sim$ 40–150 $R_{\odot}$. Bersten et al. (2018) modeled the cooling peak
and estimated the radius of the hydrogen envelope to be $\sim$ 320
$R_{\odot}$. Modelling the initial rapid rise of the light curve using Rabinak
& Waxman (2011) model, Kilpatrick et al. (2017) constrained the progenitor
radius to be $R_{*}=$ 257${}^{+389}_{-189}$. Tartaglia et al. (2017) modeled
the temperature evolution of initial peak of the light curve and estimated the
progenitor radius to be $R_{*}\sim$ 48$-$124 $R_{\odot}$.

Kilpatrick et al. (2017) detected a progenitor candidate of SN 2016gkg in the
archival HST image and estimated the luminosity and radius to be log(L/L⊙) =
5.15 and $R_{*}=$138${}^{+131}_{-103}$ $R_{\odot}$, respectively. The authors
found that single star stellar evolution models fail to reproduce the derived
progenitor properties whereas binary evolutionary tracks could reproduce them.
Tartaglia et al. (2017) identified two plausible progenitor candidates from
HST imaging analysis and suggested a range of progenitor mass 15$-$20
$M_{\odot}$ and radius (150$-$320) $R_{\odot}$. Kilpatrick et al. (2021)
present post-explosion late-time ($t\sim$ 652$-$1795 days) HST observations of
SN 2016gkg and their improved astrometric allignment between the SN and
progenitor candidate allowed them to constrain the progenitor to be a compact
yellow super-giant of radius $\sim$ 70 $R_{\odot}$ with effective temperature
($T_{\rm eff}$) $\sim$ 10800 K. Late time ($t\sim$ 300–800 days) spectroscopic
observations of SN 2016gkg reveal multi component emission lines indicating
the presence of material with different velocities possibly indicating an
asymmetric explosion (Kuncarayakti et al., 2020).

##  3 Observations and Data Reduction

###  3.1 GMRT observations

We carried out regular monitoring of SN 2016gkg with the Giant Metrewave Radio
Telescope (GMRT) during $t\sim$ 51 $-$ 1429 days at 0.33, 0.61, and 1.39 GHz.
The data were recorded with an integration time of 16 seconds in full polar
mode with a bandwidth of 33 MHz, split into 256 channels. 3C 286 and 3C 147
were used as the flux density calibrators and J0240$-$2309 was used as the
phase calibrator.

The GMRT data were inspected and calibrated using the Astronomical Image
Processing Software (AIPS; Greisen, 2003) following standard procedure (see
Nayana et al., 2017). The calibrated visibilities were imaged excluding the
short $uv$ data to minimize the contribution from the extended host galaxy
emission. The host galaxy is of flux density $\sim$ 178.3 mJy at 1.4 GHz with
an angular size $\sim$ 48.5 $\times$ 30.4 arcsec2 in the NVSS map (Condon et
al., 1998). The flux density of the SN was determined by fitting a Gaussian at
the source position using task JMFIT222http://www.aips.nrao.edu/cgi-
bin/ZXHLP2.PL?JMFIT. The SN is $\sim$ 92 arcsec away from the center of the
host galaxy, well beyond the continuum emission from the galaxy. However, we
fit a zero-level baseline while fitting the Gaussian to account for any
residual emission from the host galaxy. The details of GMRT observations are
given in Table 1. The errors on flux densities quoted in Table 1 are the sum
of map rms values and a 10% calibration uncertainty added in quadrature. The
upper limits on radio flux densities are three times the map rms at the
supernova position. The GMRT light curves are shown in Figure 1.

###  3.2 JVLA observations

The Karl G. Jansky Very Large Array (JVLA) observed SN 2016gkg during $t\sim$
8 – 492 days covering a frequency range 2–25 GHz (archival data; PI: Maria
Drout). The observations were done in standard continuum mode with a bandwidth
of 2.048 GHz split into 16 spectral windows each of 128 MHz. 3C 147 was used
as the flux density calibrator, and J0145$-$2733 was used as the phase
calibrator.

The JVLA archival data were reduced using standard packages in Common
Astronomy Software Applications (CASA; McMullin et al., 2007). We use the CASA
task TCLEAN for imaging and we exclude short baselines to minimize the host
galaxy emission while imaging. The flux density of the source is estimated by
fitting a Gaussian at the source position using task IMFIT. The details of
JVLA observations are given in Table 3. The errors on flux densities quoted in
Table 1 are the sum of map rms values and 10% calibration uncertainties added
in quadrature. Fig 1 shows the flux density evolution of SN 2016gkg at
frequencies 2$-$25 GHz.

Table 1: Details of GMRT observations of SN 2016gkg  Date of Observation | AgeaaThe age is calculated assuming 2016 Sep 20.15 (UT) as the date of explosion (Kilpatrick et al., 2017). | Frequency | Flux density
---|---|---|---
(UT) | (Day) | (GHz) | (mJy)bbThe errors on the flux densities are the sum of error from task JMFIT and a 10% calibration uncertainity added in quadrature.
2016 Nov 09.76 | 50.61 | 1.39 | 0.36 $\pm$ 0.08
2016 Dec 10.64 | 81.49 | 1.39 | 0.44 $\pm$ 0.07
2017 Apr 30.16 | 222.01 | 1.39 | 1.65 $\pm$ 0.18
2017 Sep 08.83 | 353.68 | 1.39 | 1.38 $\pm$ 0.15
2017 Nov 23.63 | 429.48 | 1.39 | 1.13 $\pm$ 0.12
2018 Jun 12.98 | 630.83 | 1.39 | 0.78 $\pm$ 0.10
2018 Sep 08.85 | 718.70 | 1.39 | 0.75 $\pm$ 0.09
2019 Jan 27.45 | 859.30 | 1.39 | 0.69 $\pm$ 0.09
2020 Aug 07.09 | 1416.94 | 1.39 | 0.59 $\pm$ 0.10
2016 Dec 12.57 | 83.42 | 0.61 |  $<$0.18
2017 Apr 29.29 | 221.14 | 0.61 | 0.44 $\pm$ 0.09
2017 Sep 08.83 | 353.68 | 0.61 | 1.11 $\pm$ 0.15
2017 Nov 23.63 | 429.48 | 0.61 | 1.67 $\pm$ 0.20
2018 Jun 10.85 | 628.70 | 0.61 | 1.51 $\pm$ 0.34
2018 Sep 09.83 | 719.68 | 0.61 | 1.60 $\pm$ 0.21
2019 Jan 28.41 | 860.26 | 0.61 | 1.43 $\pm$ 0.44
2019 Aug 20.92 | 1064.77 | 0.61 | 1.32 $\pm$ 0.23
2020 Aug 05.05 | 1414.90 | 0.61 | 1.18 $\pm$ 0.32
2020 Aug 18.82 | 1428.67 | 0.325 |  $\leq$3
Table 2: Details of JVLA observations of SN 2016gkg  Date of | AgeaaThe age is calculated using 2016 Sep 20.15 (UT) as the date of explosion (Kilpatrick et al., 2017). | Frequency | VLA | Flux density |
---|---|---|---|---|---
Observation (UT) | (Day) | (GHz) | Array | (mJy)bbThe errors on the flux densities are the sum of error from task IMFIT and a 10% calibration uncertainity added in quadrature. |
2016 Sep 28.39 | 8.24 | 8.549 | A | 0.116 $\pm$ 0.023 |
| 8.24 | 10.999 | A | 0.244 $\pm$ 0.049 |
2016 Oct 14.21 | 24.06 | 2.499 | A | 0.181 $\pm$ 0.031 |
| 24.06 | 3.499 | A | 0.234 $\pm$ 0.097 |
| 24.06 | 4.999 | A | 0.564 $\pm$ 0.063 |
| 24.06 | 7.099 | A | 1.128 $\pm$ 0.115 |
| 24.06 | 8.549 | A | 1.430 $\pm$ 0.151 |
| 24.06 | 10.999 | A | 1.749 $\pm$ 0.184 |
| 24.06 | 19.299 | A | 1.396 $\pm$ 0.184 |
| 24.06 | 24.999 | A | 0.991 $\pm$ 0.139 |
2016 Nov 08.21 | 49.06 | 2.499 | A | 0.509 $\pm$ 0.059 |
| 49.06 | 3.499 | A | 0.960 $\pm$ 0.104 |
| 49.06 | 4.999 | A | 1.708 $\pm$ 0.174 |
| 49.06 | 7.099 | A | 1.976 $\pm$ 0.201 |
| 49.06 | 8.549 | A | 1.939 $\pm$ 0.197 |
| 49.06 | 10.999 | A | 1.489 $\pm$ 0.154 |
| 49.06 | 19.299 | A | 0.581 $\pm$ 0.068 |
| 49.06 | 24.999 | A | 0.336 $\pm$ 0.063 |
2016 Dec 15.01 | 85.86 | 2.499 | A | 1.510 $\pm$ 0.033 |
| 85.86 | 3.499 | A | 2.031 $\pm$ 0.208 |
| 85.86 | 4.999 | A | 1.917 $\pm$ 0.193 |
| 85.86 | 7.099 | A | 1.427 $\pm$ 0.144 |
| 85.86 | 8.549 | A | 1.192 $\pm$ 0.123 |
| 85.86 | 10.999 | A | 0.856 $\pm$ 0.089 |
| 85.86 | 19.299 | A | 0.374 $\pm$ 0.050 |
| 85.86 | 24.999 | A | 0.282 $\pm$ 0.043 |
2017 Mar 19.84 | 299.08 | 1.749 | D |  $<$ 2.4 | -
| 299.08 | 2.499 | D | 1.940 $\pm$ 0.272 |
| 299.08 | 3.499 | D | 1.612 $\pm$ 0.220 |
| 299.08 | 4.999 | D | 1.131 $\pm$ 0.140 |
| 299.08 | 7.099 | D | 0.788 $\pm$ 0.110 |
| 299.08 | 8.549 | D | 0.567 $\pm$ 0.067 |
| 299.08 | 9.499 | D | 0.485 $\pm$ 0.064 |
| 299.08 | 13.499 | D | 0.315 $\pm$ 0.042 |
| 299.08 | 16.499 | D | 0.233 $\pm$ 0.040 |
Table 3: Details of JVLA observations of SN 2016gkg (continued)  Date of | AgeaaThe age is calculated using 2016 Sep 20.15 (UT) as the date of explosion (Kilpatrick et al., 2017). | Frequency | VLA | Flux density
---|---|---|---|---
Observation (UT) | (Day) | (GHz) | Array | (mJy)bbThe errors on the flux densities are the sum of error from task IMFIT and a 10% calibration uncertainity added in quadrature.
2017 Aug 16.44 | 330.29 | 2.499 | C | 0.992 $\pm$ 0.108
| 330.29 | 3.499 | C | 0.729 $\pm$ 0.094
| 330.29 | 4.999 | C | 0.620 $\pm$ 0.082
| 330.29 | 7.099 | C | 0.330 $\pm$ 0.054
| 330.29 | 8.549 | C | 0.290 $\pm$ 0.044
| 330.29 | 9.499 | C | 0.270 $\pm$ 0.035
2018 Jan 25.04 | 491.89 | 2.499 | B | 0.639 $\pm$ 0.075
| 491.89 | 3.499 | B | 0.436 $\pm$ 0.050
| 491.89 | 4.999 | B | 0.367 $\pm$ 0.042
| 491.89 | 7.099 | B | 0.282 $\pm$ 0.033
| 491.89 | 8.549 | B | 0.205 $\pm$ 0.027
| 491.89 | 9.499 | B | 0.200 $\pm$ 0.024
2018 Jan 26.04 | 492.89 | 1.749 | B | 1.100 $\pm$ 0.112
| 492.89 | 13.499 | B | 0.092 $\pm$ 0.014
| 492.89 | 15.999 | B | 0.082 $\pm$ 0.011
Figure 1:  The GMRT and JVLA light curves of SN 2016gkg at frequencies 0.33–25
GHz. The inverted triangle denotes 3$\sigma$ upper limits. The age of the SN
is calculated assuming the date of explosion to be 2016 Sep 20.15 UT
(Kilpatrick et al., 2017).

##  4 Radio emission Model

The general properties of radio emission from CCSNe have been discussed in
detail by Chevalier (1982, 1998); Weiler et al. (1986, 2002). The radio light
curves and spectra can be described using the “mini-shell” model (standard
model; Chevalier, 1982). According to this model, the forward shock from the
SN interacts with the ionized CSM established due to the stellar winds of the
progenitor star. At the shock, particles are accelerated to relativistic
velocities in amplified magnetic fields and emit synchrotron radiation. A
fraction of post-shock energy density is distributed into magnetic fields
($\epsilon_{\rm B}$) and relativistic electrons ($\epsilon_{\rm e}$), and this
fraction is assumed to be constant throughout the evolution of the ejecta. The
observed radio light curves/spectra will be characterized by synchrotron
radiation, where the low-frequency emission is significantly suppressed by an
absorption component. The absorption can be free-free absorption (FFA) due to
the ionized wind material along the line of sight (Weiler et al., 1986) or due
to the same relativistic electrons that generate radio emission (synchrotron
self-absorption; SSA Chevalier, 1998). The radio flux density initially rises
rapidly and then declines as a result of the combined effects of non-thermal
synchrotron emission and various absorption processes.

This standard model of hydrodynamic evolution of ejecta assumes self-similar
evolution of physical parameters across the shock (Chevalier, 1996). The shock
radius evolves as $R\sim$ $t^{m}$ where $m=(n-3)/(n-s)$, and $n$ indicates the
outer density profile of the SN ejecta ($\rho_{\rm ejecta}\propto r^{-n}$).
$s$ denotes the density profile of CSM ($\rho_{\rm CSM}\propto r^{-s}$) and
the value of $s=2$ for a wind stratified medium.

We adopt the model from Weiler et al. (1986) in a scenario where the dominant
absorption process is FFA and follow a similar method as discussed in Nayana
et al. (2018) and Nayana & Chandra (2020). In an FFA model, the spectral and
temporal evolution of radio flux densities $F(\nu,t)$ can be described as;

$F(\nu,t)=K_{1}\left(\frac{\nu}{5\,\rm
GHz}\right)^{\alpha}\left(\frac{t}{10\,\rm days}\right)^{\beta}e^{-\tau_{\rm
FFA}}$ (1)

$\tau_{\rm FFA}=K_{2}\left(\frac{\nu}{5\,\rm
GHz}\right)^{-2.1}\left(\frac{t}{10\,\rm days}\right)^{\delta}$ (2)

The multi-frequency radio flux density evolution in the case of a dominant SSA
scenario can be modeled as (Chevalier, 1998)

$F(\nu,t)=K_{1}\left(\frac{\nu}{5\,\rm GHz}\right)^{2.5}\left(\frac{t}{10\,\rm
days}\right)^{a}\left[1-e^{-\tau_{\rm SSA}}\right]$ (3)

$\tau_{\rm SSA}=K_{2}\left(\frac{\nu}{5\,\rm
GHz}\right)^{-(p+4)/2}\left(\frac{t}{10\,\rm days}\right)^{-(a+b)}$ (4)

In the above equations, $K_{1}$ and $K_{2}$ denote the flux density and
optical depth at 5 GHz on $t=$ 10 days post-explosion, respectively. In
equations 1 and 2, $\alpha$ represents the spectral index
($F_{\nu}\propto\nu^{\alpha}$) and $\beta$ represents the temporal index of
radio flux densities. The term exp(${-\tau_{\rm FFA}}$) corresponds to the
attenuation due to the absorption by a unifrormly distributed ionized CSM
external to the radio emitting region. ’$\delta$’ denotes the temporal
evolution of $\tau_{\rm FFA}$, and is related to $\alpha$ and $\beta$ as
$\delta=\alpha-\beta-3$. Assuming the CSM is created due to a steady stellar
wind ($\rho_{\rm CSM}\propto$ r-2), the shock decceleration parameter; $m$ can
be written as $m=-\delta$/3 where $m=(n-3)/(n-2)$.

In equations 3 and 4, $a$ and $b$ denote the temporal indices of flux
densities in the optically thick ($F\propto t^{a}$) and thin phase ($F\propto
t^{-b}$), respectively. $\tau_{\rm SSA}$ is the optical depth due to SSA and
$p$ is the electron energy power-law index ($N(E)\propto E^{-p}$) which is
related to $\alpha$ as $p=2\alpha-1$. In an SSA model ‘m’ is connected to a,
b, and p as $a=2m+0.5$ in the optically thick phase and $b=(p+5-6m)/2$ in the
optically thin phase (Chevalier, 1998).

We perform a two variable fit to the entire data to find the best parameter
fit to FFA model using equations 1 and 2 and SSA model using equations 3 and
4. We execute the fit adopting the Markov chain Monte Carlo (MCMC) method
using python package emcee (Foreman-Mackey et al., 2013). We choose 32 walkers
and 5000 steps to explore the parameter space to get the best-fit values (68%
confidence interval). We estimate the goodness of fit using the reduced
$\chi^{2}_{\mu}$ test. We allow the parameters $K_{1}$, $K_{2}$, $\alpha$,
$\beta$, and $\delta$ to vary freely in the FFA model, and $K_{1}$, $K_{2}$,
$a$, $b$, and $p$ to vary freely in the SSA model. The best fit values of
these parameters are listed in Table 4. The best-fit modeled curves along with
the observed data points are shown in Figures 2 and 3. The corner plots are
presented in Figure 4.

From the best-fit modeled curves and reduced $\chi_{\mu}^{2}$ values, the SSA
model ($\chi_{\mu}^{2}=$ 3.2) seems like a better representation of the
observed data compared to FFA model ($\chi_{\mu}^{2}=$ 4.4). However, there
are deviations from SSA model predictions, particularly at $t\sim$ 299 days.
The model underestimates the flux densities at multiple frequencies (see
Figure 2). This effect is clear in the spectrum of day 299 where all the flux
densities in the optically thin regime are systematically above the modeled
curve (see Figure 3). We discuss this effect in the context of a CSM density
enhancement in §5.

The best-fit value of shock deceleration parameter, $m\sim 0.6$
($m=-\delta/3$) in FFA model. This is fairly low compared to the typical $m$
values seen in SNe IIb and would imply a highly decelerating shock wave which
is unphysical at this early stages of evolution. The value of $m\sim$ 1 from
SSA model [$m=(p+5-2b)/6$], indicative of a non-decelerating blast wave. Thus
we infer SSA model to be a better representation of the data over FFA model
due to lower $\chi_{\mu}^{2}$ values and the unrealistic $m$ value implied by
the FFA model.

The low-frequency flux measurements at earlier epochs (t $\sim$ 24 and 49
days) are above the model predictions (see Figure 3). The spectral indices of
the flux densities between 0.61 and 1.39 GHz are $0.76\pm 0.45$ and $0.59\pm
0.25$ at $t\sim$ 24 and 49 days, respectively. These values are flatter than
the expected spectral indices ($\alpha=$ 2.5) in a standard SSA model. This
can be attributed to the inhomogeneities in the magnetic fields and/or
relativistic electron distribution in the emitting region (Björnsson &
Keshavarzi, 2017; Chandra et al., 2019; Ho et al., 2019; Nayana & Chandra,
2021).

Table 4: Best fit parameters for FFA and SSA models.  FFA |  SSA
---|---
$K_{1}=13.44_{-0.99}^{1.12}$ |  $K_{1}=0.06_{-0.00}^{0.00}$
$K_{2}=11.97_{-0.58}^{0.63}$ |  $K_{2}=207.04_{-20.80}^{25.00}$
$\alpha=-1.04_{-0.03}^{0.03}$ |  $a=2.33_{-0.04}^{0.04}$
$\beta=-0.93_{-0.02}^{0.02}$ |  $b=0.94_{-0.02}^{0.02}$
$\delta=-1.88_{-0.03}^{0.02}$ |  $p=3.03_{-0.05}^{0.05}$
$\chi^{2}=4.4$ |  $\chi^{2}=3.2$
Figure 2:  Radio light curves of SN 2016gkg at frequencies $\nu=$ 0.33 $-$ 24
GHz. The solid red curves represent the best-fit SSA (equations 3 and 4) model
and solid green curves represent best-fit FFA (equations 1 and 2) model (see
§4). The light green (FFA) and red (SSA) lines represent 100 random draws from
the MCMC posterior. The filled blue circles denote the JVLA flux density
measurements and the filled cyan squares are GMRT flux density measurements.
The inverted triangles denote 3$\sigma$ flux density upper limits. The fit is
a 2D fit performed by including the entire data set. The age of the SN is
calculated assuming the date of explosion to be 2016 Sep 20.15 UT (Kilpatrick
et al., 2017). Figure 3:  Radio spectra of SN 2016gkg at $t\sim$ 24, 49, 86,
299, 330, and 492 days post-explosion. The solid red curves represent the
best-fit SSA (equations 3 and 4) model and solid green curves represent best-
fit FFA (equations 1 and 2) model (see §4). The light green (FFA) and red
(SSA) lines represent 100 random draws from the MCMC posterior. The filled
blue circles denote the JVLA flux density measurements and the filled cyan
squares are GMRT flux density measurements. The fit is a 2D fit performed by
including the entire data set. The age of the SN is calculated assuming the
date of explosion to be 2016 Sep 20.15 UT (Kilpatrick et al., 2017).

Figure 4:  Corner plots that show the results of MCMC modeling of the radio
data with the SSA model (equations 3 and 4) in top panel and FFA model
(equations 1 and 2) in bottom panel as discussed in §4. The 16, 50, and 84
percentiles are marked in each of the contours and histograms.

###  4.1 Single epoch spectral analysis

To further investigate the time evolution of shock parameters, we model each
of the single epoch spectra adopting the SSA model. The functional form of the
SSA spectrum can be parametrized as (Soderberg et al., 2006)

$F_{\nu}=1.582\,F_{p}\left(\frac{\nu}{\nu_{p}}\right)^{5/2}\left(1-\rm
exp\left[-\left(\frac{\nu}{\nu_{p}}\right)^{-(5-2\alpha)/2}\right]\right)$ (5)

We model single epoch spectra by allowing peak frequency ($\nu_{p}$), peak
flux density ($F_{p}$), and $\alpha$ to vary freely and independently. The
spectra are well fitted by an SSA model with $\chi_{\mu}^{2}$ values of 0.7
$-$ 3.3 (see Figure 5 and Table 5). We attempted modeling the data with FFA
model as well and the fits resulted in higher $\chi_{\nu}^{2}$ values.

The best fit $F_{p}$, $\nu_{p}$, and $\alpha$ for SSA model are given in Table
5. The quoted errors are 1$\sigma$ errors (68% confidence interval). The
temporal evolution of $\nu_{\rm p}$ and $F_{\rm p}$ are such that $\nu_{\rm
p}\propto t^{-0.9}$ and $F_{\rm p}\propto t^{-0.02}$, respectively, consistent
with SSA model (Chevalier, 1998). We obtain the average spectral index over
six epochs as $\alpha\approx-1$. The power-law index of electrons determined
from the optically thin spectral index ($\alpha$) is $p=3$
($\alpha=-(p-1)/2$).

###  4.2 Blast-wave parameters

The shock radius ($R$) and magnetic fileds ($B$) can be estimated from
$\nu_{\rm p}$ and $F_{\rm p}$ at each epoch (Chevalier, 1998). For $p=3$, the
shock radius is given by

$\begin{split}R=8.8\times 10^{15}f_{\rm
eB}^{-1/19}\left(\frac{f}{0.5}\right)^{-1/19}\left(\frac{F_{p}}{\rm
Jy}\right)^{9/19}\\\ \left(\frac{D}{\rm
Mpc}\right)^{18/19}\left(\frac{\nu_{p}}{5\,\rm GHz}\right)^{-1}\rm
cm\end{split}$ (6)

The post-shock magnetic field is given by

$\begin{split}B=0.58f_{\rm
eB}^{-4/19}\left(\frac{f}{0.5}\right)^{-4/19}\left(\frac{F_{p}}{\rm
Jy}\right)^{-2/19}\\\ \left(\frac{D}{\rm
Mpc}\right)^{-4/19}\left(\frac{\nu_{p}}{5\,\rm GHz}\right)\rm G\end{split}$
(7)

Here, $f_{\rm eB}$ denotes the ratio of the fraction of shock energy in
relativistic electrons ($\epsilon_{e}$) to that in the magnetic fields
($\epsilon_{B}$). We assume the equipartition of energy between relativistic
electrons and magnetic fields and hence use $f_{\rm eB}=1$. $f$ is the volume
filling factor of the synchrotron emitting region, taken to be 0.5 (Chevalier,
1998). $D$ is the distance to the SN in Mpc. The mean velocity of the shock at
any epoch is $V$ $\sim$ $\frac{R}{t}$.

The mass-loss rates can be deduced using the magnetic field scaling relation
(Chevalier, 1998).

$\dot{M}=\frac{B^{2}R^{2}v_{w}}{2\epsilon_{B}V^{2}}$ (8)

We deduce the shock radius ($R$) and magnetic fileds ($B$) at multiple epochs
using the best fit $\nu_{\rm p}$ and $F_{\rm p}$ values at $t\sim$ 24, 49, 86,
299, 330, and 492 days using equations 6 and 7. The physical parameters are
presented in Table 5. The shock wave expands from $R\sim$ 0.5 $\times$ 1016 cm
to $R\sim$ 7.3 $\times$ 1016 cm during $t\sim$ 24 to 492 days. The temporal
evolution of shock radius can be described by an index $m=$ 0.80 $\pm$ 0.08,
indicating a deccelerating blast wave (see Fig 6). The shock is slightly
slowing down at $t\sim$ 299 days and later evolve consistent with the temporal
evolution as described by the previous phases ($t<$ 299 days). The temporal
index for post-shock magnetic field is found to be $\alpha_{\rm B}=$ $-$0.80
$\pm$ 0.09 ($B\propto t^{\alpha_{\rm B}}$).

The mass-loss rate of the progenitor star at different epochs are estimated
using equation 8. We assume a wind velocity $v_{\rm w}\sim$ 200 km s-1 and
$\epsilon_{\rm B}=0.33$ in the calculation. The mass-loss rates are
$\dot{M}\sim$ (2.2 $-$ 5.0) $\times$ 10-6 $M_{\odot}$ yr-1 during 15 $-$ 115
years before explosion (see Table 5). Besides, we note that the mass-loss rate
derived from the shock parameters at $t\sim$ 299 days is relatively high,
$\dot{M}\sim$ 12.6 $\times$ 10-6 $M_{\odot}$ yr-1, indicating a higher mass-
loss at $\sim$ 48 years prior explosion. Considering the uncertainty due to
the low temporal cadence of the radio observations, the timing of the enhanced
mass-loss event could be between 25$-$87 years prior explosion. Hydrodynamic
wave driven outbursts during the late nuclear burning stage can create density
enhancements in the CSM. However, these outbursts happen at 1-2 years prior
explosion in the case of SNe IIb (Fuller & Ro, 2018), which does not match
with the timescales of enhanced mass-loss in SN 2016gkg.

The $\dot{M}$ estimates will considerably vary depending on the choice of wind
velocity. We choose $v_{\rm w}\sim$ 200 km s-1 based on the physical
properties of the progenitor star ($M\sim$ 10 $M_{\odot}$, $R_{*}\sim$ 70
$R_{\odot}$, $T_{\rm eff}\sim$ 10800 K) derived from pre-explosion imaging
analysis (Kilpatrick et al., 2021). The escape velocity of a 10 $M_{\odot}$
star of radius 70 $R_{\odot}$ is $\sim$ 200 km s-1 and the range of wind
velocity of a star of $T_{\rm eff}\sim$ 10800 K is $v_{\rm w}$ $\sim$
100$-$300 km s-1 (Drout et al., 2009; Smith, 2014; Yoon et al., 2017).

Table 5: Blast wave paramters of SN 2016gkg

Agea | $F_{\rm p}$ | $\nu_{\rm p}$ | $\alpha$ | $\chi_{\mu}^{2}$ | $R$ | $B$ | $V_{\rm s}$ |  $\dot{M}$
---|---|---|---|---|---|---|---|---
(Day) | (mJy) | (GHz) | - | - | ($\times$1015cm) | (Gauss) |  $(\times$104km s${}^{-1})$ | ($\times$10-6 $M_{\odot}$ yr-1)
24 | $1.55_{-0.16}^{0.18}$ | $8.96_{-0.65}^{0.68}$ | $-0.80_{-0.22}^{0.20}$ | 1.72 | $5.09_{-0.66}^{0.68}$ | $1.03_{-0.08}^{0.08}$ | $2.46_{-0.32}^{0.33}$ | $2.20_{-0.86}^{0.90}$
49 | $2.20_{-0.14}^{0.14}$ | $5.82_{-0.25}^{0.26}$ | $-1.49_{-0.13}^{0.13}$ | 3.30 | $9.26_{-1.00}^{1.01}$ | $0.65_{-0.03}^{0.03}$ | $2.19_{-0.24}^{0.24}$ | $3.58_{-1.16}^{1.16}$
86 | $2.22_{-0.14}^{0.14}$ | $3.43_{-0.23}^{0.28}$ | $-1.25_{-0.09}^{0.09}$ | 0.66 | $15.77_{-1.90}^{2.02}$ | $0.38_{-0.03}^{0.03}$ | $2.13_{-0.26}^{0.27}$ | $3.82_{-1.40}^{1.52}$
299 | $2.28_{-0.19}^{0.19}$ | $1.79_{-0.13}^{0.12}$ | $-1.19_{-0.09}^{0.09}$ | 0.98 | $30.59_{-3.81}^{3.70}$ | $0.20_{-0.01}^{0.01}$ | $1.18_{-0.15}^{0.14}$ | $12.56_{-4.80}^{4.62}$
330 | $1.62_{-0.14}^{0.15}$ | $0.85_{-0.06}^{0.06}$ | $-0.92_{-0.07}^{0.07}$ | 0.72 | $55.08_{-6.93}^{6.88}$ | $0.10_{-0.01}^{0.01}$ | $1.93_{-0.24}^{0.24}$ | $3.68_{-1.42}^{1.42}$
492 | $1.79_{-0.14}^{0.15}$ | $0.67_{-0.06}^{0.06}$ | $-1.07_{-0.05}^{0.05}$ | 1.63 | $72.57_{-10.00}^{9.51}$ | $0.08_{-0.01}^{0.01}$ | $1.71_{-0.24}^{0.22}$ | $5.02_{-2.20}^{2.06}$

††footnotetext: $F_{\rm p}$, $\nu_{\rm p}$ and $\alpha$ denote the peak flux
density, peak frequency and optically thin spectral index in the best fit SSA
model at each epochs. $R$, $B$, $V$, and $\dot{M}$ denote shock radius,
magnetic field, mass-loss rate and mean shock velocity, respectively.

##  5 Non-uniforn density of the CSM

The overall evolution of radio light curves and spectra of SN 2016gkg is best
modeled by a self-absorbed synchrotron emission that arises due to the
interaction of SN shock with the CSM created by a uniform mass-loss from the
progenitor. However, there are some deviations from the smooth light
curve/spectra evolution as prescribed by the standard model. There is a
fractional increase by a factor of $\sim$ 2 in flux densities (in the
optically thin phase) at different frequencies on day 299 above the model
prediction (see Fig 2 and 3). This abrupt rise in flux densities could be due
to the interaction of the forward shock with density enhancements in the CSM
at a radius $R\sim$ 3.1 $\times$ 1016 cm. These density fluctuations could be
either due to the non-uniform mass-loss rate of a single star progenitor via
stellar winds and/or due to the mass stripping by a binary companion
(Soderberg et al., 2006). There are several observational pieces of evidence
from supernova remnants and massive stars that support complex mass-loss
events happening towards the end stages of stellar evolution (Soderberg et
al., 2006). In the case of a binary scenario, the strength and position of CSM
density enhancement will be influenced by the binary parameters (Podsiadlowski
et al., 1992). A build-up of CSM material can happen at particular spatial
scales due to the modulation of progenitor stellar wind depending on the
orbital period of the binary companion (Weiler et al., 1992) and the
eccentricity of the binary orbit. A binary scenario with a period of 4000
years has been attributed to the periodic modulations in the radio light
curves of SN 1979C (Weiler et al., 1992; Montes et al., 2000).

Multiple episodes of mass-loss events have been attributed to periodic light
curve bumps of modest (factor of $\sim$ 2) flux density fluctuations in SN
2003bg (Soderberg et al., 2006) and SN 2001ig (Ryder et al., 2004). Both these
SNe showed variations in the light curve during a period $t=$120$-$300 days
(see Fig 8), indicating a radial distance of 4$\times$1016 cm to 8
$\times$1016 cm from the explosion center, similar to that of SN 2016gkg. We
note that the temporal cadence of the follow-up observations of these SNe are
good enough to map the periodic bumps in their light curves (see Fig 8).

In the case of SN 2016gkg, we see flux density enhancement at $t\sim$ 299
days, that correspond to a stellar evolution phase of $\sim$ 48 years prior
explosion for $v_{\rm w}\sim$ 200 km s-1. Kilpatrick et al. (2017) suggested a
binary progenitor model for SN 2016gkg where the initial period of the binary
orbit is $\sim$ 1000 days (2.7 years). This periodicity will be seen as
another flux density enhancement at $t\sim$ 313 days for the derived shock
velocities. This epoch is not sampled in the observations. Thus even if there
is a binary companion and related periodicity in the density distribution of
CSM, the radio data do not have a temporal cadence to probe those
fluctuations, and we cannot rule out the binary scenario.

The radio luminosity is related to the density of CSM as $L\propto$ $\rho_{\rm
CSM}^{(p-7+12m)/4}$, where $\rho_{\rm CSM}\propto(\dot{M}/v_{\rm w})$ (Ryder
et al., 2004). In the case of SN 2016gkg, the best fit value of $p=3$ and
$m=0.8$ indicates that a factor of 2 increase in radio flux density indicates
$\sim$70% increase in the CSM density. The effect of density enhancement is
reflected in the evolution of shock radii and magnetic field as well (see Fig
6). The magnetic field is slightly enhanced compared to its expected regular
temporal evolution and there is a slight change in the expansion of shock wave
at the same time. The magnetic field is enhanced by a factor of 1.6 in
comparison to the extrapolation of its evolution in the previous phases. The
additional thermal energy produced due to these density enhancements increases
the $B$. A similar increase in $B$ value by a factor of 1.3 is seen in SN
2003bg during its first light curve bump (Soderberg et al., 2006). To
summarize, there is a close resemblance between SN 2016gkg, and SN 2003bg, and
SN 2001ig in terms of the timescale and strength of flux density enhancement
in late-time radio light curves.

##  6 Discussion

The shock velocities of SN 2016gkg derived from SSA modeling is $v\sim$ 24600
km s-1 ($v\sim$ 0.1 c) at $t\sim$ 24 days. The velocities estimated from
optical lines at $t\sim$ 21 days is $v\sim$ 12,200 km s-1 (Tartaglia et al.,
2017). Thus the shock wave is traveling with a velocity a factor of 2 faster
than the material in the photosphere. The SSA-derived shock velocity is
greater than the velocities from optical lines and indicates that FFA is not
contributing much in defining the peak of the light curve/spectra. We also
note that the FFA models were resulting in higher $\chi_{\mu}^{2}$ values
while modeling. The radio data being inconsistent with the FFA model and
relatively higher shock velocity of $v\sim$ 0.1 c are indicative of a compact
progenitor star with faster stellar winds.

The temporal evolution of shock wave radius is best fitted by a power law
$R\propto$ $t^{0.80\pm 0.08}$. The $m$ value will be 1, for a non-
deccelerating blast wave and the derived $m$ value indicates a decelerating
blast-wave. The temporal index of post-shock magnetic field is found to be
$\alpha_{\rm B}$ = $-$0.8$\pm$0.1 ($B\propto t^{\alpha_{\rm B}}$).
$\alpha_{\rm B}$ can be connected to the CSM density as $\alpha_{\rm
B}=[m(2-s)/2]-1$, which gives $s=1.5$. Thus the CSM density is slightly
flatter than the one created by a steady stellar wind. The derived values of
$m$ and $s$ indicate an ejecta density profile of $\rho_{\rm ejecta}\propto
r^{-9}$ ($m=n-3/n-s$), consistent with low-mass compact progenitor (Chevalier,
1998). Assuming $\epsilon_{\rm B}=$ 0.33 and $v_{\rm w}=$ 200 km s-1, the
mass-loss rate of the progenitor is in the range $\dot{M}\sim$ (2.2–5.0)
$\times$ 10-6 $M_{\odot}\,yr^{-1}$, at $t\sim$ 24, 49, 86, 330, and 492 days.
The mass-loss rate corresponding to $t\sim$ 299 days is $\dot{M}\sim$ 12.6
$\times$ 10-6 $M_{\odot}\,yr^{-1}$, a factor of three higher than the
$\dot{M}$ values at other epochs. This is suggestive of an enhanced phase of
mass-loss at $\sim$ 48 years prior explosion for the assumed wind speed as
discussed in §5.

Figure 5:  Radio spectra of SN 2016gkg at $t\sim$ 24, 49, 86, 299, 330, and
492 days post-explosion. The filled circles denote the JVLA flux density
measurements and the filled squares are near-simultaneous flux density
measurements from the GMRT. The solid curves are the best-fit SSA models at
each epoch. Figure 6:  The shock radii (green filled circles) and magnetic
fields (blue filled circles) at $t\sim$ 24, 49, 86, 299, 330, and 492 days
post-explosion. The green solid line denotes a power-law fit to the shock
radii; $R\sim$ $t^{0.80\pm 0.08}$ and blue solid line denotes a power-law fit
to the magnetic field; $B\sim$ $t^{-0.80\pm 0.09}$. Figure 7: The peak radio
spectral luminosities versus time of peak of the light curves are plotted for
all well-observed SNe IIb in the literature. Each SNe is denoted by the last
two digits of the year of explosion and last two letters of the name. The
dotted lines represent the mean velocity of the radio emitting shell if SSA is
the dominant absorption process that defines the peak of the light curve
(Chevalier, 1998) and for $p=$ 3\.  Table 6: Properties of well-studied type
IIb supernovae

SN | Distance | $\nu_{p}$ | $t_{p}$ | $F_{p}$ | $L_{\rm p}$ | $\dot{M}$ | References
---|---|---|---|---|---|---|---
- | (Mpc) | (GHz) | (days) | (mJy) | (erg s-1 Hz-1) | ($M\odot yr^{-1}$) | -
SN 1993J | 3.6 | 5 | 133 | 96.9 | 1.5 $\times$ 1027 | (2$-$6) $\times$10-5 | 1
SN 1996cb | 9.1 | 5 | 19.4 | 1.8 | 1.8 $\times$ 1026 | – | 2
SN 2001gd | 17.5 | 5 | 173 | 8.0 | 2.9 $\times$ 1027 | 3$\times$10-5 | 3
SN 2001ig | 11.5 | 5 | 74 | 22 | 3.5 $\times$ 1027 | 8.6$\times$10-5 | 4
SN 2003bg | 19.6 | 22.5 | 35 | 85 | 3.9 $\times$ 1028 | 6.1 $\times$ 10-5 | 5
SN 2008ax | 6.2 | 4.86 | 15.34 | 3.54 | 1.6 $\times$ 1026 | (1-6) $\times$ 10-6 | 6
SN 2008bo | 19.1 | 8.5 | 11.6 | 0.52 | 2.3 $\times$ 1026 | - | 7
SN 2010P | 44.8 | 5 | 464 | 0.52 | 1.2 $\times$ 1027 | (3.0$-$5.1)$\times$10-5 | 8
SN 2010as | 27.4 | 9 | 34.3 | 2.43 | 2.2 $\times$1027 | - | 7
SN 2011dh | 7.9 | 4.7 | 35.3 | 7.3 | 5.4 $\times$ 1026 | 6 $\times$10-5 | 9
SN 2011hs | 26.4 | 5 | 59 | 2.0 | 1.7 $\times$ 1027 | 2 $\times$ 10-5 | 10
SN 2013df | 16.6 | 5 | 67.3 | 1.55 | 5.1 $\times$1026 | 8 $\times$ 10-5 | 11
SN 2016bas | 42.4 | 5.5 | 148.8 | 4.33 | 9.3 $\times$1027 | - | 7
SN 2016gkg | 26.4 | 5 | 85.86 | 1.92 | 1.6 $\times$1027 | 3.8$\times$10-6 | 12

he listed $\dot{M}$ values are taken from the literature as cited in the last
column of the table. These values are strongly dependent on the assumed wind
velocities.

References. — (1) Fransson et al. (1996), (2) Weiler et al. (1998), (3)
Stockdale et al. (2003), (4) Ryder et al. (2004), (5) Soderberg et al. (2006),
(6) Roming et al. (2009), (7) Bietenholz et al. (2021), (8) Romero-Cañizales
et al. (2014), (9) Soderberg et al. (2012), (10) Bufano et al. (2014), (11)
Kamble et al. (2016), (12) This work

Note. — Among the listed SNe, five of them have progenitor detections from
archival images. They are SN 1993J (Aldering et al., 1994), SN 2013df (Van Dyk
et al., 2014), SN 2008ax (Folatelli et al., 2015), SN 2011dh (Van Dyk et al.,
2013; Arcavi et al., 2011) and SN 2016gkg (Kilpatrick et al., 2017, 2021;
Tartaglia et al., 2017).

Note. — T

Figure 8: 5 GHz radio light curves of well-observed SNe IIb. Flux measurements are taken from Weiler et al. (2007) for SN 1993J, Soderberg et al. (2006) for SN 2003bg, Ryder et al. (2004) for SN 2001ig, Roming et al. (2009) for SN 2008ax, Bietenholz et al. (2021) for SN 2016bas, Kamble et al. (2016) for SN 2013df, Krauss et al. (2012) for SN 2011dh, Bufano et al. (2014) for SN 2011hs, Stockdale et al. (2007) for SN 2001gd, Romero-Cañizales et al. (2014) for SN 2010P, and This work for SN 2016gkg. Table 7: Compilation of the radius estimates of the progenitor of SN 2016gkg from the literature Radius ($R_{\odot}$) | Method | Reference
---|---|---
138${}^{+131}_{-103}$ | Archival HST imaging analysisaaThe age is calculated using 2016 Sep 20.15 (UT) as the date of explosion (Kilpatrick et al., 2017). | Kilpatrick et al. (2017)
150 $-$ 320 | Archival HST imaging analysisaaKilpatrick et al. (2017) and Tartaglia et al. (2017) used Keck and Very Large Telescope + Nasmyth adaptive optics systems, respectively to perform relative astrometry. Kilpatrick et al. (2017) considered one progenitor candidate and Tartaglia et al. (2017) considered two progenitor candidates.  | Tartaglia et al. (2017)
$\sim$ 70 | Archival HST imaging analysisbbPost-explosion HST imaging of the field containing SN 2016gkg was used to perform relative astrometry. | Kilpatrick et al. (2021)
257${}^{+389}_{-189}$ | Analytical shock cooling modelccKilpatrick et al. (2017) and Tartaglia et al. (2017) modeled the luminosity (up to $t\sim$ 1.5 days) and temperature evolution (up to $t\sim$ 5 days), respectively using Rabinak & Waxman (2011) model. | Kilpatrick et al. (2017)
48 $-$ 124 | Analytical shock cooling modelccKilpatrick et al. (2017) and Tartaglia et al. (2017) modeled the luminosity (up to $t\sim$ 1.5 days) and temperature evolution (up to $t\sim$ 5 days), respectively using Rabinak & Waxman (2011) model. | Tartaglia et al. (2017)
40 $-$ 150 | Analytical shock cooling modelddAnalytic shock cooling models (Piro, 2015; Nakar & Piro, 2014; Sapir & Waxman, 2017). | Arcavi et al. (2017)
180 $-$ 260 | Numerical shock cooling model | Piro et al. (2017)
$\sim$ 320 $R_{\odot}$ | Numerical shock cooling model | Bersten et al. (2018)

###  6.1 Comparison with other radio bright SNe IIb

We view the radio properties of SN 2016gkg in comparison with other radio
bright SNe IIb events in Figure 7. We compile all SNe IIb with well sampled
light curves/spectra that defines $F_{\rm p}$, $\nu_{\rm p}$ and $t_{\rm p}$
in the $L_{\rm p}$ \- $t_{\rm p}$ diagram (also see Table 6). The dotted lines
indicate the mean shock velocities in an SSA scenario for $p=3$ assuming
equipartition of energy between relativistic particles and magnetic fields
($\epsilon_{\rm B}$ = $\epsilon_{\rm e}$). This plot is an updated version of
a similar plot presented in Chevalier & Soderberg (2010). The authors compiled
radio properties of a sample of SNe IIb and divided them into two populations
based on their position in the $L_{\rm p}$ \- $t_{\rm p}$ diagram, one is SNe
IIb with compact progenitors (SNe cIIb) and the other with extended
progenitors (SNe eIIb). The SNe cIIb group consists of SN 2008ax, SN 2003bg
and SN 2001ig which has faster shocks, less dense CSM and compact progenitor
in comparison with the SNe eIIb (like SN 1993J and SN 2001gd). SNe eIIb have
slower shocks owing to its denser CSM from slow stellar winds of extended
progenitors.

SNe IIb show peak spectral luminosities that span 2 orders of magnitude, and
the peak time vary over a factor of $\sim$ 40\. This broad distribution in
peak spectral luminosities and rise times indicates the variety in the
intrinsic properties of their progenitors. There are only 5 SNe IIb (including
SN 2016gkg) that have progenitors identified from pre-explosion images. It is
important to tie up the radio properties of these SNe with the inferences from
pre-explosion imaging analysis. The progenitor of SN 1993J is identified to be
a YSG of radius $\sim$ 600 $R_{\odot}$ (Aldering et al., 1994) from pre-
explosion images. Similarly, an extended progenitor of radius $\sim$ 545 $\pm$
65 $R_{\odot}$ was identified as the progenitor of SN 2013df (Van Dyk et al.,
2014). A slightly less extended star of radius $\sim$ 200 $R_{\odot}$ was
identified as the progenitor of SN 2011dh (Van Dyk et al., 2013). However, the
presence of a binary companion has been speculated for this SN (Van Dyk et
al., 2013), and Arcavi et al. (2011) suggested the progenitor of SN 2011dh to
be a relatively compact star from the analysis of a series of spectra and
bolometric light curve. The authors also argued that the larger radius
($R\sim$ 1013 cm) derived from pre-explosion HST images could be due to the
identification of a blended source. The progenitor of SN 2008ax was identified
to be of radius $\sim$ 30$-$50 $R_{\odot}$ (Folatelli et al., 2015). These
estimates of progenitor radius from pre-SN image analysis is roughly
consistent with the radio-derived properties. SN 1993J and SN 2013df to have
more extended progenitors whereas SN 2008ax and SN 2011dh to have relatively
compact progenitors with higher shock velocity (see Figure 7). In light of
these results, one can also argue that the classification of SNe IIb
progenitors into two categories as eSNe IIb and cSNe IIb (Chevalier &
Soderberg, 2010) is rather simplistic and the progenitor properties could be a
continuum between these two. The position of SN 2016gkg in this diagram is
among SNe cIIb, towards the right of SN 2008ax. This could imply that the
progenitor of SN 2016gkg is a relatively compact progenitor with a radius
slightly more than that of SN 2008ax, i.e., 50 $R_{\odot}$ .

###  6.2 Inferences on the progenitor

Multiple pieces of evidence from radio modeling are in favor of a compact
progenitor for SN 2016gkg. The broad agreement of multi-frequency radio data
with an SSA model indicates that the CSM is relatively rarer created due to
faster stellar winds from a compact star. The mean shock velocities ($v\sim$
0.1 c) derived from SSA formulation are difficult to incorporate in the
framework of a shock breakout from an extended progenitor (Nakar & Sari,
2010).

A correlation between $\dot{M}$ values and progenitor radius of SNe IIb is
proposed by Maeda et al. (2015) and Kamble et al. (2016), where the extended
progenitors experience stronger mass-loss towards their end-of-life compared
to compact progenitors. The progenitor mass-loss rate of SN 2016gkg derived
from the shock parameters at $t\sim$ 24 days is $\dot{M}\sim$ 2.2 $\times$
10-6 $M_{\odot}$ yr-1, comparable to the $\dot{M}$ values derived for SN
2008ax (Roming et al., 2009). These $\dot{M}$ values are an order of magnitude
lower than that of SN 1993J and SN 2013df (see Table 6) which are known to
have extended progenitors from direct detection efforts. Thus the $\dot{M}$
estimates also imply a relatively compact progenitor.

The best-fit value of electron power-law index is $p\sim$ 3, typically found
for SNe Ibc which are presumed to have compact WR stars as progenitors. Late-
time variability in the radio light curves is an important observational
characteristic of SNe cIIb progenitors (e.g., SN 2001ig, SN 2003bg Ryder et
al., 2004; Soderberg et al., 2006). All SNe cIIb except SN 2011dh with well-
sampled radio light curves exhibit fluctuations indicative of density
modulations in the CSM. These density fluctuations could be due to the
variability in the stellar winds of compact stars or due to the influence of a
binary companion. The radio observations of SN 2011dh probe a radius up to 1.5
$\times$ 1016 cm that translates to $\sim$ 5 years prior explosion for a
$v_{\rm w}\sim$ 1000 km s-1, which could be shorter for any substantial wind
variability (Krauss et al., 2012). We see similar late-time flux density
enhancement in the radio light curves of SN 2016gkg. The position of SN
2016gkg in the $L_{\rm p}-t_{\rm p}$ diagram is in the contour of SNe cIIb,
indicating a progenitor radius slightly more than that of SN 2008ax (i.e., $>$
50 $R_{\odot}$).

The radius estimates of the progenitor from shock cooling models span a wide
range $\sim$ 40 – 646 $R_{\odot}$ (Bersten et al., 2018; Kilpatrick et al.,
2017; Tartaglia et al., 2017; Piro et al., 2017; Arcavi et al., 2017). One can
argue that the radio derived constraints of a compact progenitor is broadly in
agreement with the results from shock cooling models. The large range of
radius estimates from these models will be in agreement with an extended
progenitor model as well.

Kilpatrick et al. (2017) and Tartaglia et al. (2017) determined the progenitor
radii to be $R=$138${}^{+131}_{-103}$ $R_{\odot}$ and $R=$150$-$320
$R_{\odot}$, respectively, from the pre-explosion HST imaging analysis of the
field containing SN 2016gkg. Relative astrometry was done using Keck and VLT
adaptive optics system in these studies. The late-time HST imaging of the
field of SN 2016gkg using Advanced Camera for Surveys (ACS) and Wide Field
Camera 3 (WFC 3) resulted in superior resolution and improved astrometric
alignment between the SN and progenitor candidate (Kilpatrick et al., 2021).
The updated photometric analysis suggests the progenitor to be a yellow
supergiant of mass 10${}^{+2}_{-1}$ $M_{\odot}$ and radius $\sim$ 70
$R_{\odot}$. Thus the inferences on the progenitor star from radio analyses is
in agreement with that derived from pre-explosion imaging analysis (Kilpatrick
et al., 2021).

##  7 Conclusions

We present long-term ($t\sim$ 8 $-$ 1429) radio monitoring of SN 2016gkg over
a frequency range $\nu\sim$ 0.3 $-$ 24 GHz to investigate the properties of
its progenitor and CSM. The inferences from our observations and modeling can
be summarized as follows

  * •

The radio data is best represented by a self-absorbed synchrotron emission
that arises due to the interaction of an SN shock wave of $v\sim$ 0.1 c
propagating into a CSM created due to the mass-loss of the progenitor star.

  * •

The CSM density is found to have moderate density fluctuation at a distance
$R\sim$ 3.1 $\times$ 1016 cm, likely due to enhancement in the progenitor
mass-loss rate or due to the effect of a binary companion. Assuming a stellar
wind velocity $v_{\rm w}\sim$ 200 km s-1, this corresponds to a stellar
evolution phase $\sim$ 48 years prior explosion.

  * •

We estimate the average mass-loss rate to be $\dot{M}\sim$ 3.7 $\times$ 10-6
$M_{\odot}\,yr^{-1}$ during 8 to 115 years before explosion, with a factor of
$\sim$ 3 higher $\dot{M}$ at $\sim$ 48 years prior explosion.

  * •

The radio data being consistent with SSA model, shock velocities of $v\sim$
0.1 c, the position of SN 2016gkg in the region of SNe cIIb in $L_{\rm
p}$-$t_{\rm p}$ diagram, and late time modest variability in radio flux
densities are suggestive of a compact progenitor star.

We thank the staff of the GMRT that made these observations possible. GMRT is
run by the National Centre for Radio Astrophysics of the Tata Institute of
Fundamental Research. The National Radio Astronomy Observatory is a facility
of the National Science Foundation operated under cooperative agreement by
Associated Universities, Inc. Nayana A.J. acknowledges DST- INSPIRE Faculty
Fellowship (IFA20-PH-259) for supporting this research. P.C. acknowledges the
support of the Department of Atomic Energy, Government of India, under project
no. 12-R&D-TFR-5.02-0700.

## References

  * Aldering et al. (1994) Aldering, G., Humphreys, R. M., & Richmond, M. 1994, AJ, 107, 662. doi:10.1086/116886
  * Arcavi et al. (2011) Arcavi, I., Gal-Yam, A., Yaron, O., et al. 2011, ApJ, 742, L18. doi:10.1088/2041-8205/742/2/L18
  * Arcavi et al. (2017) Arcavi, I., Hosseinzadeh, G., Brown, P. J., et al. 2017, ApJ, 837, L2. doi:10.3847/2041-8213/aa5be1
  * Baron et al. (1993) Baron, E., Hauschildt, P. H., Branch, D., et al. 1993, ApJ, 416, L21. doi:10.1086/187061
  * Bersten et al. (2012) Bersten, M. C., Benvenuto, O. G., Nomoto, K., et al. 2012, ApJ, 757, 31. doi:10.1088/0004-637X/757/1/31
  * Bersten et al. (2018) Bersten, M. C., Folatelli, G., García, F., et al. 2018, Nature, 554, 497. doi:10.1038/nature25151
  * Björnsson & Keshavarzi (2017) Björnsson, C.-I. & Keshavarzi, S. T. 2017, ApJ, 841, 12. doi:10.3847/1538-4357/aa6cad
  * Bietenholz et al. (2021) Bietenholz, M. F., Bartel, N., Argo, M., et al. 2021, ApJ, 908, 75. doi:10.3847/1538-4357/abccd9
  * Bufano et al. (2014) Bufano, F., Pignata, G., Bersten, M., et al. 2014, MNRAS, 439, 1807. doi:10.1093/mnras/stu065
  * Cappa et al. (2004) Cappa, C., Goss, W. M., & van der Hucht, K. A. 2004, AJ, 127, 2885. doi:10.1086/383286
  * Chandra et al. (2019) Chandra, P., Nayana, A. J., Björnsson, C.-I., et al. 2019, ApJ, 877, 79. doi:10.3847/1538-4357/ab1900
  * Chevalier (1982) Chevalier, R. A. 1982, ApJ, 259, 302. doi:10.1086/160167
  * Chevalier (1982) Chevalier, R. A. 1982, ApJ, 258, 790. doi:10.1086/160126
  * Chevalier (1992) Chevalier, R. A. 1992, ApJ, 394, 599. doi:10.1086/171612
  * Chevalier (1996) Chevalier, R. A. 1996, Radio Emission from the Stars and the Sun, 93, 125
  * Chevalier (1998) Chevalier, R. A. 1998, ApJ, 499, 810. doi:10.1086/305676
  * Chevalier & Soderberg (2010) Chevalier, R. A. & Soderberg, A. M. 2010, ApJ, 711, L40. doi:10.1088/2041-8205/711/1/L40
  * Condon et al. (1998) Condon, J. J., Cotton, W. D., Greisen, E. W., et al. 1998, AJ, 115, 1693. doi:10.1086/300337
  * Crockett et al. (2008) Crockett, R. M., Eldridge, J. J., Smartt, S. J., et al. 2008, MNRAS, 391, L5. doi:10.1111/j.1745-3933.2008.00540.x
  * Dalton & Sarazin (1995) Dalton, W. W. & Sarazin, C. L. 1995, ApJ, 448, 369. doi:10.1086/175968
  * Drout et al. (2009) Drout, M. R., Massey, P., Meynet, G., et al. 2009, ApJ, 703, 441. doi:10.1088/0004-637X/703/1/441
  * Ergon et al. (2014) Ergon, M., Sollerman, J., Fraser, M., et al. 2014, A&A, 562, A17. doi:10.1051/0004-6361/201321850
  * Filippenko (1988) Filippenko, A. V. 1988, AJ, 96, 1941. doi:10.1086/114940
  * Filippenko (1997) Filippenko, A. V. 1997, ARA&A, 35, 309. doi:10.1146/annurev.astro.35.1.309
  * Folatelli et al. (2015) Folatelli, G., Bersten, M. C., Kuncarayakti, H., et al. 2015, ApJ, 811, 147. doi:10.1088/0004-637X/811/2/147
  * Foreman-Mackey et al. (2013) Foreman-Mackey, D., Conley, A., Meierjurgen Farr, W., et al. 2013, Astrophysics Source Code Library. ascl:1303.002
  * Fransson et al. (1996) Fransson, C., Lundqvist, P., & Chevalier, R. A. 1996, ApJ, 461, 993. doi:10.1086/177119
  * Fuller & Ro (2018) Fuller, J. & Ro, S. 2018, MNRAS, 476, 1853. doi:10.1093/mnras/sty369
  * Greisen (2003) Greisen, E. W. 2003, Information Handling in Astronomy - Historical Vistas, 285, 109. doi:10.1007/0-306-48080-8_7
  * Ho et al. (2019) Ho, A. Y. Q., Phinney, E. S., Ravi, V., et al. 2019, ApJ, 871, 73. doi:10.3847/1538-4357/aaf473
  * Horesh et al. (2013) Horesh, A., Stockdale, C., Fox, D. B., et al. 2013, MNRAS, 436, 1258. doi:10.1093/mnras/stt1645
  * Hummel et al. (1987) Hummel, E., Jorsater, S., Lindblad, P. O., et al. 1987, A&A, 172, 51
  * Hummel & Jorsater (1992) Hummel, E. & Jorsater, S. 1992, A&A, 261, 85
  * Kamble et al. (2016) Kamble, A., Margutti, R., Soderberg, A. M., et al. 2016, ApJ, 818, 111. doi:10.3847/0004-637X/818/2/111
  * Kilpatrick et al. (2017) Kilpatrick, C. D., Foley, R. J., Abramson, L. E., et al. 2017, MNRAS, 465, 4650. doi:10.1093/mnras/stw3082
  * Kilpatrick et al. (2021) Kilpatrick, C. D., Coulter, D. A., Foley, R. J., et al. 2021, arXiv:2112.03308
  * Krauss et al. (2012) Krauss, M. I., Soderberg, A. M., Chomiuk, L., et al. 2012, ApJ, 750, L40. doi:10.1088/2041-8205/750/2/L40
  * Kuncarayakti et al. (2020) Kuncarayakti, H., Folatelli, G., Maeda, K., et al. 2020, ApJ, 902, 139. doi:10.3847/1538-4357/abb4e7
  * Maeda et al. (2014) Maeda, K., Katsuda, S., Bamba, A., et al. 2014, ApJ, 785, 95. doi:10.1088/0004-637X/785/2/95
  * Maeda et al. (2015) Maeda, K., Hattori, T., Milisavljevic, D., et al. 2015, ApJ, 807, 35. doi:10.1088/0004-637X/807/1/35
  * Maund et al. (2004) Maund, J. R., Smartt, S. J., Kudritzki, R. P., et al. 2004, Nature, 427, 129. doi:10.1038/nature02161
  * Maund et al. (2011) Maund, J. R., Fraser, M., Ergon, M., et al. 2011, ApJ, 739, L37. doi:10.1088/2041-8205/739/2/L37
  * McMullin et al. (2007) McMullin, J. P., Waters, B., Schiebel, D., et al. 2007, Astronomical Data Analysis Software and Systems XVI, 376, 127
  * Moffat et al. (1986) Moffat, A. F. J., Lamontagne, R., Shara, M. M., et al. 1986, AJ, 91, 1392. doi:10.1086/114116
  * Montes et al. (2000) Montes, M. J., Weiler, K. W., Van Dyk, S. D., et al. 2000, ApJ, 532, 1124. doi:10.1086/308602
  * Nakar & Sari (2010) Nakar, E. & Sari, R. 2010, ApJ, 725, 904. doi:10.1088/0004-637X/725/1/904
  * Nakar & Piro (2014) Nakar, E. & Piro, A. L. 2014, ApJ, 788, 193. doi:10.1088/0004-637X/788/2/193
  * Nasonova et al. (2011) Nasonova, O. G., de Freitas Pacheco, J. A., & Karachentsev, I. D. 2011, A&A, 532, A104. doi:10.1051/0004-6361/201016004
  * Nayana et al. (2017) Nayana, A. J., Chandra, P., Roy, S., et al. 2017, MNRAS, 467, 155. doi:10.1093/mnras/stx044
  * Nayana et al. (2018) Nayana, A. J., Chandra, P., & Ray, A. K. 2018, ApJ, 863, 163. doi:10.3847/1538-4357/aad17a
  * Nayana & Chandra (2020) Nayana, A. J. & Chandra, P. 2020, MNRAS, 494, 84. doi:10.1093/mnras/staa700
  * Nayana & Chandra (2021) Nayana, A. J. & Chandra, P. 2021, ApJ, 912, L9. doi:10.3847/2041-8213/abed55
  * Pastorello et al. (2008) Pastorello, A., Kasliwal, M. M., Crockett, R. M., et al. 2008, MNRAS, 389, 955. doi:10.1111/j.1365-2966.2008.13618.x
  * Piro (2015) Piro, A. L. 2015, ApJ, 808, L51. doi:10.1088/2041-8205/808/2/L51
  * Piro et al. (2017) Piro, A. L., Muhleisen, M., Arcavi, I., et al. 2017, ApJ, 846, 94. doi:10.3847/1538-4357/aa8595
  * Podsiadlowski et al. (1992) Podsiadlowski, P., Joss, P. C., & Hsu, J. J. L. 1992, ApJ, 391, 246. doi:10.1086/171341
  * Rabinak & Waxman (2011) Rabinak, I. & Waxman, E. 2011, ApJ, 728, 63. doi:10.1088/0004-637X/728/1/63
  * Romero-Cañizales et al. (2014) Romero-Cañizales, C., Herrero-Illana, R., Pérez-Torres, M. A., et al. 2014, MNRAS, 440, 1067. doi:10.1093/mnras/stu430
  * Roming et al. (2009) Roming, P. W. A., Pritchard, T. A., Brown, P. J., et al. 2009, ApJ, 704, L118. doi:10.1088/0004-637X/704/2/L118
  * Ryder et al. (2004) Ryder, S. D., Sadler, E. M., Subrahmanyan, R., et al. 2004, MNRAS, 349, 1093. doi:10.1111/j.1365-2966.2004.07589.x
  * Sapir & Waxman (2017) Sapir, N. & Waxman, E. 2017, ApJ, 838, 130. doi:10.3847/1538-4357/aa64df
  * Sahu et al. (2013) Sahu, D. K., Anupama, G. C., & Chakradhari, N. K. 2013, MNRAS, 433, 2. doi:10.1093/mnras/stt647
  * Smith & Conti (2008) Smith, N. & Conti, P. S. 2008, ApJ, 679, 1467. doi:10.1086/586885
  * Smith (2014) Smith, N. 2014, ARA&A, 52, 487. doi:10.1146/annurev-astro-081913-040025
  * Smith & Tombleson (2015) Smith, N. & Tombleson, R. 2015, MNRAS, 447, 598. doi:10.1093/mnras/stu2430
  * Soderberg et al. (2006) Soderberg, A. M., Chevalier, R. A., Kulkarni, S. R., et al. 2006, ApJ, 651, 1005. doi:10.1086/507571
  * Soderberg et al. (2012) Soderberg, A. M., Margutti, R., Zauderer, B. A., et al. 2012, ApJ, 752, 78. doi:10.1088/0004-637X/752/2/78
  * Stockdale et al. (2003) Stockdale, C. J., Weiler, K. W., Van Dyk, S. D., et al. 2003, ApJ, 592, 900. doi:10.1086/375737
  * Stockdale et al. (2007) Stockdale, C. J., Williams, C. L., Weiler, K. W., et al. 2007, ApJ, 671, 689. doi:10.1086/522584
  * Swartz et al. (1993) Swartz, D. A., Clocchiatti, A., Benjamin, R., et al. 1993, Nature, 365, 232. doi:10.1038/365232a0
  * Tartaglia et al. (2017) Tartaglia, L., Fraser, M., Sand, D. J., et al. 2017, ApJ, 836, L12. doi:10.3847/2041-8213/aa5c7f
  * Taubenberger et al. (2011) Taubenberger, S., Navasardyan, H., Maurer, J. I., et al. 2011, MNRAS, 413, 2140. doi:10.1111/j.1365-2966.2011.18287.x
  * van der Hucht et al. (1988) van der Hucht, K. A., Hidayat, B., Admiranto, A. G., et al. 1988, A&A, 199, 217
  * Van Dyk et al. (2013) Van Dyk, S. D., Zheng, W., Clubb, K. I., et al. 2013, ApJ, 772, L32. doi:10.1088/2041-8205/772/2/L32
  * Van Dyk et al. (2013) Van Dyk, S. D., Zheng, W., Clubb, K. I., et al. 2013, ApJ, 772, L32. doi:10.1088/2041-8205/772/2/L32
  * Van Dyk et al. (2014) Van Dyk, S. D., Zheng, W., Fox, O. D., et al. 2014, AJ, 147, 37. doi:10.1088/0004-6256/147/2/37
  * van Moorsel et al. (1996) van Moorsel, G., Kemball, A., & Greisen, E. 1996, Astronomical Data Analysis Software and Systems V, 101, 37
  * Weiler et al. (1986) Weiler, K. W., Sramek, R. A., Panagia, N., et al. 1986, ApJ, 301, 790. doi:10.1086/163944
  * Weiler et al. (1992) Weiler, K. W., van Dyk, S. D., Pringle, J. E., et al. 1992, ApJ, 399, 672. doi:10.1086/171959
  * Weiler et al. (1998) Weiler, K. W., Van Dyk, S. D., Montes, M. J., et al. 1998, ApJ, 500, 51. doi:10.1086/305723
  * Weiler et al. (2002) Weiler, K. W., Panagia, N., Montes, M. J., et al. 2002, ARA&A, 40, 387. doi:10.1146/annurev.astro.40.060401.093744
  * Weiler et al. (2007) Weiler, K. W., Williams, C. L., Panagia, N., et al. 2007, ApJ, 671, 1959. doi:10.1086/523258
  * Wheeler et al. (1993) Wheeler, J. C., Clocchiatti, A., & Another, D. 1993, IAU Circ., 5756
  * Woosley et al. (1994) Woosley, S. E., Eastman, R. G., Weaver, T. A., et al. 1994, ApJ, 429, 300. doi:10.1086/174319
  * Yoon et al. (2010) Yoon, S.-C., Woosley, S. E., & Langer, N. 2010, ApJ, 725, 940. doi:10.1088/0004-637X/725/1/940
  * Yoon et al. (2017) Yoon, S.-C., Dessart, L., & Clocchiatti, A. 2017, ApJ, 840, 10. doi:10.3847/1538-4357/aa6afe


11institutetext:  CITEC, Bielefeld University, Bielefeld, Germany
11email<EMAIL_ADDRESS>22institutetext: Recommendy
UG, Bielefeld, Germany
22email<EMAIL_ADDRESS>
# Explainable Artificial Intelligence for Improved Modeling of
Processes††thanks: This research was supported by the research training group
“Dataninja” (Trustworthy AI for Seamless Problem Solving: Next Generation
Intelligence Joins Robust Data Analysis) funded by the German federal state of
North Rhine-Westphalia, and supported by the European Commission Horizon for
ICU4COVID project, and the VW-Foundation for the project IMPACT.

Riza Velioglu() 1122 Jan Philip Göpfert() 22 André Artelt() 11 Barbara
Hammer() 11

###### Abstract

In modern business processes, the amount of data collected has increased
substantially in recent years. Because this data can potentially yield
valuable insights, automated knowledge extraction based on process mining has
been proposed, among other techniques, to provide users with intuitive access
to the information contained therein. At present, the majority of technologies
aim to reconstruct explicit business process models. These are directly
interpretable but limited concerning the integration of diverse and real-
valued information sources. On the other hand, Machine Learning (ML) benefits
from the vast amount of data available and can deal with high-dimensional
sources, yet it has rarely been applied to being used in processes. In this
contribution, we evaluate the capability of modern Transformer architectures
as well as more classical ML technologies of modeling process regularities, as
can be quantitatively evaluated by their prediction capability. In addition,
we demonstrate the capability of attentional properties and feature relevance
determination by highlighting features that are crucial to the processes’
predictive abilities. We demonstrate the efficacy of our approach using five
benchmark datasets and show that the ML models are capable of predicting
critical outcomes and that the attention mechanisms or XAI components offer
new insights into the underlying processes.

###### Keywords:

machine learning process mining transformer XAI.

##  1 Introduction

Data is collected on anything, at any time, and in any location. A study by
IBM [6] found that in 2020, 40 trillion gigabytes (40 zettabytes) were
generated. The majority of data in the digital realm, however, is
unstructured, making it impossible for humans to process it in its entirety
and leaving businesses struggling to manage such massive amounts of data. As a
result, one of today’s major issues for businesses is extracting information
and value from data contained in information systems.

Process Mining (PM) [2] is a relatively new area of research that lies between
process modeling and analysis, machine learning, and data mining. It is an
umbrella term for a family of techniques to facilitate the analysis of
operational processes by extracting knowledge from event logs in symbolic
form, typically in the form of process models. An event log serves as the
input to PM algorithms, which is essentially a database of events where each
event has (1) a case id: a unique identifier for a particular process
instance, (2) an activity: description of the event that is occurring, and (3)
a timestamp: timestamp of the activity execution. Resources, expenses, and
other event-related attributes may also be integrated by some techniques. PM
techniques mostly address three tasks: (1) process discovery: transform the
event log into a process model which describes the underlying process
generating such data: common techniques include the alpha algorithm, heuristic
miner, or inductive miner [1], (2) conformance checking: analyse discrepancies
between an event log and an existing process model to detect deviations, for
quality checking or risk management [22], (3) process enhancement: improve the
existing model’s performance in terms of specific process performance metrics
[34].

The initial focus of PM was the analysis of historical event data to generate
a process model. These monitoring systems are reactive in nature, allowing
users to detect a violation only after it has occurred. In contrast,
Predictive Process Mining (PPM) or Predictive Process Monitoring [18] aims for
forward-looking forms of PM with predictive qualities. It is a field that
combines machine learning with process mining, aiming for specific tasks such
as predicting next activity, activity suffix, next timestamp, remaining time,
attributes, attribute suffix, and outcome of running cases. Most tasks can be
modeled as classification problems, except next timestamp and remaining time
prediction tasks, which are regression problems. Some approaches which have
been used in this realm are based on recurrent neural networks (RNNs) or Long-
short term memory (LSTMs) [11, 29, 23]. Alternatives are based on Autoencoders
[20] and Convolutional neural networks [24]. More recently, Transformer [32]
models have gained a lot attention due to their overwhelming success in
computer vision [10] and natural language processing [7]. Inspired by their
success, recent approaches proposed a Transformer-based model for process data
[5]. A more detailed overview on different models for the aforementioned tasks
is covered in the survey [27].

Previous work in PPM faces a few challenges: unclear training/test set splits
and pre-processing schemes, which leads to noticeably different findings and
challenges w.r.t reproducibility [33]. Since datasets often contain duplicates
which are not respected by the evaluation schemes, results are often
confounded. In this work, we investigate the capability of important classical
ML technologies as well as modern Transformers on a variety of benchmark sets
with clear train/test splits, duplicates respected, and different types of
representation for classical schemes. More importantly, we do not treat the
methods as opaque schemes, but rather rely on feature relevance determination
and attention mechanisms to highlight the most important factors for each
model. This is in line with previous approaches such as [13], which visualizes
the feature relevances of a gated graph neural network or [14] which enhances
next activity prediction using LSTMs with counterfactual explanations for
wrongly classified samples. Yet both approaches evaluate the behavior on a
single and specific dataset only. In our approach, we systematically
investigate relevant ML models for different approaches and representations
under the umbrella of explainability. In this way, we demonstrate how, in many
cases, the presence of specific events can easily fool the system into relying
on ‘trivial’ signals. When these are removed, ML models reveal not only high
quality behavior but also intriguing insights into the relevance of more
subtle signals or sequences.

##### Contributions

Our contribution is twofold: a systematic study of the behavior of diverse ML
models with train-test split respecting the peculiarities of datasets in PM
and appropriate pre-processing, which prevent troubles due to data leakage.
Second, we utilize eXplainable Artificial Intelligence (XAI) techniques to
compute and visualize the most crucial features. Our code is open source and
available at https://github.com/rizavelioglu/ml4prom.

The remainder of this work is structured as follows: In Section 2, we present
the peculiarities of the datasets and our methodology. In Section 3, we
present experiments and results. In Section 4, we conclude our contribution
and discuss potential directions for future work.

##  2 Methodology

In this section we describe the datasets and how they are pre-processed to be
utilized for binary classification, and highlight a few domain-specific
peculiarities which have to be taken into account. Then we introduce both
classical ML models as well as the state-of-the-art Transformer model used
subsequently. Lastly, we present the XAI techniques applied.

###  2.1 Data

We focus on five widely used datasets which come from a variety of domains
including loan applications, road traffic management, and healthcare. These
datasets are benchmarks from PM and are also accessible within the well-known
ProM tools [31], with the exception of healthcare, which contains personal
data. Each dataset has been labeled in order to assess explainability and put
them in the ML context. The variety of datasets enables us to check the
robustness of our methodology across a wide range of domains. To apply ML
classifiers to event data, we first transform the data into a format that
resembles a binary classification problem. We present five real-life event
logs, each with its own evaluation criterion for a prediction task (referred
to as the positive/negative subpart of the logs), and thereafter explain how
they are transformed for the task.

##### BPIC 2017_application+offer [8]:

a loan application process of a Dutch financial institute that covers
$31\,509$ cases. Each case represents a loan application and every application
has an outcome: positive if the offer is accepted and negative if rejected.

##### BPIC 2018_application [9]:

an agricultural grant-application process of the European Union that contains
$43\,809$ applications over a period of three years together with their
duration.

##### Traffic Fine Management [16]:

an event log of an information system managing road traffic fines by the
Italian government, with fines being payed or not

##### COVID [26]:

a dataset of COVID-19 patients hospitalized at an intensive care unit
consisting of 216 cases, of which 196 are complete cases (patients have been
released either dead or alive) and 20 ongoing cases (partial process traces)
that were being treated in the COVID unit at the time data was exported.

##### Hospital Billing [19]:

an event log of the billing of medical services that have been provided by a
regional hospital. Each trace in the event log keeps track of the actions
taken to bill a group of medical services. A randomly selected sample of
$100\,000$ process instances make up the event log.

###  2.2 Encoding

Table 1 shows the rules used to label the samples as well as important
statistics. One characteristic of typical PM datasets is that they can contain
a large number of duplicates, i.e. observations of the same process. For
example, Hospital dataset contains $41\,343$ and $58\,653$ traces in positive
and negative classes, respectively. Out of those traces only $306$ and $884$
are unique. In addition, the top-$10$ most frequent traces make up the $95.3$%
and $89.7$% of the whole traces. Hence such data are challenging for ML due to
a narrow variety. In addition, duplicates need to be accounted for in
evaluation to avoid information leakage.

Table 1: Datasets statistics and the rules used to generate binary classification task (L+/L- represent positive/negative logs, respectively). The statistics are the number of traces, the number of unique traces, the percentage of unique traces, and lastly the cumulative percentage of the top-10 most frequent traces. Dataset | Class | Classification Criteria | traces | uniq. | uniq./% | top-10%
---|---|---|---|---|---|---
BPIC17 |  L+ | without activity ‘A_incomplete’ | $16\,506$ | $529$ | $3.2$ | $83.5$
L- | with activity ‘A_incomplete’ | $15\,003$ | $2101$ | $14.0$ | $51.5$
Traffic |  L+ | end activity ‘Payment’ | $67\,201$ | $122$ | $0.2$ | $99.1$
L- | end activity not ‘Payment’ | $83\,169$ | $109$ | $0.1$ | $99.2$
COVID |  L+ | with activity ‘Discharge alive’ | $136$ | $33$ | $24.0$ | $74.3$
L- | with activity ‘Discharge dead’ | $60$ | $20$ | $33.0$ | $83.3$
BPIC18 |  L+ | duration less than 9 months | $27\,966$ | $1081$ | $3.9$ | $57.7$
L- | duration more than 9 months | $15\,843$ | $477$ | $3.0$ | $82.1$
Hospital |  L+ | duration less than 3 months | $41\,343$ | $306$ | $0.7$ | $95.3$
L- | duration more than 3 months | $58\,653$ | $884$ | $1.5$ | $89.7$

Data consist of sequences of symbolic events of different length. For each of
these datasets we compute $n$-grams for $n\in\\{1,2,3\\}$ to encode traces in
vectorial form for classic ML models. Unigrams (n=1) encode occurrence of
events, bigrams (n=2) encode two subsequent events and trigrams (n=3) encode
three subsequent events. Transformer models can directly deal with sequences.
The network learns a vector embedding for each event within a trace. This has
the advantage of avoiding problems caused by the high-dimensionality of one-
hot encoding method, for example.

###  2.3 Models

Using these representations, we train a variety of models, including Logistic
Regression (LR), Decision Tree (DT), Random Forest (RF), Gradient Boosting
(GB), and Transformer models. We only present the LR, DT and Transformer
models because the findings do not vary significantly. We selected LR and DT
due to their widespread use and overall high performance for classification
problems [12], and we selected Transformer models due to their propensity for
learning complex patterns. The type of encoding used in this case limits the
amount of information that is accessible because only a portion of the
sequential structure is represented. Conversely, transformer models can easily
handle sequential data. We use models as proposed in the works [32, 5].

###  2.4 XAI

Since our primary concern is not the classification performance of our trained
models, but rather whether ML models can capture underlying regularities, we
employ XAI techniques to gain insight into which relationships between events
in the data are used by the models [21]. Unlike many explanation approaches
such as LIME [28], LRP [3], we are interested in global explanations rather
than explanations of single decisions. This is because our goal is not to
comprehend the ML model, but rather to gain understanding of the key
components of the PM as a whole, which can enable users to enhance processes.
As an example, one may anticipate altered process behavior if they come across
a specific activity, such as ‘Discharge dead’ in COVID dataset, which was
identified as crucial for the ML model. As a result, we strive for global
explanations and, more specifically, we employ a number of well-established
feature selection methods for both linear and non-linear settings, including
LR with LASSO as a linear model with strong mathematical guarantees [30], DT
as a non-linear model with efficient Mean Decrease in Impurity (MDI) and
permutation importance [17] that takes into account non-linear relations
between features to determine the relevances, and Transformer model equipped
with an attention mechanism that highlights complex relationships. Attention
mechanisms are technically local explanations. We utilize them to see if
highly flexible non-linear explanations provide more complex relations than
global feature importance measures.

As a first result of the analysis, feature selection methods immediately
reveal a potential source of information leakage: Leaving the datasets as they
are, the algorithms rely on attributes that directly encode the class label
but provide no additional insight into the process. Table 2 shows the detected
features that are removed from the logs to avoid such trivial outcomes. All of
the listed events exist only in their respective class, e.g. ‘A_Incomplete’ is
present only in L- but not in L+, except ‘Payment’ where the event is present
in both classes. Therefore, we removed the event only if it is the last event
in a trace from L+, as that is the source of leakage (see Table 1).

Table 2: Datasets with features that leak label information. Dataset | Class | Biased feature
---|---|---
BPIC17 |  L- | ‘A_Incomplete’
Traffic |  L- | ‘Send for Credit Collection’
L+ | ‘Payment’
COVID |  L- | ‘Discharge dead’
L+ | ‘Discharge alive’

##  3 Experimental Results

In this section we present the data pre-processing pipeline used to transform
all five datasets for training. Then, we present the evaluation metric used as
well as the model design and training on each dataset. Finally, we present the
scores of different approaches on the binary classification task and the
resulting feature relevances. To save space, we limit ourselves to the BPIC17
dataset; results for the other datasets can be found in the GitHub repository.

###  3.1 Data Split and Pre-Processing

For simplicity, we only consider event names to encode traces and no other
attributes, e.g. timestamp, or resource of an event. After computing unique
traces we randomly sample from them to construct train/test sets with a ratio
of $70$%/$30$%, respectively. As the datasets are highly imbalanced, we sample
class-wise from data–preserving the proportion of classes both in train and
test sets. To account for frequency of traces, we keep the duplicate traces in
train set while removing the ones in test set.

We apply minimal pre-processing to data at hand. First, we remove the biased
features from traces. Then we add <start> and <end> tokens to the input to
explicitly define the beginning and the end of traces. To have a fixed-length
input, the traces that are shorter than the longest trace in an event log are
padded with the <pad> token. We then build a token dictionary consisting of
unique event names in an event log and the aforementioned special tokens.
Finally, the tokens in the input are replaced by their unique integer values
stored in the dictionary.

###  3.2 Evaluation Metric

Because the datasets are highly imbalanced, we measure performance using the
area under the receiver operating characteristic curve (AUROC) [4]. ROC
analysis does not favor models that outperform the majority class while under-
performing the minority class, which is desirable when working with imbalanced
data [20, p. 27]. The area under the curve of a binary classifier is the
probability that a classifier would rank a randomly chosen positive instance
higher than a randomly chosen negative one, which is given by the following
equation:

$\texttt{AUROC}=\int_{x=\infty}^{-\infty}\texttt{TPR}(T)\texttt{FPR}^{{}^{\prime}}(T)dT$
(1)

AUROC ranges from $0$ to $1$, with an uninformative classifier (random
classifier) producing a result of $0.5$.

###  3.3 Model Design and Training

For ML models we utilize scikit-learn [25] and initialize the models with
default parameters. For logistic regression, we employ regularization by a L1
penalty term(penalty hyper-parameter) with regularization strength (C hyper-
parameter). We use repeated stratified k-fold cross-validation (CV) to
evaluate and train a model across a number of iterations. The number of
repeats is $50$, and the number of splits $\textit{k}=5$, yielding $250$
models being trained. The reported metric is then the average of all scores
computed.

For Transformer model, following [5], we chose the embedding dimension as
$36$, i.e. each trace is represented by a point in $36$-dimensional space.
Since Transformer disregards the positional information of events in traces,
we add positional encoding to token embedding which have the same dimensions.
During training, the model learns to pay attention to input embedding as well
as positional encoding in an end-to-end fashion. The embedding outputs are
then fed to a multi-head attention block with $h=6$ heads. On the final layer
of the attention block, we aggregate features using a global average pooling
followed by a dropout at a rate of $0.1$. Then, we employ a dense layer with
ReLU activation of $64$ hidden units and a dropout at a rate of $0.1$.
Finally, we use a dense layer with sigmoid activation that outputs a value
between $0$ and $1$, which is interpreted as the “probability” of a trace
belonging to positive (desirable) class. We train the model for $50$ epochs
with ADAM optimizer [15], a learning rate of $1\text{\times}{10}^{-3}$, and
batch size of $16$.

###  3.4 Results

####  3.4.1 Predictive Accuracy

We report the experimental results for the binary classification task. Table 3
reports the AUROC scores of models on BPIC17, Traffic, and COVID datasets,
where there are some features leaking the label information. Here the perfect
score is achieved as the models correctly discover the biased features, as
expected. However, the Transformer model and the ML models with $1$-gram on
Traffic dataset do not achieve the perfect score on the test set. This is due
to the fact that the biased feature–‘Payment’–appears in both of the classes.
Therefore, $1$-gram encoding is not capable of leaking the information. The
Transformer model also fails to discover the bias. On the other hand, models
achieve the expected results with the $2$-gram and $3$-gram encoding because
those encoding types explicitly define the biased feature: as <end> token
added to traces, the feature (Payment, <end>) in $2$-gram and (Payment, <end>,
<PAD>) in $3$-gram would leak the label information (positive class if last
event is ‘Payment’, negative class otherwise).

Table 3: AUROC scores on datasets with and without biased features using
different encoding methods and models. The values inside the parenthesis
represent the score on the hold-out test set, whereas others represent the
score on the training set.

|  | with biased features | without biased features
---|---|---|---
Dataset | Encoding | LR | DT | Transformer | LR | DT | Transformer
BPIC17 | integer | - | - | 100.(100.) | - | - | 97.3(97.3)
1-gram | 100.(100.) | 100.(100.) | - | 89.0(81.3) | 89.5(76.0) | -
2-gram | 100.(100.) | 100.(99.9) | - | 97.4(97.9) | 97.3(93.5) | -
3-gram | 100.(100.) | 99.9(99.8) | - | 97.4(97.8) | 97.3(91.7) | -
Traffic | integer | - | - | 100.(98.4) | - | - | 63.9(53.5)
1-gram | 100.(93.3) | 100.(92.2) | - | 61.5(36.8) | 61.5(52.7) | -
2-gram | 100.(100.) | 100.(100.) | - | 63.8(45.5) | 63.9(49.8) | -
3-gram | 100.(100.) | 100.(100.) | - | 63.8(48.1) | 63.9(52.3) | -
COVID | integer | - | - | 90.9(100.) | - | - | 74.8(94.2)
1-gram | 100.(100.) | 100.(100.) | - | 67.3(75.0) | 69.6(68.3) | -
2-gram | 100.(100.) | 99.9(100.) | - | 89.3(61.7) | 85.4(85.0) | -
3-gram | 99.9(96.7) | 98.1(66.7) | - | 90.4(48.3) | 85.5(76.7) | -
BPIC18 | integer | - | - | - | - | - | 98.6(81.2)
1-gram | - | - | - | 97.5(77.9) | 97.5(78.5) | -
2-gram | - | - | - | 98.4(87.3) | 98.2(84.5) | -
3-gram | - | - | - | 98.4(87.3) | 98.2(81.5) | -
Hospital | integer | - | - | - | - | - | 92.4(78.7)
1-gram | - | - | - | 91.8(73.9) | 92.3(56.2) | -
2-gram | - | - | - | 92.5(70.5) | 92.7(52.5) | -
3-gram | - | - | - | 92.6(66.0) | 92.6(48.5) | -

Table 3 also reports the AUROC scores of models on all of the five real-life
event logs, where the biased features are removed from BPIC17, Traffic, and
COVID datasets. For BPIC17 and COVID datasets the scores worsen but the
results are still promising, hinting that the processes, sequences of events,
maintain valuable information for the task even though the distinct/biased
features are removed. The different amount of information in different
encoding is mirrored by the results: $1$-gram achieves the worst results
compared to $2$-gram and $3$-gram, as it only incorporates single events,
whereas $2$-gram and $3$-gram incorporates event pairs. None of those methods
integrate order information, while Transformer model learns this information
during training. We observe that Transformer model successfully captures the
relations between events as well as the order information and outperforms
other models in Traffic and COVID datasets, whereas in BPIC17 it receives a
comparable result. On the other hand, we observe that the scores for BPIC17
and COVID datasets do not fluctuate as much as it does for Traffic dataset
when compared to biased scores represented in Table 3. This is due to the fact
that after removing the biased feature in Traffic dataset, some traces in both
classes become identical.

####  3.4.2 Feature Relevances

Figure 1: Relevance of features as considered by the logistic regression
model. A high positive/negative value indicates a considerable contribution
towards predicting the positive/negative label, i.e. desirable/undesirable
event trace.

We present the relevances of features taken into account by our trained
logistic regression model in Figure 1. We find that ‘A_Validating’ contributes
the most towards predicting an undesirable trace, whereas ‘O_Refused’
contributes most to predict a desirable trace. Other events have significantly
lower impact on the predictions. Interestingly, the special tokens, i.e.
<pad>, <start>, and <end> affect the predictions when they should not, despite
the fact that their influence is negligible. In addition, some features have
no effect on the predictions, e.g. ‘A_Denied’, ‘A_Cancelled’.

Figure 2: Relevance of features as considered by decision tree model. Left:
Relevances according to MDI: A high value indicates a considerable
contribution towards deciding between the positive and negative label. Right:
Relevances according to the permutation importance calculated on the test set,
i.e. how much shuffling a given feature negatively impacts the performance of
the decision tree.

In Figure 2 we show the feature relevances that our trained decision tree
model considers. Based on the MDI value (left figure) the most relevant
feature is ‘A_Validating’, followed by ‘O_Refused’, ‘A_Submitted’, and
‘O_Cancelled’ with significantly lower impact. Based on the permutation
feature importance (right figure) the most relevant features are
‘A_Validating’, and ‘O_Refused’, where the rest of the features have no effect
on the prediction performance, which aligns well with the relevances of LR
model.

Figure 3:  Left: Relevance of features as considered by transformer model for
one randomly selected event trace. The thickness of a line connecting two
features indicates the intensity of attention in between, whereas the color
represents one of the six different attention heads. Right: Normalized
attention scores which are averaged among all attention heads and all traces
in test set.

Figure 3 visualizes the attention scores of six attention heads for a given
trace, as well as normalized attention scores over all attention heads and
test samples. We observe that, regardless of where it appears in the trace,
the events mostly attend to ‘A_Validating’ event for the aforementioned trace.
In addition, some events focus on ‘O_Sent’ in some heads, which differs from
the other models. The normalized attention scores, however, demonstrates that
not all traces exhibit this behavior. The plot also highlights features whose
importances overlap with other models’ results. In summary, all models agree
on the most relevant features.

##  4 Conclusion

We have demonstrated how to prepare process data in such a way that it can be
used to train classic and modern – state-of-the-art – ML classifiers. All our
trained models exhibit high classification performance, i.e. they are capable
of learning the underlying regularity of the observed processes, whereby
Transformers benefit from the fact that the full sequence information is
available, unlike e.g. 1-gram representations. XAI technologies prevent
pitfalls such as information leakage by explicit encoding of the predicted
event, and reveal insights into the relevance of events or sequences of
events, respectively. These insights enable a further exploration of crucial
aspects of the processes, which is useful e.g. for the improvement or
correction of undesired process outcomes.

Future research could investigate the effects of various model architectures
and encoding schemes on the outcomes of feature relevances. Another potential
research direction might be to study the impact of adding further event
attributes on the learnt representations.

## References

  * [1] van der Aalst, W., Weijters, T., Maruster, L.: Workflow mining: Discovering process models from event logs. IEEE Trans Knowl Data Eng (2004)
  * [2] van der Aalst, W.: Process Mining: Data Science in Action. Springer (2016)
  * [3] Binder, A., Montavon, G., Bach, S., Müller, K., Samek, W.: Layer-wise relevance propagation for neural networks with local renormalization layers. CoRR abs/1604.00825 (2016)
  * [4] Bradley, A.P.: The use of the area under the roc curve in the evaluation of machine learning algorithms. Pattern recognition 30(7), 1145–1159 (1997)
  * [5] Bukhsh, Z.A., Saeed, A., Dijkman, R.M.: Processtransformer: Predictive business process monitoring with transformer network. CoRR (2021)
  * [6] Dailey, D.: Netezza and ibm cloud pak for data: A knockout combo for tough data. https://ibm.co/3xvK4MG, [accessed 17-June-2022]
  * [7] Devlin, J., Chang, M.W., Lee, K., Toutanova, K.: Bert: Pre-training of deep bidirectional transformers for language understanding. arXiv:1810.04805 (2018)
  * [8] van Dongen, B.: BPI Chall. 2017 (2 2017). https://doi.org/jcmn
  * [9] van Dongen, B., Borchert, F.: BPI Chall. 2018 (3 2018). https://doi.org/jcmm
  * [10] Dosovitskiy, A., Beyer, L., Kolesnikov, A., Weissenborn, D., Zhai, X., Unterthiner, T., Dehghani, M., Minderer, M., Heigold, G., Gelly, S., et al.: An image is worth 16x16 words: Transformers for image recognition at scale. arXiv:2010.11929 (2020)
  * [11] Evermann, J., Rehse, J.R., Fettke, P.: Predicting process behaviour using deep learning. Decis Support Syst 100, 129–140 (2017)
  * [12] Fernández, M., Cernadas, E., Barro, S., Amorim, D.: Do we need hundreds of classifiers to solve real world classification problems? JMLR 15(1) (2014)
  * [13] Harl, M., Weinzierl, S., Stierle, M., Matzner, M.: Explainable predictive business process monitoring using gated graph neural networks. J Decis Syst 29(sup1), 312–327 (2020)
  * [14] Hsieh, C., Moreira, C., Ouyang, C.: Dice4el: interpreting process predictions using a milestone-aware counterfactual approach. In: ICPM. pp. 88–95. IEEE (2021)
  * [15] Kingma, D.P., Ba, J.: Adam: A method for stochastic optimization (2014)
  * [16] de Leoni, M.M., Mannhardt, F.: Road Traffic Fine Management Process (2 2015). https://doi.org/jcmk
  * [17] Louppe, G.: Understanding random forests: From theory to practice (2014). https://doi.org/jcms
  * [18] Maggi, F.M., Di Francescomarino, C., Dumas, M., Ghidini, C.: Predictive monitoring of business processes (2013)
  * [19] Mannhardt, F.: Hospital Billing - Event Log (8 2017). https://doi.org/gm85w4
  * [20] Mehdiyev, N., Evermann, J., Fettke, P.: A novel business process prediction model using a deep learning method. Business&inf. sys. eng. 62(2), 143–157 (2020)
  * [21] Molnar, C.: Interpretable Machine Learning. Lulu. com, 2 edn. (2022)
  * [22] Munoz-Gama, J.: Conformance checking and diagnosis in process mining, vol. 270. Springer (2016)
  * [23] Nguyen, A., Chatterjee, S., Weinzierl, S., Schwinn, L., Matzner, M., Eskofier, B.: Time matters: time-aware lstms for predictive business process monitoring. In: ICPM. pp. 112–123. Springer (2020)
  * [24] Pasquadibisceglie, V., Appice, A., Castellano, G., Malerba, D.: Using convolutional neural networks for predictive process analytics. In: ICPM. IEEE (2019)
  * [25] Pedregosa, F., Varoquaux, G., Gramfort, A., Michel, V., Thirion, B., Grisel, O., Blondel, M., Prettenhofer, P., Weiss, R., Dubourg, V., et al.: Scikit-learn: Machine learning in python. JMLR 12, 2825–2830 (2011)
  * [26] Pegoraro, M., Narayana, M.B.S., Benevento, E., van der Aalst, W.M., Martin, L., Marx, G.: Analyzing medical data with process mining: A covid-19 case study. In: Int. Conf. on Business Inf. Sys. pp. 39–44. Springer (2022)
  * [27] Rama-Maneiro, E., Vidal, J., Lama, M.: Deep learning for predictive business process monitoring: Review and benchmark. IEEE Trans. Serv. Comput. (2021)
  * [28] Ribeiro, M.T., Singh, S., Guestrin, C.: ”why should I trust you?”: Explaining the predictions of any classifier. In: ACM SIGKDD. pp. 1135–1144 (2016)
  * [29] Tax, N., Verenich, I., Rosa, M.L., Dumas, M.: Predictive business process monitoring with lstm neural networks. In: CAiSE. pp. 477–492. Springer (2017)
  * [30] Tibshirani, R.: Regression shrinkage and selection via the lasso. J. R. Stat. Soc. 58(1), 267–288 (1996). https://doi.org/gfn45m
  * [31] Van Dongen, B.F., de Medeiros, A.K.A., Verbeek, H., Weijters, A., van Der Aalst, W.M.: The prom framework: A new era in process mining tool support. In: Int. Conf. on App. and Theory of Petri Nets. pp. 444–454. Springer (2005)
  * [32] Vaswani, A., Shazeer, N., Parmar, N., Uszkoreit, J., Jones, L., Gomez, A.N., Kaiser, L., Polosukhin, I.: Attention is all you need (2017)
  * [33] Weytjens, H., De Weerdt, J.: Creating unbiased public benchmark datasets with data leakage prevention for predictive process monitoring (2021). https://doi.org/jcmp
  * [34] Yasmin, F., Bukhsh, F., De Alencar Silva, P.: Process enhancement in process mining: A literature review. CEUR (2018)