Radio Analysis of SN 2004C Reveals an Unusual CSM Density Profile as a Harbinger of Core Collapse

We present extensive multi-frequency VLA and VLBA observations of the radio-bright supernova (SN) IIb SN 2004C that span $\sim(40-2793)$ days post-explosion. We interpret the temporal evolution of the radio spectral energy distribution (SED) in the context of synchrotron self-absorbed (SSA) emission from the explosion's forward shock as it expands in the circumstellar medium (CSM) previously sculpted by the mass-loss history of the stellar progenitor. VLBA observations and modeling of the VLA data point to a blastwave with average velocity $\sim0.06c$ that carries an energy of $\sim 10^{49}$ erg. Our modeling further reveals a flat CSM density profile $\rho_{\rm{CSM}} \propto R^{-0.03 \pm0.22}$ up to a break radius $R_{br} \approx (1.96 \pm 0.10) \times 10^{16}$ cm, with a steep density gradient following $\rho_{\rm{CSM}} \propto R^{-2.3 \pm 0.5}$ at larger radii. We infer that the flat part of the density profile corresponds to a CSM shell with mass $\sim0.021 M_{\odot}$, and that the progenitor's effective mass-loss rate varied with time over the range $(50-500) \times 10^{-5} M_{\odot} \rm{yr}^{-1}$ for an adopted wind velocity $v_w =1000$ km $s^{-1}$ and shock microphysical parameters ${\epsilon}_e = 0.1, {\epsilon}_B = 0.01$. These results add to the mounting observational evidence for departures from the traditional single-wind mass-loss scenarios in evolved, massive stars in the centuries leading up to core collapse. Potentially viable scenarios include mass loss powered by gravity waves and/or interaction with a binary companion.


INTRODUCTION
All stars lose mass to their environments throughout their lifetimes, critically affecting both the star's evolution and fate in addition to enriching the Universe with heavy elements. If the star dies in a supernova (SN) explosion, shock interaction with recent mass lost to the circumstellar medium shapes the appearance of the SN display.
In recent years, observations of supernovae (SNe) have uncovered evidence that evolved massive stars live complex lives with a rich and diverse mass-loss history, particularly in the final moments before their deaths. Large optical datasets that sampled the pre-explosion phase of SNe (such as the targets observed by the Palomar Transient Factory -PTF-and the Zwicky Transient Facility -ZTF-; Law et al. 2009;Bellm et al. 2019;Graham et al. 2019) have unveiled evidence of pre-explosion eruptions arXiv:2203.07388v1 [astro-ph.HE] 14 Mar 2022 in > 50% of progenitors of Type-IIn SNe , later updated by Strotjohann et al. 2021), which are SNe that explode into dense H-rich CSM and show narrow hydrogen lines in their spectra as a consequence (Filippenko 1997). Among Type-IIn SNe, the precursor event of SN 2010mc (Ofek et al. 2013) and the long history of pre-SN eruptions of SN 2009ip (Smith et al. 2010;Foley et al. 2011) were among the first events inferring the presence of a dense CSM shell extending ∼ 5×10 14 − ∼ 4 × 10 16 cm from their explosion sites. The existence of these shells points to a mass-loss mechanism different from the traditional line-driven winds, and hence exposes gaps in our understanding of massive stellar evolution (Fraser et al. 2013;Mauerhan et al. 2013Mauerhan et al. , 2014Pastorello et al. 2013;Prieto et al. 2013;Margutti et al. 2014;Smith & Arnett 2014). Interestingly, evidence of large mass-loss events (0.01− several M lost to the environment) occurring in the centuries to years before death extends to the H-rich progenitors of superluminous SNe (SLSNe-II), a class assigned to SNe orders of magnitude more luminous than normal ones that also show H in early spectra. SN 2008es is one such example, which interacted with 2 − 3M material 0.5-1.6 years prior to explosion (Bhirombhakdi et al. 2019).
Radio observations of SNe played a key role in our understanding of "unusual" CSM density profiles around SNe (e.g., Chevalier & Fransson 2017 for a recent review). In a SN explosion, while optical observations sample the slowly expanding ejecta (v ≤ 10 4 km s −1 ) radio observations probe the fastest ejecta (v ≥ 0.1c). Radio synchrotron emission originates from the interaction of the fastest SN ejecta with the local CSM deposited by the progenitor star prior to explosion. By monitoring the temporal evolution of the spectral peak frequency and flux density, one may directly constrain physical properties of the system, such as the mass-loss rateṀ of the progenitor star. Therefore, radio observations allow us to map the density profile of the CSM around the explosion site. Previous works in the literature have built a foundation of using this approach to measure pre-explosion mass loss of H-stripped progenitors at radio frequencies as the blastwave interacts with the surrounding CSM (see e.g., SNe 2001em, 2003bg, 2003gk, 2007bg, 2008ax, and PTF11qcj;Soderberg et al. 2006a;Bietenholz et al. 2014;Salas et al. 2013;Schinzel et al. 2009;Roming et al. 2009;Corsi et al. 2014). Using the standard techniques of radio analysis and radio data modeling (as reviewed in Weiler et al. 2002;Chevalier & Fransson 2017), we have found that SN 2004C is one of such H-poor explosions with a shockwave that encounters a dense CSM with a profile that is not consistent with the traditional single stellar wind models.
Observations across the electromagnetic spectrum of SNe thus point to the presence of strongly timedependent mass loss predicating core collapse in a diversity of explosion types, ranging from stripped-envelope SNe Ib/c to SNe II and extending to SLSNe. The fact that these mass-loss ejections occur across different types of SNe implies an origin of the outbursts that is independent of progenitor size. The exact nature of the underlying physical process causing this time-dependent mass loss remains the subject of debate, as well as how these processes inform the evolutionary path of the progenitor (e.g., Smith 2014).
In this work we analyze the radio data of SN 2004C, a SN Type IIb that exploded December 15, 2003 (MJD 52988, see §2.1) located in the galaxy NGC 3683, 35.9 Mpc away Tully et al. 2009Tully et al. , 2013. Our radio dataset spans 40 − 2793 days post explosion and represents one of the most extensive radio datasets of a SN thus far, which we model with synchrotron self-absorption. Furthermore, Very Long Baseline Array (VLBA) imaging of SN 2004C at 1931 days post explosion offers a direct measurement of the blastwave radius, making SN 2004C among the ∼ 10 SN with published VLBA constraints to date. With the use of our generalized equations from Chevalier (1998) (hereafter C98), we find the explosion site of SN 2004C contains a dense shell of CSM material and we speculate on its physical origin.
Our paper is organized as follows. The VLA and VLBA data of SN 2004C are presented in §2, and our generalization of the equations following the formalism of C98 is in §3. Our modeling of the data and analysis of the physical parameters of the system are presented in §4. Our comparison to plausible mass-loss mechanisms are presented in §5, and we summarize our results in §6.
2. OBSERVATIONS 2.1. Explosion date and SN typing SN 2004C was discovered January 12.5, 2004 UT in near-infrared K'-band (2.1µm) imaging with relatively poor constraints on the time of explosion . The most recent prior observation of that region was taken May 18, 2003 and published photometry are insufficient to robustly constrain the time of maximum light and potentially extrapolate an explosion date. The only known attempt was made by Meikle et al. (2004), who estimated that SN 2004C was 3 weeks past maximum light as of January 16, 2004 by comparing its g − r and r − i colors to those of Type Ic supernovae. That estimate, however, is problematic given that years later Shivvers et al. (2017) reported that the spectral classification of SN 2004C was in fact Type IIb, and not Type Ic as originally reported (Matheson et al. 2004).
To better estimate the date of explosion, we supplemented these reports and the updated spectroscopic classification with our own investigation of archival optical spectra of SN 2004C retrieved from WISeREP (Yaron & Gal-Yam 2012). Spectra spanned January 15, 2004 to March 17, and were inspected for their likeness to other Type IIb spectra tagged with time since explosion and/or maximum light. Our comparisons were aided with use of SNID (Blondin & Tonry 2007); and we heavily weighted a high quality Lick 3m spectrum of SN 2004C obtained January 17, 2004 and originally published in Shivvers et al. (2017) in our manual comparisons. Among all Type IIb supernova spectra inspected, those of SN 1993J and SN 2000H displayed the closest similarity in spectral evolution with the limited available epochs of SN 2004C. Comparisons converged on the time of maximum light of SN 2004C to be around January 2, 2004. Given an average Type IIb V-band rise-topeak time of approximately 18 days from explosion (see, e.g., Milisavljevic et al. 2013), we estimate the explosion date to be December 15, 2003 with an uncertainty of ∼ 2 weeks dominated by ambiguity between fits with other epochs. In the following we adopt MJD 52988 (corresponding to December 15, 2003) as the explosion date. This uncertainty has no impact on our major conclusions.

Karl G. Jansky Very Large Array
We observed SN 2004C with the NSF's Karl G. Jansky Very Large Array (VLA) beginning on January 23.65, 2004 UT, δt ≈ 30 days after the initial discovery under programs AK0575 and AS0796 (PI Soderberg). We found a radio source at α = 11 h 27 m 29.72 s , δ = +56 • 52 48.2 , coincident with the optical SN po-sition with a flux density of F ν = 1.6 ± 0.05 mJy at a frequency of ν = 8.5 GHz. We continued to monitor this source for the next δt ≈ 2800 days.
Radio data were taken in the 4.9, 8.5, 15.0, and 22.0 GHz bands, and we give the measured flux densities in Table 3. Observations were taken in continuum observing mode with 2 × 50 MHz bandwidth, except that in August 2011, which was taken with the upgraded EVLA system, with an increased 2 GHz bandwidth. We used observations of J1127+568 to calibrate the time dependent complex gains of the instrument and 3C286 to calibrate the absolute flux scale and the bandpass response for all observations. Data were reduced using the Astronomical Image Processing System (AIPS) by fitting a elliptical Gaussian model to the radio SN in each observation to measure the integrated flux density. We note that that SN 2004C has a peak luminosity density at 100 days of ∼ 10 28 erg s −1 Hz −1 and is at the high end of the radio luminosity function for SNe (see, e.g., Bietenholz et al. 2021). A subsection of data points showcasing the SN evolution is shown in Figure 1 and the full dataset is plotted in Appendix Figure 10. In Figure 2, we show the lightcurves for SN 2004C at 4.9, 8.5, 15.0, and 22.0 GHz bands. Temporally, the radio flux density at 8.5, 15.0, and 22.0 GHz rises slowly over the first 100 days and then steadily declines. At 4.9 GHz, the radio flux density increases for approximately the first 200 days, peaking at 10 mJy, and then plateaus for several hundred more days until beginning to decline. We see evidence for a low frequency turnover in our spectral energy distribution, which is readily apparent in Figure 10 of the Appendix, as well as in the raw radio lightcurves in Figure 2.

Very Long Baseline Array
We observed SN 2004C with the NRAO Very Long Baseline Array (VLBA) as part of program BS192 (PI Soderberg) on March 29, 2009 at 8.4 GHz. The observations were performed with eight frequency bands of 8 MHz bandwidth each in dual circular polarization, resulting in a total data rate of 512 Mbps. ICRF J112813.3+592514, located at α 2000 = 11 h 28 m 13.3406 s and δ 2000 = +59 • 25 14.798 (Fey et al. 2004) was used for phase-referencing at both frequencies. We used cycle time of ∼4 min on SN 2004C and 1 min on ICRF J112813.3+592514, and the total length of the observations was 8 hours.
The data were correlated at the VLBA Array Operations Center in Socorro, New Mexico and calibrated using AIPS and ParselTongue (Kettenis et al. 2006  Note that for SN 2004C the evolution of the peak flux density and peak frequency does not follow a power law with time and our modeling is designed to adapt to this non-standard evolution. Solid lines: best-fitting model. Solid points: radio observations of SN 2004C. Hollow points: intercepts of the asymptotic power-law segments, which mark the location of ν brk and F brk . These values are used to estimate the physical properties of the SN outflow and its environment.  2004C data at 4.9, 8.5, 15, and 22 GHz (grey points). Overlaid is our model fit to the data; shaded regions represent the spread of the 1σ uncertainty on the broken power law fit parameters (Table 1). plitude calibration used system temperature measure-ments and standard gain curves. We performed a "man-ual phase-calibration" using the data from one scan of J112813.3+592514 to remove instrumental phase offsets among the frequency bands. We then fringe-fitted the data from ICRF J112813.3+592514. Finally, the calibration was transferred to SN 2004C. The data were imaged in AIPS using robust weighting (with AIPS ROBUST=0). We clearly detect the source with an integrated flux density of F ν = 0.74 ± 0.13 mJy at coordinates α = 11 h 27 m 29.728309 s ± 0.000012 s , δ = 56 • 52 48.12937 ± 0.00014 . The beam size of our observation is ≈ 3 mas, and the source is clearly resolved in our data with a size of 0.707 +0.120 −0.178 mas (3.8 +0.6 −1.0 ×10 17 cm at 35.9 Mpc). Our flux density of F ν = 0.74 ± 0.13 mJy is consistent with our VLA flux density measurement taken on April 5.22, 2009 UT (7 days after the VLBA measurement) of 0.72 ± 0.07 mJy (Table 3).

RADIO MODELING
3.1. Radio Synchrotron Emission in SNe: emitting radius R, post-shock magnetic field B, and energy U In a SN explosion, the propagation of the fastest ejecta into the environment creates a double-shock structure: one shock propagates outward (the forward shock, or FS) and the other propagates inward in mass coordinates into the SN ejecta (the reverse shock, or RS). The two shocked regions of CSM and SN ejecta are separated by a contact discontinuity.
Radio emission in young non-relativistic SNe is dominated by synchrotron radiation from electrons accelerated at the FS (e.g., Chevalier 1998;Weiler et al. 2002;Chevalier & Fransson 2017). Relativistic electrons gyrate in a magnetic field B with pitch angle θ in the shocked CSM and radiate photons as non-thermal synchrotron emission (Rybicki & Lightman 1986). The acceleration of these electrons creates a power law distribution in electron Lorentz factor γ e , where dN (γ e )/dγ e = K 0 γ −p e for γ e ≥ γ m , where γ m is the minimum Lorentz factor. Here, N is the number of electrons per unit volume, K 0 is the normalization constant, and p is the power-law index of the electron distribution. 1 Therefore, the typical electron Lorentz factor of the powerlaw distribution is γ e ≈ γ m . Radio observations of SNe typically indicate p ≈ 3 (e.g., Weiler et al. 2002;Berger et al. 2002;Soderberg et al. 2006aSoderberg et al. , 2004Soderberg et al. , 2006cSoderberg et al. ,b, 2007Soderberg et al. , 2010bSoderberg et al. ,a, 2012Chevalier & Fransson 2006;Corsi et al. 2014Corsi et al. , 2016Kamble et al. 2016, for observational results, and a theoretical basis is given in Diesing & Caprioli 1 We note that the power-law variable p is sometimes expressed as γ in the radio SN literature. However, we wish to avoid confusion with the electron Lorentz factor, and adopt p instead. 2021). Diversity among the values of p in SN explosions might be a manifestation of different physical properties in the shocks where particle acceleration happens. The resulting radio flux density spectrum, F ν , is a series of broken power laws with spectral breaks that occur at frequencies and flux densities set by the physical properties of the system (e.g., B, electron density n e , etc.). Notable break frequencies are the synchrotron cooling frequency, ν c , the synchrotron characteristic frequency, ν m , and the synchrotron self-absorption frequency, ν sa . The slopes of F ν power-law segments depend on the order of ν c , ν m , and ν sa (e.g., Fig. 1 of Granot & Sari 2002).
We further define the shock microphysical parameters B and e as the fraction of the energy density 2 1+i ρ CSM v 2 that goes into the post-shock magnetic field energy density E B and the relativistic electron energy density E e . We introduce the numerical factor i to parametrize the different definitions of the reference energy density that appear in the literature (see Appendix C for details). For our modeling, we keep the dependence on i explicit. For reference, the traditional parametrization of C98 corresponds to i = 1. We note that in the SN literature there are at least three other definitions of the shock microphysical parameters that we discuss in detail in Appendix C. It is important to keep in mind that these different definitions lead to different values of inferred parameters (most notably the mass loss,Ṁ ) even when the same shock microphysical parameters values are assumed ( §3.2), as we also show in Appendix C. The energy densities are defined as follows: and where ρ CSM is the unshocked CSM density and v ≡ dR/dt is the FS velocity.
The frequency above which electrons cool efficiently through synchrotron radiation over a dynamical timescale t dyn is ν c , defined as (Rybicki & Lightman 1986): where σ T is the Thomson cross section, m e is the mass of the electron, e the charge of the electron, and c is the speed of light, all in c.g.s. units. For our purposes t dyn can be identified as the time since the SN explosion. In radio SNe, B ≡ B(t), which adds another dependency of ν c on time.
The single-electron synchrotron spectrum F ν peaks at ν = γ 2 e ( eB 2πmec ) (e.g., Rybicki & Lightman 1986), which allows us to define the characteristic synchrotron frequency: Following , under the assumptions that particle kinetic energy is conserved across the shock discontinuity and a pure H composition: whereγ is the average electron Lorentz factor and m p is the proton mass. For a dN/dγ e power-law distribution of electronsγ ≈ γ m , which, substituting into Equation 4, leads to: For typical radio SN parameters (e.g., from Soderberg et al. 2006a), B ∼ 1 G, v ∼ 0.1 c, it follows from Equations 3 and 6 that ν c 100 GHz and ν m 1 GHz. Taking SN 2003bg as an example, at t * ≈ 35 d since explosion Soderberg et al. (2006a) infer v(t * ) ≈ 0.13 c and B(t * ) ≈ 0.9 G, which implies ν m (t * ) ≈ 2.5( e /0.1) 2 × 10 −2 GHz and ν c (t * ) ≈ 2.6 × 10 2 GHz. It follows that for young radio SNe ν m < ν c and the only relevant spectral break frequency that typically enters the 0.1 − 100 GHz radio spectrum is ν sa .
In order to derive physical quantities from a single SED, we follow Rybicki & Lightman (1986) and C98. The specific intensity I ν can be expressed as follows: where S ν (ν sa ) = c 5 c 6 (B sin θ) −1/2 ν sa 2c 1 and y ≡ ν νsa . The constant c 1 and the parameters c 5 (p), and c 6 (p) are defined in Pacholczyk (1970) and provided in full in Appendix A. As these parameters are functions of p, it is typical in the radio SN literature to approximate c 5 ∼ c 5 (p = 3) and c 6 ∼ c 6 (p = 3), irrespective of what the best-fitting p is for the SED. In this work we self-consistently account for the dependency of c 5 and c 6 on the electron power-law index p. Additionally, we leave B sin θ explicit in our derivation (as opposed to approximating B sin θ ≈ B), as further physical quantities will not maintain the same exponential relationships between B and sin θ. 2 In Equations 7, 8, and 10 above, the radio emitting region is assumed to be a disk in the sky of radius R, thickness s, and volume V = πR 2 s.
The filling factor f is the fractional volume of a sphere of radius R such that the emitting volume is Canonically, a value f = 0.5 is adopted, but we leave f explicit in this work. The number density distribution of relativistic electrons per unit volume in energy space can be expressed as dN/dE = N 0 E −p , with N 0 as the normalization. This relates to e dγ e and E l is the lowest energy of the electrons accelerated to a power law, E l = γ m m e c 2 . For p > 2, the energy density in relativistic electrons of Equation 2 can be also expressed as We may now express the way Equation 12 relates to B: The flux density F ν of a source of radial extent R and constant I ν over the subtended solid angle, Ω = Considering the asymptotic limits ν ν sa (optically thick) and ν ν sa (optically thin), and using the equations above, C98 derives: 2 See Yang & Zhang (2018) for a recent discussion regarding the flux density dependence on the pitch angle θ distribution of the radiating electrons.
which scales as F ν,thick ∝ ν 5/2 , and the optically thin flux density F ν,thin ∝ ν −(p−1)/2 is: (15) We define a break frequency parameter ν brk as the spectral frequency at which the asymptotically thick and thin flux densities (Equations 14 and 15) meet: F ν,thin (ν brk ) = F ν,thick (ν brk ) ≡ F brk . After combining Equations 13, 14, and 15 and solving for B and R, we find the results below. We provide coefficients containing several digits to avoid rounding errors if a reader were to use these equations as provided.
For direct comparison to C98, we evaluate B and R at typical values of the physical parameters for radio SNe, i.e., p = 3, E l = m e c 2 , and sin θ ≈ 1 (e.g., Weiler et al. 2002;Berger et al. 2002;Soderberg et al. 2006a;Chevalier & Fransson 2006;Corsi et al. 2014 We note that for the case of a monatomic gas and a strong shock, the post-shock thermal energy is U th = 9 8 V ρ CSM v 2 (e.g., Petropoulou et al. 2016, see Appendix C for a comparison). Substituting p = 3, E l = m e c 2 , and sin θ ≈ 1 as above, our Equation 21  3.2. Radio Synchrotron Emission in SNe: environment density ρ CSM and pre-explosion mass-loss rateṀ A physical quantity of interest is the density of matter in the pre-shocked circumstellar medium, which directly traces the mass-loss history of the progenitor system before the star's death. We derive the density of the unshocked CSM, ρ CSM , by rearranging Equation 1 and substituting for B and R (Equations 16 and 17). Additionally, we follow C98 in assuming the evolution of the shock radius with time follows a power law R ∝ t q (either globally or locally), such that v ≡ dR/dt = qR/t (see Appendix C, and R(t) in Figure 3 for an example). Thus, For a shock that originates from the interaction between the outer layers of dense SN ejecta (which has a power-law distribution of density in radius ρ ej ∝ R −n ) and a power-law density distribution of CSM (ρ CSM ∝ R −s ), the shock radius evolution is described by the self-similar solutions in Chevalier (1982). 3 R(t) ∝ t n−3 n−s , which implies q = n−3 n−s . We note this relationship only holds in the regime of the self-similar solution, i.e., so long as n > 5 and s < 3 (Chevalier 1982). For a core-collapse SN originating from a compact, massive, stripped progenitor such as a Wolf-Rayet (WR) star n ≈ 10, and the CSM created from massive stellar winds scales with s = 2, which leads to a value of q = 7/8 that is typically applied in the literature (see e.g., Chevalier & Fransson 2006). From Equation 22, and a radio spectrum, the inferred pre-shock CSM density scales as ρ CSM ∝ 1+i q 2 . Introducing the mean molecular weight per electron parameter µ e and proton mass m p for the CSM, the electron number density in the unshocked CSM is n e = ρCSM mpµe , or: 3 The power-law density of the unshocked medium, s, is not to be confused with the model smoothing parameter s in Section 3.3. We use this notation for ease of comparison to Chevalier (1982).
During its lifetime, a star loses material to its surroundings and sits in a polluted environment of its own debris. For stars that lose mass through a constant mass-loss rateṀ and wind velocity v w , the CSM density is related toṀ as ρ CSM =Ṁ 4πR 2 vw . Combining Equations 17 and 22, the mass-loss rate at a radius probed by the forward shock thus satisfies the relation: Similarly to ρ CSM , both n e andṀ scale as ∝ i+1 q 2 . We emphasize that radio observations constrain the ra-tioṀ /v w by providing ρ CSM and that an independent estimate of v w , from optical spectroscopy for example, is needed to resolve for the model degeneracy. For CSM environments that deviate from a wind profile, Equation 24 provides an estimate of the effective mass-loss rate in units of v w .
When we evaluate our Equation 23 instead at q = 1 and i = 5/3 and assume the material is fully ionized hydrogen (such that µ e = 1), we recover Equation 16 Equation 24 reduces to: , which is Equation 23 in Ho et al. (2019) (in terms of flux and distance, rather than luminosity) once one accounts for the facts that: (i) Equation 23 in Ho et al. (2019) reports the wrong scaling for the filling factor f ; (ii) the reported normalization is also not correct ( 10 times larger than the correct value). See also footnote 4 of Yao et al. (2021). Our Equation 27 above gives the correct parameter scalings and normalization. We conclude with considerations regarding the impact different values of i and q will have on pre-shock CSM density and mass-loss rates. In the radio SN literature, a diversity of parameterization choices map to values of i and q that we explore in Appendix C. It follows that inferred ρ CSM values are not directly comparable to one another, even when the same microphysical parameters e and B are chosen. For example, using i = 5/3 and q = 1 (as in Ho et al. 2019) leads to a factor of  (2020) is equivalent to i = −1/2 and q = 1; this definition leads to 1+i q 2 = 1/2. Amongst a range of assumed i and q values in the literature (Appendix C), the ratio 1+i q 2 spans values 1/2 to 8/3, leading to reported values of ρ CSM (as well as n e andṀ ) that are systematically different up to a factor ≈ 5. We thus caution the reader regarding direct comparisons ofṀ (or n e or ρ CSM ) values from independent works, even when the same shock microphysical parameters are used.

Radio SED fitting
In this section we describe the fitting procedure used on the radio SEDs to extract the physical parameters of the emitting source based on the physics described in the previous sections. We model the radio SED produced by synchrotron emission described in §3.1 with a smoothed broken power law (BPL) of the form: This functional form is widely used in the radio SN and non-jetted TDE literature (e.g., Soderberg et al. 2006a, Soderberg et al. 2006cAlexander et al. 2020 and references therein). The smoothing parameter s represents how abruptly the SED turns over from the "optically thick" slope (α 2 ) to the "optically thin" slope (α 1 ). For an SSA spectrum with ν m < ν sa α 1 = −(p − 1)/2 and we expect α 2 = 5/2 (Equations 14 and 15). Sharper transitions between slopes are associated with smaller s values. For the spectral regime of interest, ν brk can be identified as ν sa . F brk represents the peak flux density where the asymptotic power-laws describing the optically thin and the optically thick spectrum meet (i.e., the F brk parameter in Equation 28 is the break flux density that appears in the equations of the previous sections).
In general, ν brk ≡ ν brk (t) and F brk ≡ F brk (t). For a radio SN, the typical behavior is to observe a migration of ν brk to lower frequencies as time progresses and the emitting region expands over a larger volume, leaving a broader range of higher frequencies optically thin to synchrotron radiation. We fit each SED without imposing a pre-defined temporal evolution of the parameters (in the SN literature a power-law evolution of ν brk (t) and F brk (t) is typically assumed, (e.g., Soderberg et al. 2006a)) allowing us to reconstruct the actual density profile of the CSM created by the star before explosion, and hence the true mass-loss history of the progenitor. An example of this procedure on our data is represented in Figures 1 and 2, and the full SED evolution is Figure  10 in the Appendix.
The SED power-law slopes α 1 and α 2 are free parameters that can, in principle, take different values in each of the SEDs. Our fitting procedure provides the flexibility to jointly fit all the available SEDs with common values of α 1 and α 2 , along with the smoothing parameter s. At late times when the peak frequency has migrated below our observable frequency range, we fit the SED to a single-sided power law and include this information in a joint-fit of the slopes. The information from singlesided SEDs provide lower and upper limits on physical parameters. Our fitting procedure thus does not assume any value for p; instead p can be derived from the bestfitting optically thin slope α 1 . We use the BPL function of Equation 28 with lmfit and perform a least squares minimization.
We point out that an important source of confusion in the radio SN literature is that the F brk of Equations 16−23 is often incorrectly identified as the peak of the smoothed broken power-law flux density, as opposed to the intersection of the two asymptotic power-laws. This is important because this implies that the true F brk is underestimated by a factor of 2 s that propagates into all estimates of the physical parameters (see Equations 16−23).

Radio SED fitting of SN 2004C
The radio data of SN 2004C are reported in Table 3. We consider data acquired within date ranges δt/t < 0.025 as part of the same SED. This procedure identifies 28 individual radio SEDs spanning the range δt = 41.41 − 1938.22 days since explosion. We modeled each SED with the BPL function of Equation 28, and we jointly fit all the SEDs with the same α 1 , α 2 and s value, while allowing F brk and ν brk to evolve from one SED to the other. We explored with our software four cases: (i) fixing both the thick slope α 2 = 5/2 and s = −1 ("BPL1"), (ii) fixing only s = −1 ("BPL2"), (iii) fixing only the thick slope α 2 = 5/2 ("BPL3"), and (iv) allowing both slopes and s vary ("BPL4").
We compare these models using an F-Test and report our results in Appendix B, Table 2. We adopt the standard p-value threshold of 5% to reject a given model, and χ 2 when the degrees of freedom are equal. We identify BPL1 (s = −1, α 2 = 5/2) as the favored model by these statistics. The best-fitting model parameters for BPL1 are reported in Table 1. Figures 1, 2, and 10 show the best fitting model and the data.

Inferred physical parameters of SN 2004C: p, R(t), B(R)
We use the best-fitting parameters of Table 1, their uncertainties, and covariance information to estimate the physical parameters of the shock launched by SN 2004C and the circumstellar environment with which it interacts. First, the optically thin slope α 1 = −(p − 1)/2 implies a power-law index of the electron distribution p =2.92 ± 0.06, in-line with typical values inferred from SN from the radio modeling of SNe p ≈ 3 (e.g., Chevalier & Fransson 2006;Soderberg et al. 2006a) and recently interpreted by Diesing & Caprioli (2021). Values of p show consistency amongst all our fit models, falling within ∼ 0.1 of one another. Next, we estimate R(t) and B(t) using Equations 17 and 16, respectively. In our analysis we adopt i = 1 and we infer q = 0.77 ± 0.04 (Figure 3). We present our results for the case of equipartition ( e = 1/3, B = 1/3) and for our fiducial case e = 0.1 and B = 0.01. We note that the radius and magnetic field scale as R ∝ , ρ CSM (r) ∝ Table 1. SED best-fitting parameters for our baseline model (BPL1), which has an assumed smoothing parameter and optically thick slope s = −1, α2 = 5/2 and an optically thin slope α1 = −0.96 ± 0.03. "N/A" entries show SEDs containing only one point and could not be fit to a one-sided power law.   13+2p) . Changes to the distance by a factor N moderately affect our results; for p ≈ 3, R will change by a factor of N 0.9 anḋ M by a factor of N −0.8 . Figure 3 shows the evolution of the forward shock radius and shock velocity with time, implying mild deceleration of the shock as it interacts with material in the ambient medium. We find R(t) ∝ t 0.77 ± 0.04 over the distance range (8.1 ± 0.5) × 10 15 − (3.4 ± 0.4) × 10 16 cm using observations acquired 41 − 1938 days since explosion. With the VLBA we inferred R = 3.8 +0.6 −1.0 × 10 17 cm at 1931 days since explosion. The extrapolation of our best-fitting FS radius derived from VLA observations to the time of the VLBA observations is R = (1.5 ± 0.5) × 10 17 cm, which is consistent, but slightly lower, than the interferometry measurement at the 3 σ level (Fig. 3). We remark that this is an extrapolation over almost one order of magnitude in time. Including the VLBA measurement in the R(t) fit we find q = 0.81±0.03, which is consistent within uncertainties with the radius expansion rate that we infer at < 242 days since explosion. The forward shock has an average shock velocity of ≈ 0.06c throughout the duration of our observations. For a power-law evolution of the shock radius with time, R(t) ∝ t q the velocity is v = q R t . Instead of assuming a specific progenitor type and CSM density profile to compute q using the self-similar solutions by Chevalier (1982), we allow our code to be as agnostic and self-consistent as possible: we derive q by taking the slope of R(t) (solid line in Figure 3). Other fitting methods that hold the SED slope constant for synchrotron emission or allow all BPL fit parameters to vary lead to reasonably similar values. Appendix C demonstrates the effect this has on derived physical quantities. Figure 3 presents the inferred evolution of the magnetic field with shock radius B(R): we find evidence for a clear change in the B evolution at R br ≈ (1.96±0.10)× 10 16 cm, with a flat profile B ∝ R −0.13 ± 0.08 at R < R br evolving to B ∝ R −1.10 ± 0.07 at R > R br . As a reference, SN shocks expanding in wind-like media where e and B are not time-varying quantities are expected to show B ∝ R −1 (Chevalier 1998). The location of this "break radius" is largely independent from the choice of microphysical parameters (Figure 3), due to the very mild dependence of R on e and B .
We end this section by noting that in the context of SSA radio spectra for which the data provide constraints on only one break frequency ν sa , the parameters p, R and B are the three physical properties that can be independently constrained from the data. Any other parameter, such as those discussed in the next section, are derived from combinations of p, R and B. In our calculations of the physical parameters that follow, we utilize the value of q associated with the powerlaw slope of R vs. t. We use Equations 20, 22, and 24 to calculate U (t), ρ CSM (R) andṀ (R). We find that in equipartition the blastwave energy U (t) increases with time from (1.60 ± 0.25) × 10 47 erg at 41 days since explosion to (1.9 ± 0.4)10 48 erg at 242 days since explosion (Figure 4, green points). This behavior is consistent with the progressive conversion of ejecta kinetic energy into shock internal energy, as the stratified ejecta decelerates into the environment. These values are lower limits on the true shock internal energy because of the likely deviation from equipartition. However, as long as the shock microphysical parameters are constant with time, the ratio between the final and initial U (t) is preserved, as is the temporal behavior of U (t) or ρ CSM (R).
With reference to Figure 5, our analysis highlights the presence of an internal region with roughly constant density extending to R br ≈(1.96 ± 0.10) × 10 16 cm, followed by a rapidly declining density profile with ρ CSM ∝ R −2.3 ± 0.5 at R > R br . For a nominal ejection velocity of the CSM material by the progenitor of v w = 1000 km s −1 , R br corresponds to a look-back time of ≈ 10 3.36 days (6.2 yr) prior to stellar collapse. Importantly, the ρ CSM (R) profile that we derived from the radio observations violates the expectation from a single wind density profile, which predicts ρ CSM ∝ R −2 . Our results point instead at the presence of a shell-like density structure ∼ 10 16 cm from the explosion site of SN 2004C not unlike other stripped-envelope SNe such as SN 2004dk, ASASSN15no, SN 2017dio, and SN 2014C (Margutti et al. 2017;Benetti et al. 2018;Kuncarayakti et al. 2018;Stroh et al. 2021;Balasubramanian et al. 2021, Brethauer, in prep.). Interestingly, the phenomenology of mass lost to the environment in the last moments of evolution extends all the way to SLSNe (i.e., iPTF15esb, iPTF16bad, SN 2017ens), which presumably have different progenitors (Yan et al. 2017;Chen et al. 2018).
Radio observations constrain the ratio ρ CSM ∝Ṁ /v w but do not independently constrainṀ and/or v w . It follows that the shell-like density profile around SN 2004C can be the result of a time-varying mass-loss rateṀ (t) or a time-varying progenitor wind v w (t), or both.   value, these results imply that in this case the very last phases of evolution of the progenitor at t <10 3.36 days (6.2 yr) have been characterized by a rapidly decreasinġ M with time, withṀ reaching a peak value ofṀ ≈ (460 ± 80) × 10 −5 M yr −1 (for e = 0.1 and B = 0.01).
We note thatṀ vw ∝ e For the complementary case of a constantṀ ≈ 10 −5 M yr −1 (Figure 4), the observed density profile would instead imply a significant change in the wind velocity in the last thousands of years preceding core-collapse, with a rapid increase of one of magnitude in the final ∼ 3000 years (from 2 km s −1 to 18 km s −1 for our assumeḋ M unit value). These winds' velocities are clearly not consistent with a compact, single WR star.
In reality, neither of these two limiting cases might apply. The shell-like structure might be the result of a "wind-wind interaction", where faster, lighter winds interact with mass lost to the environment through slower winds of the preceding phase of stellar evolution, for example during the Red Supergiant (RSG) → WR transition, (Heger et al. 2003;Dwarkadas et al. 2010;van Loon et al. 2005;Mauron & Josselin 2011). We explore the "wind-wind interaction" scenario and additional explanations in the following section.

MASS LOSS DISCUSSION
The major result of modeling our radio observations is a radial profile in ρ CSM (R) that is robust against the choice of SED fit model (Figures 5 and 9). Most importantly, this profile reveals a structure in the CSM that is flat (ρ ∝ R −0.03 ± 0.22 ) out to (1.96±0.10)×10 16 cm, followed by an immediate, steep decline of ρ ∝ R −2.3 ± 0.5 .
The maximum mass-loss rate that can be supported by line-driven winds depend on metallicity Z and the luminosity L of the progenitor star as follows (e.g., Smith where L 6 ≡ 10 6 L . For solar composition, the hydrogen, helium, and metal fractions by mass are X = 0.7381, Y = 0.2485, Z = 0.0134, respectively (Asplund et al. 2009 impliesṀ max,lines ≈ 7 × 10 −5 M yr −1 , but SN 2004C shows a maximum mass-loss rate of (460 ± 80) × 10 −5 M yr −1 (assuming v w = 1000 km s −1 , B = 0.01, and e = 0.1), two orders of magnitude greater than the limit of single line-driven winds. We conclude that a single fast wind (v w = 1000 km s −1 ) from a compact progenitor cannot sustain the large mass-loss rates inferred for SN 2004C. We now compare SN 2004C to mass-loss rates observed in massive stars. The densities measured for SN 2004C overlap in theṀ vs. v w phase space with RSG and Luminous Blue Variable (LBV) progenitors, because RSGs have significantly slower winds and LBVs are not powered by line driven winds ( Figure 6). How- Figure 5. CSM density probed by the forward shock (upstream) as revealed by synchrotron emission at radio frequencies using our re-derived Equations 17 and 22. Purple points refer to values calculated using e = 0.1 and B = 0.01; and green points refer to values calculated in equipartition, e = B = 1/3; they will have the same uncertainties as the purple points. Diagonal lines: comparison to constant mass loss assuming vw = 1000 km s −1 . The mass loss-rate of SN 2004C is much greater than can be explained by line-driven winds (red dashed line,Ṁ = 10 −4 M yr −1 ) even at extreme luminosities L = 10 6 L . The yellow dashed line denotingṀ = 10 −2.5 M yr −1 is to assist the reader's estimate of the slope. Top abscissa: time before explosion material was deposited in the CSM, assuming vw = 1000 km s −1 and free expansion. ever, both cases expect copious H in the ejecta, which was not observed for SN 2004C (Shivvers et al. 2017). Assuming a WR wind speed of 1000 km s −1 the inner density profile of SN 2004C corresponds to a mass loss rate of (54 ± 8) × 10 −5 M yr −1 at radius R =(8.1 ± 0.5) × 10 15 cm (with B = 0.01 and e = 0.1; Equation 24, Figure 4). WRs typically span a range of mass-loss ratesṀ ≈ 10 −5.6 − 10 −4.4 M yr −1 , which is significantly lower than observed in SN 2004C for the same wind velocity (Crowther 2007 Figure 6). Similarly, the region of SN 2004C that overlaps with the parameter space partitioned to RSGs can be ruled out, as RSGs are H-rich and we did not see H in the spectra of SN 2004C (Shivvers et al. 2017). This comparison to known stellar mass-loss rates strengthens our previous conclusion that a single line-driven wind cannot produce the CSM material surrounding SN 2004C; and instead we must invoke slower winds or multiple line-driven winds ( Figure 6).
As an alternative, continuum-driven winds can power larger mass-loss rates than line-driven winds due to their dependence on the electron scattering cross-section rather than the cross-sections of metals in the stellar atmosphere. In principle, this implies that continuumdriven winds are not as sensitive to metallicity. However, even if continuum-driven winds could manage to supply a constant mass-loss rate at the appropriate magnitude, we find the single-wind mechanism as a whole is in tension with the observed CSM density profile surrounding SN 2004C. A single wind (line-or continuumdriven) with a constant mass-loss rate would create an environment density that scales with ρ ∝ R −2 , rather than the measured ρ ∝ R −0.03 ± 0.22 that transitions to ρ ∝ R −2.3 ± 0.5 ( Figure 5).
Instead, the density profile that we infer from radio observations requires us to consider more exotic forms of mass loss or wind-wind interaction scenarios. We investigate the possibility the CSM could have been sculpted by the interaction of multiple winds or time-varying winds in §5.1. We explore the possibility that this shell of material was ejected by a single event (i.e., a shell ejection scenario) in §5.2.

Wind-Wind Interaction and Time-Varying Continuum Winds
Wind-wind interaction can occur in the circumstellar medium where a faster-moving wind encounters a slower-moving wind launched by the star at an earlier date. For reference, free expansion of a wind would produce a density profile ρ ∝ R −2 , which is not too dissimilar from the outer shell density profile we observe in SN 2004C, where ρ ∝ R −2.3 ± 0.5 ( Figure 5). Figure 6. Wind velocity vs. mass-loss rate phase space (see §5 for full discussion). Blue Shaded Region: Mass loss rates of main sequence stars and subgiants. Rates are taken from Seaquist et al. (1993). Light Green Shaded Region: Mass loss rates of Wolf-Rayet stars from Crowther (2007). Brown Shaded Region: Mass loss rates from LBV continuum winds (Smith 2014;Terreran et al. 2019, and references therein) Teal Shaded Region: Wind speeds and mass loss rates of red supergiant stars (Seaquist et al. 1993). Gold Shaded Region: Extreme winds possible in red supergiants (de Jager et al. 1988;Marshall et al. 2004;van Loon et al. 2005). Purple Shaded Region: Observed mass losses of LBV eruptions (Smith & Owocki 2006). Orange Diagonal-Dashed Lines demarcate the parameter space constrained by our observations ofṀ /vw in SN 2004C. Note that although RSG winds fall within the realm of possibility, they have not been identified as SN Type IIb progenitors.
In the SN literature, wind-wind interaction models have been invoked to explain the "anomalous" density profiles found around SN 1996cr (Dwarkadas et al. 2010) and SN 2001em (Chugai & Chevalier 2006). Chugai & Chevalier (2006) speculate that the shell of material around SN 2001em resulted from a combination of latestage mass transfer in a binary system followed by a RSG → WR transition. In this instance, both stages could could have polluted the environment and subsequently swept up the CSM as winds accelerated, then pushed the shell out to a radius of ∼ 10 16 cm. In the case of a single star, the transition to a compact WR phase (which is associated with a faster wind) is expected to take place ≈ 0.5 − 1 Myr before explosion, much earlier than the approximate shell age of 1000 years in SN 2004C (Dwarkadas et al. 2010;Margutti et al. 2017). However, for some binaries, the envelope removal of the primary star by binary interaction can happen much closer to the epoch of core collapse, as it was proposed for SN 2014C (Margutti et al. 2017;Milisavljevic et al. 2015). A similar connection between progenitor envelope and binary-induced mass loss was also proposed for Type IIb SNe by Maeda et al. 2015. It is possible that SN 2004C resembles this scenario, requiring a binary system and a shorter WR phase than what single-star evolution allows. In this case, faster winds sculpt slower-moving winds in the CSM and the star consequently explodes as a stripped-envelope SN.
Irrespective of the physical origin of the CSM shell, sustained mass loss in the final decades before explosion can significantly modify the stellar progenitor structure at the time of explosion. It is thus instructive to review our current knowledge of stellar progenitors of Type IIb SNe. Based on pre-explosion imaging with HST (Smartt 2015), the prominence of the cooling-envelope phase in the optical light-curve (e.g., Arcavi et al. 2011), and the dynamic of the fastest outflows launched by the explosions as constrained by radio observations (Soderberg et al. 2012) 5 , the stellar progenitors of Type IIb SNe show diversity in their structure, from extended yellow hypergiants like the recent SN 2016gkg (Kilpatrick et al. 2017;Tartaglia et al. 2017;Piro et al. 2017;Orellana et al. 2018;Sravan et al. 2018;Kilpatrick et al. 2021) to more compact progenitors with R star < 50 R , as in the case of Type IIb SN 2008ax (e.g., Folatelli et al. 2015;Yoon et al. 2017;Groh et al. 2013;Roming et al. 2009;Crockett et al. 2008). Based on the current constraints on masses and radii of Type IIb progenitors (compiled from Sravan et al. 2020;Gilkis & Arcavi 2022 and references in both), the inferred escape velocities lie in the range 45 − 300 km s −1 , which maps into typical stellar winds velocities of 135 − 900 km s −1 (order of magnitude estimate). We note that the current sample of Type IIb SNe with detected progenitors is biased toward the more extended stars, and that compact stars that are more difficult to detect are associated with faster winds. In this sense, the estimates above should be viewed as a likely lower limits on the true wind velocity of SN 2004C just before explosion. We thus conclude that for SN 2004C, v w 100 km s −1 at the time of the explosion, which implies large mass-loss ratesṀ 10 −4 M yr −1 .

Shell Ejection Scenarios
Unlike steady winds, shell ejection scenarios mark an abrupt, sudden, and violent expulsion of mass from the progenitor, often triggered by an internal deposition of energy (Quataert & Shiode 2012;Woosley & Heger 2015;Woosley 2017;Leung et al. 2019;Leung & Fuller 2020;Wu & Fuller 2021). In SN 2004C, we observe a structure in the CSM between (8.1 ± 0.5) × 10 15 cm and (3.4 ± 0.4) × 10 16 cm. Assuming a constant wind speed of 1000 km s −1 , our model suggests this structure was produced over the course of 10 2.97 days (2.6 yr) -10 3.59 days (10.8 yr) before core collapse ( Figure 5). Therefore, we find evidence pointing to a shift in the behavior of the progenitor star in the final 10 −5 − 0.01% of its lifetime.
Traces of violent, substantial outbursts depositing material ∼ 10 16 cm from the star at the time of explosion have been reported across multiple SN types. Hydrogenpoor SNe Types Ib/c, IIb (such as SN 2004C), and even SLSNe-I show evidence of CSM shells at this distance as we continue to monitor these events at late-time with radio observations (Yan et al. 2017;Lunnan et al. 2018;Mauerhan et al. 2018;Benetti et al. 2018;Kuncarayakti et al. 2018;Kundu et al. 2019;Gomez et al. 2019;Prentice et al. 2020;Jin et al. 2021;Gutiérrez et al. 2021;Maeda et al. 2021;Jacobson-Galán et al. 2021b;Stroh et al. 2021;Dong et al. 2021 Auchettl in prep., and the sampled population in Hosseinzadeh et al. 2021).
An extreme case of continuum-driven winds is an LBV outburst (Humphreys et al. 2017(Humphreys et al. , 2019. These hot, blue stars quickly ionize stellar material, creating a large population of free electrons. In turn, the electrons serve as cross sections for radiative absorption and form metallicity-independent continuum-driven winds (Smith 2014). Outbursts become progressively more extreme the closer a star gets to its Eddington limit, consequently providing a spectrum of mass-loss events and outbursts. Such outbursts have been posited to source quasi-periodic mass loss observed in SNe at radio wavelengths (e.g. Kotak & Vink 2006). For example, η Car's Great Eruption in the nineteenth century involved a loss of about 12 − 20 M of stellar material. The majority of this material was lost over five years, but a lower limit of aboutṀ = 0.5 yr −1 is estimated for the average over the ≈20 year duration of the eruption (Smith & Owocki 2006). Another example is the 1600 eruption of P-Cygni, wherein 0.1 M was lost, implying a rate ofṀ = 10 −2 M yr −1 due to continuum winds (Owocki et al. 2004;Smith & Owocki 2006). These specific ejected mass and mass-loss rates are significantly larger than the M ≈ 0.021 M CSM shell observed in SN 2004C (Figures 5 and 6). While an LBV progenitor was suggested by quantitative modeling of the very early spectra of the Type IIb SN 2013cu (Groh 2014) and the transitional Type Ibn/IIb SN 2018gjx (Prentice et al. 2020), SN 2004C, shows a clear difference from the historical LBV eruptions because its progenitor underwent core-collapse shortly after shell ejection. Furthermore, mass-loss rates of > 10 −3 M are rare and have not been observed in living LBVs within the last few decades (Davidson 2020). In the following we explore shell-ejection mechanisms that are tied to the requirement of imminent core collapse.
In the case of "Pulsational Pair-Instability" SNe (PPISNe), stars with He cores M > 30 M can undergo electron-positron pair production, removing pressure in their cores that would otherwise support the star against the crush of gravity. Consequential pulses rippling through the star can launch enormous amounts of mass at the stellar surface, or even the H envelope itself (0.3 − several M , Woosley 2017). Because PPISNe exhibit multiple eruptions of unstable mass loss that mimic SNe, there is no way to know for certain if an explosion is a PPISN without observing a second eruption (Woosley 2017). Additionally, pre-SN core masses large enough for PPISNe require large Zero-Age Main Sequence (ZAMS) masses (e.g., 70-260 M ), greatly favoring low-metallicity (sub-solar) environments (Woosley 2017). Because the host galaxy of SN 2004C is about solar-metallicity (with possibly a higher metallicity concentration at the site of the explosion (Dittman, private communication) we rule out the PPISN mechanism as an underlying explanation to our observations. Another mechanism we consider is gravity waves, or wave heating. Proposed by Quataert & Shiode 2012, hydrodynamic waves during late-stage nuclear burning can excite violent convection. At the onset of core carbon burning, neutrino cooling is unable to fully dissipate the generated energy and as a result, the extreme convection launches internal gravity waves with energies ∼ 10 7 L (Fuller & Ro 2018). Barring large amounts of nonlinear wave breaking and damping throughout the body of the star (explored thoroughly in Wu &Fuller 2021 andRo 2018), the wave energy that survives to the surface can unbind mass, accelerate it above the escape velocity, and result in observable mass loss.
Fuller & Ro 2018 explore the wave heating mechanism in stripped and partially stripped progenitors at solar metallicity with 1D MESA models. Their models result in Type Ib/c and IIb SNe and place CSM material too nearby (within 6 × 10 14 cm) to apply to SN 2004C. Similarly, Wu & Fuller (2021) conclude that stars of M < 30 M could exhibit gravity-wave driven mass-loss outbursts in the years to decades before explosion, which is too late in the life of a star to be relevant for SN 2004C. Specifically, their H-poor models with M ZAMS = 36 − 40 M end their lives as M ∼ 15 M WR stars. In this case, the energy generated by gravity waves was sufficient to eject 10 −2 M , but the ejected mass only reached r ∼ 10 14 cm before core collapse, significantly smaller than the inner radius of (8.1 ± 0.5) × 10 15 we observe in SN 2004C. Their model with the most distant CSM material (r ∼ 10 15 cm) was an 11 M yellow supergiant (YSG) progenitor, which lost 1 M over the 10 years prior to core collapse. We conclude that if related to wave-driven instabilities, the phenomenology that we observe in SN 2004C requires this mechanism to be active earlier than nuclear burning stages.
As a result of pre-explosion imaging coupled with modeled interactions, binary star systems are one of the favored progenitor systems of SNe IIb (Sravan et al. 2019;Yoon et al. 2017;Ryder et al. 2018;Fox et al. 2014;Zapartas et al. 2017). The amount of stripping the primary star undergoes from the secondary companion strongly depends on the initial masses of the system, the initial separation of the two bodies, and their metallicities (Smith & Arnett 2014;Yoon et al. 2017). It is believed that about 2/3 of massive stars are members of binary systems that interact before stellar death, and both early and late Case B mass transfer (when the donor star is evolving into an RSG) can induce primary stars to become blue hypergiants, yellow hypergiants, and red supergiants -all progenitors of SNe Type IIb (Smith & Arnett 2014;Yoon et al. 2017;Ryder et al. 2018;Fox et al. 2014). We find that we cannot rule binary stripping and its possible effects on mass-loss winds throughout the lifetime of SN 2004C, because binary interaction is capable of ejecting the entire hydrogen envelope of a star over a variety of timescales (Yoon et al. 2017;Ryder et al. 2018;Fox et al. 2014).

CONLCUSIONS
In this paper, we began with the synchrotron derivation of the radio emission from a SN shock by Chevalier (1998) and present here the expressions for the postshock magnetic field B, forward shock wave radius R, shock internal energy U , circumstellar density ρ CSM , and associated mass-loss rateṀ , generalized for any value of the power-law index of the electron density distribution, p. As part of our generalized formalism, we reviewed the different definitions of the shock microphysical parameters e and B that are used in the radio SN literature and we show in detail how hidden assumptions can lead to systematically different inferences of the shock and environment parameters, even when the same microphysical parameters values are adopted. These discrepancies are summarized and quantified in Appendix C. We applied our generalized formalism of Equations 16 -24 to a collection of 7.7 years of multi-frequency radio data of SN 2004C and we showed how this formalism leads to inferences on the forward shock dynamics that are in agreement with our direct VLBI measurements. Our modeling revealed the presence of a dense shell of CSM material containing a mass of 0.021 M at a distance of about 10 16 cm from the explosion. Our major results can be summarized as follows: • The shock microphysical parameters e and B have been defined in different ways in the radio SN literature. In this work we quantitatively show the impact of the most common definitions of these parameters on the inferred shock and environment parameters; and we present a parametrization that will enable a direct comparison of the inferences from different works in the literature. We emphasize that in the specific case of the mass-loss rates, different, typically hidden, definitions lead toṀ values that can differ up to a factor of ∼ 5 -even for the same choice of value for the shock microphysical parameters. These are compared in Appendix C.
• Our radio modeling reveals a density profile with a flat structure in the CSM (ρ ∝ R −0.03 ± 0.22 ) out to (1.96±0.10)×10 16 cm, and an abrupt transition to ρ ∝ R −2.3 ± 0.5 at larger radii. The most impor-tant result from the modeling of our radio observations is that a radial profile density ρ CSM (R) in SN 2004C that exists independently of the fitting model ( Figures 5 and 9).
• If the 0.021 M of CSM material around SN 2004C was created by a shell ejection event, the implied ejection epoch is 128 − 2.5 years prior to explosion for an ejection velocities between 20−1000 km s −1 .
• Of the mass-loss mechanisms we discuss in §5, we assert a single, line-driven wind with constant speed and mass-loss rate could not have sculpted the structure observed in the CSM of SN 2004C ( Figure 6). Instead, multiple line-driven winds, a single time-varying wind, or a time-varyingṀ could be responsible. We cannot rule out binary interaction as the underlying cause. While the extreme limits of gravity waves/wave heating and LBV mass-loss rates marginally overlap with the uncertainties of our calculated rates, current models predict material either too close to the star or released too slowly.
The connection between the shell characteristics, mass-loss mechanism, and timing of mass-loss events events highly suggests a fundamental aspect of stellar evolution in evolved massive stars at work. Other SNe across different explosion types exhibit similarly substantial overdensities in their environments at r ∼ 10 16 cm from their explosion site (Margutti et al. 2017;Hosseinzadeh et al. 2021;Benetti et al. 2018;Yan et al. 2017;Leung et al. 2021a;Jacobson-Galán et al. 2021b;Kuncarayakti et al. 2018;Chen et al. 2018;Jacobson-Galán et al. 2021a;Stroh et al. 2021), possibly suggesting a shared mechanism. By observing changes in radio emission from the SN shock wave to probe the CSM density, we gather observational evidence that a rich and complex mass-loss history exists in stars, implying episodic mass loss in their final moments while they are on the verge of collapse. Consequently, the inability to conclusively point to any particular mechanism as the underlying cause highlights the need to model the final moments of late-stage core burning in H-poor progenitors. ACKNOWLEDGMENTS L.D. heartily thanks Avni Coda for her contribution to all data visualization and gratefully acknowledges the help of tablesgenerator.com for importing .csv files as L A T E X tables with minimal headache. L.D. is grateful for the partial financial support of the IDEAS Fellowship, a research traineeship program funded by the National Science Foundation under grant DGE-1450006. The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. The Margutti team at UC Berkeley and Northwestern is supported in part by the National Science Foundation under Grant No. AST-1909796 and AST-1944985, Table 2 is a statistical comparison of the performance of different SED models. We first employ the F-Test, which begins by calculating the F-statistic, generally expressed as: .
Model 1 has a larger number of degrees of freedom (d.o.f.) and fewer free parameters (M params ). Our Null Hypothesis (H 0 ) is that Model 1 is the preferred fit, following the argument of Occam's Razor. We reject H 0 if the p-value is p < 5%.
Our radio data set of SN 2004C has a number of data points N data = 94. The model names are abbreviated with "BPL", referring to their broken power law shape (Equation 28) and they represent a joint fit of 28 radio SEDs which have 56 free parameters (i.e. ν pk , F pk , one for each SED) with an additional 1 to 3 shared parameters between them (slope α 1 , slope α 2 , and smoothing parameter s). A full discussion of fit parameters and methods is provided in §3.3 and 4.1. For the comparison of BPL2 to BPL3, the number of degrees of freedom is the same, and we use the χ 2 statistic. We identify BPL1 (s = −1, α 2 = 5/2) as the favored model by these statistics.
In the Table, Model 1 refers to BPL1, Model 2 refers to BPL2 etc. As an example of how to read the table, for BPL2 (which is Model 2), the column "F − Stat 1" reports the value of the statistics with respect to Model 1, while the column "p − value 1" reports the corresponding p-value. We also present Figures 7 and 8 as extensions of Figures 3 and 4, respectively; and we provide a version of Figure  5 in Figure 9. These figures contain all four models as well as the equipartition case of BPL1 for ease of comparison. We note that the choice in model can cause a spread amongst values of physical parameters, but the overall evolution of the inferred parameters remains the same. Additionally, we highlight that the largest variation is associated with the choice of micro-physical parameter values . Figure 7. The same as Figure 3, but with all models compared in Section 3.3 and Appendix B. Note that the normalization differs between models, but the overall evolution remains the same. Purple pentagons: Favored model where smoothing parameter is fixed at s = −1 and the slope α2 = 5/2 ("BPL1"). Yellow triangles: only smoothing parameter is fixed s = −1 ("BPL2"). Blue squares: only slope is fixed α2 = 5/2 ("BPL3"). Pink diamonds: All slopes and s are free to vary ("BPL 4"). Green points: BPL1 in equipartition ( B = e = 1/3). Top, left panel; blue line: Power-law fit R ∝ t q , where q = 0.77 ± 0.04 and appears in our generalized formulae as the variable q. Figure 8. The same as Figure 4, but with all models compared in Section 3.3 and Appendix B. Note that the normalization differs between models, but the overall shape of the evolution remains the same. Purple pentagons: Favored model where smoothing parameter is fixed at s = −1 and the slope α2 = 5/2 ("BPL1"). Yellow triangles: only smoothing parameter is fixed s = −1 ("BPL2"). Blue squares: only slope is fixed α2 = 5/2 ("BPL3"). Pink diamonds: All slopes and s are free to vary ("BPL 4"). Green points: BPL1 in equipartition ( B = e = 1/3). Figure 9. The same as Figure 5, but with all models compared in Section 3.3 and Appendix B. Note that the normalization differs between models, but the overall shape of the evolution remains the same. Purple pentagons: Favored model where smoothing parameter is fixed at s = −1 and the slope α2 = 5/2 ("BPL1"). Yellow triangles: only smoothing parameter is fixed s = −1 ("BPL2"). Blue squares: only slope is fixed α2 = 5/2 ("BPL3"). Pink diamonds: All slopes and s are free to vary ("BPL 4"). Green points: BPL1 in equipartition ( B = e = 1/3).

C. DEFINITIONS OF THE SHOCK MICROPHYSICAL PARAMETER B (AND E )
In the radio SN literature the parameter B (and e ) 6 has been defined in at least four different ways: where B is the magnetic field in the shocked CSM, v is the FS velocity, ρ CSM (ρ CSM,sh ) is the pre-shock (post-shock) CSM density, Γ is the adiabatic index, E th is the post-shock thermal energy density, and P 2 is the post-shock gas pressure. In other words, B is defined as a fraction of different energy densities that in this paper we parametrize as 2 1+i ρ CSM v 2 . The corresponding energy is U (Equation 20).
In the limit of a strong shock scenario the traditional Rankine-Hugoniot jump conditions apply: ρ CSM,sh = Γ+1 Γ−1 × ρ CSM . For an ideal monatomic gas, Γ = 5/3, which leads to ρ CSM,sh = 4×ρ CSM , E th = 9 8 ρ CSM v 2 , and E th = P2 Γ−1 = 3 2 P 2 . In this paper we use B 2 8π = 2 i+1 B ρ CSM v 2 as a generic definition, where the parameter i provides flexibility between the different definitions above. When one chooses definition C8, i in our equations becomes conveniently equivalent to the adiabatic index, Γ. Chevalier (1998) The obvious consequence of these different definitions of the microphysical parameters is that the inferred values of the physical parameters of the system (e.g., B, R, ρ CSM ) are not directly comparable among different works, even when the same B (and e ) values are assumed or derived. Derived values ofṀ can differ by as much as a factor of ∼ 5 in this way.
In addition, there are two approaches used in the literature for the FS velocity v ≡ dR/dt. The first approach utilizes an averaged FS shock velocity since explosion: A second approach assumes a (local or global) power-law evolution of the FS radius with time R(t) ∝ t q : At early times the typical non-relativistic SN FS is slowly decelerating and the assumption of a linear evolution of the shock radius with time implicit in Equation C9 is approximately correct. For SNe in the interaction phase expanding in a power-law density medium for which the self-similar solutions of Chevalier (1982) apply, Equation C10 provides the exact solution. Following Chevalier (1982), q = n−3 n−s , where ρ CSM ∝ r −s and the SN ejecta outer density profile is ρ ej ∝ r −n . For a wind density profile s = 2, and a compact massive star (i.e., n ≈ 10) such as typically assumed for the progenitors of H-stripped SNe, q = 7 8 . Definitions of Equations C5-C8 have been combined with either Equation C9 or C10 to solve for the parameter of interest ρ CSM (and, by extension, n e andṀ ). For example, Chevalier & Fransson (2006) adopt Equation C10 for q = 7 8 , while Ho et al. (2019) uses the definition of Equation C8 with the averaged FS velocity v avg of Equation C9. In this paper we assume a local power-law evolution of the FS radius with time without assuming values for n or s. Instead, we estimate q from a local fit of the FS R(t) inferred from radio data following Equation 17. In our equations we leave q as a variable, and note that Equation C10 reduces to Equation C9 for the special case q = 1.

D. RADIO DATA AND CALCULATIONS
We provide the full radio dataset of SN 2004C in Table 3. Similarly, we provide in Tables 4 and 5 the physical parameters calculated using the best fit parameters inferred from the modeling of the 28 SEDs modeled with BPL1 (Table 1) and Equations 16 -24. These results are illustrated in Figures 3, 4, and 5. Figure 10 shows the best fit model BPL1 applied to the full data set.   Table 3). Solid lines: our favored model BPL1, a broken power law with smoothing parameter s fixed to −1 and optically thick slope α2 fixed to 5/2. Hollow points: the asymptotic intercept of the thick and thin slopes of the broken power law, identified as F br and ν br in Table 1 and used in our calculations of the physical parameters. At later times, the SED peak migrates to frequencies too low to capture both sides of the SED. In this case, a single power law is fit to the SED and the peak becomes a flux density (frequency) lower limit (upper limit) propagated through the physical parameters.