Introduction

The study of astrophysical transients is growing rapidly through a combination of new instruments, observing strategies, and theoretical advances. Extragalactic transients, such as supernovae (SNe), gamma-ray bursts (GRBs), and tidal disruption events (TDEs), are especially luminous and are typically produced during stellar death. New classes of extragalactic transient continue to be recognized (Inserra 2019).

High-energy and optical telescopes have traditionally dominated the discovery of energetic transients (Gehrels et al. 2004;Bellm et al. 2019). However, radio measurements are emerging as a valuable platform for transient discovery because radio wavelengths are sensitive to shocks formed by fast ejecta that may be expelled in such events (Chevalier 1982a;Mészáros & Rees 1997;Frail et al. 2001a). Radio synchrotron emission formed in shocks has a luminosity that is proportional to the total kinetic energy (Frail et al. 2005;Metzger et al. 2015). This fact has motivated a new generation of radio telescopes and surveys designed to be sensitive to transients (e.g., Murphy et al. 2013;Fender et al. 2016;Shimwell et al. 2017;Lacy et al. 2020).

Early efforts to search for extragalactic radio transients were limited by lack of sensitivity or sky coverage (Levinson et al. 2002;Gal-Yam et al. 2006;Croft et al. 2010;Thyagarajan et al. 2011;Bannister et al. 2011;Mooley et al. 2013;Bell et al. 2015).

New surveys can robustly detect sources brighter than ∼1 mJy over ten thousand square degrees. At this scale, the surveys are sensitive to radio spectral luminosities of L ν = 10 30 erg s -1 Hz -1 over a volume of ∼1 Gpc 3 , and can discover GRBs, TDEs, and more (Metzger et al. 2015). Other factors that have traditionally limited radio transient discovery are issues related to correlator software (specifically for the legacy Very Large Array, e.g., phase-center noise), source significance statistics and lack of supporting multiwavelength measurements (e.g., Gal-Yam et al. 2006;Thyagarajan et al. 2011;Frail et al. 2012). These issues have been relieved though improvements to radio software/data analysis pipelines and better integration with follow-up observing resources (Mooley et al. 2016(Mooley et al. , 2019;;Driessen et al. 2020;Pintaldi et al. 2021).

FIRST J141918.9+394036 (hereafter J1419+3940; Ofek 2017; Law et al. 2018) is an example showing the potential for radio discovery of extragalactic transients. Law et al. (2018) identified the source to be bright in the VLA Faint Images of the Radio Sky (FIRST) survey (Becker et al. 1995) in 1993, undetected by the VLA Sky Survey (VLASS; Lacy et al. 2020) in 2017, and a sub-milliJansky radio source in archival data steadily declining in flux density between 2010 and 2018. The transient was associated with a host galaxy, SDSS J141918.80+394035.9, at z = 0.01957 (Ahn et al. 2012), implying a radio spectral luminosityfoot_0 of at least 2 × 10 29 erg s -1 Hz -1 . This makes J1419+3940 more luminous and longer lived than most SNe, including those associated with long gamma-ray bursts (LGRBs; Corsi et al. 2016). The volume over which this source could have been detected is small (out to ∼100 Mpc), which implies a volumetric rate that is in tension with many of the known radio transient populations. J1419+3940 is luminous, nearby, and at least three decades old, which makes it either a highly fortunate discovery or a prototype of a class of transient not well probed by past radio surveys.

No gamma ray, X-ray or optical counterparts of J1419+3940 have been found, and multiple origin models have been proposed. The luminosity, timescale, and host galaxy are consistent with an afterglow of a LGRB (Levinson et al. 2002). If so, the explosion occurred around 1993 with a total energy E j ∼ 10 51 erg that is interacting with a density of n ∼ 10 cm -3 (Law et al. 2018) of circum-burst medium (CSM). The radio evolution and lack of a gamma-ray counterpart suggests that the event was an off-axis LGRB (also known as an orphan afterglow; e.g., Ghirlanda et al. 2014), the first of its kind. Very-long-baseline interferometry (VLBI) observations measured an expansion speed of 0.1c, which is consistent with that hypothesis (Marcote et al. 2019). Lee et al. (2020) proposed that J1419+3940 could be the sign of interaction between ejecta from a neutron star merger and its surrounding interstellar material (Nakar & Piran 2011). This model can explain the early light curve shape and is more consistent with the high volumetric rate implied by J1419+3940.

Other models for J1419+3940 include new classes of transient powered by central engines. High-cadence optical surveys have defined a new class of engine-driven transients, akin to the prototype AT2018cow. AT2018cow-like events are a subclass of fast blue optical transients (FBOTs; Drout et al. 2013) having luminous radio emission. Radio observations of AT2018cow-like transients have shown synchrotron emission from mildly relativistic outflows (Ho et al. 2019(Ho et al. , 2020;;Margutti et al. 2019;Coppejans et al. 2020). Potentially related is another new class of radio transient hypothesized to be associated with newborn magnetars (Metzger et al. 2015;Murase et al. 2016). This class is potentially frequent (Prajs et al. 2017) and radio luminous, but difficult to identify due to their long evolution timescale (Margalit & Metzger 2018; but also see Ofek 2017). Magnetar engines for transients are especially interesting, because they may leave a magnetar remnant long after the supernova (SN). If so, this could tie the events to millisecond transients such as fast radio bursts or highly luminous off-nuclear radio sources (Ofek 2017; Law et al. 2019;Eftekhari et al. 2020).

Here, we present new radio (including VLBI) and optical observations of transient J1419+3940, and a detailed interpretation of the transient source nature. New VLA and Low Frequency Array (LOFAR) data define a quasi-simultaneous radio spectrum from 0.15-10 GHz in two epochs and extend the time baseline of measurements at 1.4 GHz to 26 yr. We also describe Very Long Baseline Array (VLBA) observations and reprocessing of previous VLA and European VLBI Network (EVN) data sets (Section 2). These measurements allow new analysis of the synchrotron blastwave energetics and a comparison to the direct measure of the shock expansion measured through VLBI observations (Section 3). An analysis of the host galaxy and comparison with the host galaxies of known transients is presented in Section 4. The properties of the radio transient and the host galaxy together with the event rate support the initial interpretation from Law et al. (2018) that J1419+3940 is likely associated an with off-axis LGRB. Section 5 gives possible explanations of a broad emission feature observed in the optical spectrum of J1419+3940. We present a summary and discussion of results in Section 6, and end with the conclusions in Section 7.

Observations, Data Processing and Initial Analysis

All the radio data (new and revised from previously published studies) used in this work are tabulated in Table 1.

VLA

We carried out observations with the National Science Foundation's Karl G. Jansky Very Large Array (VLA; under project code 19A-393; P.I.: Law), on 2019 May 18. Standard 8-bit wideband interferometric digital architecture (WIDAR) correlator setups were used for the P (300-500 MHz), L (1-2 GHz) and S (2-4 GHz) bands and 3-bit setups for the C (4-8 GHz) and X (8-12 GHz) bands to obtain the full frequency coverage possible with the VLA up to 12 GHz. 3C48 and PKS J0118-2141 were used as the flux and phase calibrators, respectively. The data were processed using the NRAO Common Astronomy Software Applications (CASA) pipeline (default version for CASA 5.6.1), and each band was split into an independent measurement set using CASA split. Each measurement set was then imaged using CASA tclean with natural weighting, pixel sizes chosen so as to resolve the synthesized beam with 4 pixels, image sizes appropriate to cover the primary beam full-width at half-maximum (FWHM), and a CLEAN stopping threshold of 3× the thermal noise (as estimated from the VLA exposure time calculator).

VLBA

We observed J1419+3940 with the VLBA at 2.3 GHz on 2019 April 12 (project BL266). All ten VLBA antennas were used. A standard continuum observing mode was used, with eight spectral windows (each of 32 MHz width), and the observations were conducted in a phase-referenced manner with ICRF J141946.6+382148 used as the phase reference calibrator (angular separation of 1.3 degrees). The phase reference cycle consisted of alternating scans of duration approximately 4 min. on the science target J1419+3940 and approximately 45 s on the phase-reference calibrator.

The data were calibrated in two different manners as a consistency check, first using the rPICARD pipeline (Janssen et al. 2019;CASA v. 5.5), second using AIPS (v. 31DEC19). Because of considerable radio frequency interference (RFI), only approximately half of the total bandwidth of 256 MHz was suitable for use.

We made an image using Briggs weighting in CASA and robust = 0.5, which had image background rms values of 31 μJy beam -1 and a synthesized beam of 5.9 mas × 3.1 mas at 12°position angle. J1419+3940 appeared largely unresolved in the image. To get a more precise idea of the source size, we fitted a single circular Gaussian component directly to the visibilities by weighted least-squares using the AIPS task OMFIT. The best-fit Gaussian had a flux density of 415 ± 75 μJy, where the uncertainty includes an assumed 15% uncertainty on the flux-density scale.

For marginally resolved sources, the source size can be significantly correlated with any residual antenna amplitude miscalibration. We therefore determined the uncertainty on the FWHM size by including two components, added in quadrature: first, the statistical contribution from the fit, and, second, the scatter obtained from a small Monte Carlo trial where the antenna gains were randomized by 10% and the resulting visibility data refitted. In this case the statistical component dominated. The measurement suggests a best-fit FWHM size of 1.5 mas, but a completely unresolved source is excluded only at the 1.7σ level. The size measurement can also be compared with the FWHM value of 3.9 ± 0.7 mas (Marcote et al. 2019) obtained with the EVN. In Section 2.5.3 we carry out an independent analysis of that EVN data.

However, for the radio source we are interested in an outer radius. For a circular Gaussian model the formal outer radius is infinity, so the Gaussian FWHM itself is not am appropriate estimator of the radio source size. For marginally resolved sources, the fit is only very weakly dependent on the choice of model. The outer diameters of more physically appropriate models are related to the Gaussian FWHM as follows: uniform disk 1.60 × FWHM, optically thin shell 1.81 × FWHM. Motivated by SNe, we consider an optically thin, uniform spherical shell model where the shell thickness is 20% of the outer radius (see Bietenholz et al. 2021b, and discussion therein).

We found the best-fit outer diameter of such a shell to be 2.5 mas (and, again, a completely unresolved source is excluded only at the 1.7σ level). The detailed model fit results are given in Table 2.

LOFAR/LoTSS

Observations were carried out with the LOFAR (van Haarlem et al. 2013) on 2015 July 28 (P214+40) and 2019 April 13 (P213+37) and as part of the ongoing 120-168 MHz LOFAR Two-meter Sky Survey (LoTSS; Shimwell et al. 2019Shimwell et al. , 2017)). The observations were processed following the current standard imaging procedures as described by Tasse et al. (2021). The data were first calibrated to remove directionindependent effects (van Weeren et al. 2016; Williams et al. To remove the remaining severe ionospheric and beam model errors, the data were then calibrated using the directiondependent self-calibration pipeline DDF pipelinefoot_3 that uses kMS (Tasse 2014;Smirnov & Tasse 2015) to derive direction- , except epochs 1993.87, 1994.31, 1995.32, and 2008.54 Note. The fitted source sizes are given in mas for the VLBA and EVN data sets, with 1σ uncertainties. All values were determined by fitting geometrical models directly to the visibilities by least-squares. The uncertainties include the statistical component and a systematic one derived by allowing for 10% uncertainty in the amplitude calibration of the individual antennas. In the case of the EVN, where the range of data weights is high, we used the square root of the data weights in the fitting, which improves convergence at the expense of a small loss of statistical efficiency. At an angular-diameter distance of 85 Mpc, 1 mas = 0.41 pc = 1.3 × 10 16 cm.

dependent calibration solutions and DDFacet to apply these while imaging (Tasse et al. 2018). The flux scale of the final images was refined using the procedure outlined by Hardcastle et al. (2021). At a FWHM resolution of 6″, the final images have background rms noise levels of 260 μJy beam -1 and 80 μJy beam -1 at the pointing centers of P214+40 and P213 +37, respectively, where the large discrepancy in background rms levels is due to unusually poor conditions during the P214 +40 observation.

Keck/LRIS

The host galaxy of J1419+3940 was observed with the Low-Resolution Imaging Spectrometer (LRIS; Oke et al. 1995) at Keck Observatory on 2019 April 5 (P.I.: Hallinnan). The blueside spectrum (1 × 1200 s) was obtained using the 400/3400 grism and the red-side spectrum (2 × 550 sec) was obtained using the 400/8500 grating. Observations were reduced using LPIPE (Perley 2019). Spectra were flux calibrated using an observation of Feige 34 and the absolute scaling was adjusted to match an archival spectrum (taken in 2004) of the galaxy from the Sloan Digital Sky Survey (SDSS; also shown in Figure 1). The red-side CCD was incorrectly windowed during the observation, producing a small gap in the wavelength coverage.

The reduced spectrum, shown in Figure 1, is that of a strongly star-forming, metal-poor galaxy (Section 4). However, a single broad feature is also evident underlying the narrow [O III] lines. If interpreted as a broad component of [O III], the inferred velocity is v ∼ 3000 km s -1 , orders of magnitude in excess of the escape velocity of the low-mass host galaxy. The origin of this feature is currently unclear. No similar components are seen under any other narrow lines (or elsewhere in the spectrum).

A zoom-in on the [O III] and Hα lines is shown in Figure 2. We estimate the broad component [O III] flux to be about 5 × 10 -16 erg cm -2 s -1 , corresponding to a luminosity of 5 × 10 38 erg s -1 , and a 3σ upper limit of 1 × 10 38 erg s -1 for any broad Hα component. The SDSS spectrum also shows some hints of a broad [O III] component, but it is difficult to ascertain its significance due to the larger noise and coarser resolution compared to the Keck spectrum. Nevertheless, we find that if the broad component is present in 2004 then its luminosity is significantly lower than that measured in 2019.

Reprocessing of Archival Data

There are significant discrepancies in the radio spectral index and source size found between the results we present above and those previously reported in the literature (Law et al. 2018;Marcote et al. 2019). In order to investigate the potential sources of these discrepancies and to obtain improved estimates of the fitted values and/or uncertainties if possible, we reprocessed the archival data, and present the results below.

VLA/FIRST

VLA observations were performed on 1994 August 14, 19, and 20 as part of the FIRST survey (Becker et al. 1995), in two adjacent frequency bands centered on 1364.9 and 1435.1 MHz, with a bandwidth of 21.9 MHz in each band. Observations were calibrated in MIRIAD (Sault et al. 1995); data were flagged for RFI, calibrated for flux using observations of 3C 286, and had their time-dependent antenna gains estimated using observations of 3C 286 and QSO B1504+377 (observed around an hour before and after the target field, respectively). The calibrated UVFITS data were then imported into CASA. Each epoch and each channel was independently imaged using clean using a pixel size of 1 arcsec, an image size of 4000 pixels, Briggs weighting with the CASA robust parameter set to 0.5, a single Taylor term and a CLEAN stopping threshold in the range 0.7-0.9 mJy (roughly 4-5 times the thermal noise). The three maps (each with a slightly different pointing center) for each frequency were then combined in the image plane using AIPS FLATN (with the appropriate primary beam correction defined in AIPS PBCOR).

Figure 1. Spectroscopy of the host galaxy SDSS J141918.80+394035.9. The SDSS spectrum (from 2004) is shown in gray and the Keck spectrum (from 2019) in black. The upper and lower panels show different y-axis scalings of the data (spectra in the upper panel have also been convolved to a resolution of 12 Å). The host is an intensely starbursting galaxy. A broad component is visible under the [O III]λ4959,5007 line.

VLA/15A-033

In order to verify the spectral index of the optically thin part of the radio spectrum at epoch 2015.36, we reprocessed the VLA/15A-033 data set (P.I.: J. Farnes; galactic magnetic fields project), observed in the 1.5 GHz and 3 GHz bands. The raw data were put through the NRAO pipeline built into CASA 5.6.1. The processed data were then split into three frequency bins using CASA split and imaged using CASA clean using a suitable pixel size to sample the synthesized beam with 4 pixels, an image size suitable for imaging the FWHM primary beam, Briggs weighting with robust 0.5, two Taylor terms and a CLEAN stopping threshold roughly 3× the thermal noise.

EVN

As described by Marcote et al. (2019), EVN observations of J1419+3940 were performed in September 2018 under project code RM015. We downloaded the correlator data products from the EVN archive and reprocessed them using a ParselTongue pipeline (Kettenis et al. 2006), which was adapted from that described by Mooley et al. (2018a) with an additional step using the APCAL task to load the a priori amplitude calibration corrections for EVN data.

.We edited data during time periods affected by RFI, and for our final processing we also deleted all the data from the following stations: the Sardinia radio telescope (SRT), Cambridge, Deffin, and Knockin. The SRT solutions displayed a high residual phase rate, while the other three telescopes exhibited phase rate discontinuities. The SRT issue is thought to arise from a position error, while the discontinuities affecting the other three telescopes are believed to result from fiber delay corrections introduced by the WIDAR correlator (B. Marcote 2021, private communication) However, we found that the inclusion or exclusion of these four antennas did not substantially bias the resulting size constraints. Regardless of this issue, the sparse uv sampling of the data set leads to challenges. The longest baselines are primarily to just two stations: Tianma and Hartebeesthoek, and there are few baselines of intermediate length, meaning that, as noted by Marcote et al. (2019), the gain calibration for these two stations can considerably affect the fitted size.

After editing and calibration, we imaged the data using AIPS IMAGR (robust = 0; uvtaper 40 Mλ) and found the fitted synthesized beam to be 5.50 × 4.86 mas at a position angle of 71°(the synthesized beam represents just the narrow inner lobe; sidelobes up to ∼80% in amplitude are present across the broad plateau caused by the abundant shorter baselines). The source appears to be marginally resolved in the image plane, but we note that the negative extremum in the image is at 50% the peak brightness, which is 379 μJy beam -1 (rms noise is 33 μJy beam -1 ).

We noted that the calibrated EVN data have a large range of data weights, with a small fraction of the baselines, in particular those to a single antenna, Effelsberg, having much higher weights than the remainder. To reduce the dominance of this small fraction of the baselines, we used the square root of the original data weights for all baselines in the uv-plane model fittingfoot_4 . Next, following a similar procedure adopted for the VLBA data (Section 2.2), we fit both a spherical shell model and a Gaussian model in the uv plane (using AIPS OMFIT) to find the source size. 25 We find the best-fit values for the shell and Gaussian models in the uv plane are 1.8 mas (outer radius, corresponding to a source diameter of 3.6 mas) and 2.3 mas (FWHM), respectively, as given in Table 2. These measurements are consistent, given the Gaussian FWHM to shell radius conversion factor of 0.8 (Section 2.2). However, we find that these nominal best-fit values have large uncertainties (∼50%-100%, asymmetrical error bars), especially when compared with the relatively precise Gaussian FWHM measurement reported by Marcote et al. (2019) (3.9 ± 0.7 mas; further discussed below).

In an attempt to reproduce the results of Marcote et al. (2019), we undertook a model-fitting procedure similar to the one used in that paperfoot_6 , again fitting a circular Gaussian directly to the visibilities by least-squares, but using Difmap (Shepherd 1997) rather than AIPS OMFIT, and not taking the square root of the nominal data weights. Using Difmap, we obtained a best-fit FWHM of 2.8 mas. As with OMFIT, we found the uncertainty range asymmetric, with a larger uncertainty toward larger sizes. Formally the Difmap 2σ range was 2.4-4.4 mas, approximately consistent with the value of 3.9 ± 0.7 mas published by Marcote et al. (2019) but smaller than the corresponding range we obtained using OMFIT. The non-Gaussian image-plane errors (with, as noted, the negative image extremum being at <-10σ) strongly suggest that there are residual calibration errors, which are probably not Gaussian distributed, and, being antenna-based, would introduce correlations between the visibility measurements. Since the least-squares fits assumed that the errors in the visibility measurements are Gaussian distributed and independent, this could be the cause of the discrepant uncertainty ranges between OMFIT and Difmap. In any case, given the low signal-to-noise ratio and the likely presence of non-Gaussian-distributed errors, we consider that the larger uncertainty range obtained from OMFIT is likely more realistic for the size measurement from the EVN data, and we use these results (as reported in Table 2) from this point onwards.

Modeling

In this section we estimate the physical parameters based on the observational data. In view of the VLBI measurements presented in the previous section, we believe that the most appropriate measurement to use is the outer radius from the shell model. Considering the corresponding VLBA and EVN values listed in Table 2, we take the weighted mean to find the resulting outer radius of - + 1.3 0.6 0.3 mas. 27 This corresponds to a physical radius of

at an angular-diameter distance of 85 Mpc.

Power-law Fits to the Radio Light Curve and Spectra

We fit the data points from early times (obtained in 1993-95) and late-time data points at ∼1.4 GHz with power laws F ∝ t b . We do not know the time of explosion, so we try two fiducial values, 100 days and 1000 days (motivated by the arguments presented by Law et al. 2018), for the age of the transient at the time of the first radio detection at epoch 1993.87. Assuming age 100 days (1000 days) postexplosion, we find F ν ∝ t -0.3 and F ν ∝ t -2.3 (F ν ∝ t -1.0 and F ν ∝ t -2.7 ; Figure 3) for the early and late time, respectively.

Simple power-law fits, with F ∝ ν α, where α is the spectral index, to the optically thin spectra at epochs 2015.36 and 2019.46 give values of α of -0.8 ± 0.2 and -1.06 ± 0.10, respectively. The optically thin spectral index is therefore consistent with α thin = -1, suggesting that the electron power-law index is p ≈ 3 (where we have assumed that the GHz spectrum lies between the synchrotron self-absorption and cooling frequencies). Fixing this spectral index and using α thick = +2.5, we fit the epoch 2019.5 radio spectrum with a smoothly broken power law (SBPL) 28 of the form described by Beuermann et al. (1999) and Mooley et al. (2018b), using Markov Chain Monte Carlo to find 29 the peak flux density, 2.25 ± 0.52 mJy, the peak frequency, 0.30 ± 0.04 GHz, and the smoothness parameter, -+ 0.39 0.33 0.45 . The radio spectral evolution with these fits is shown in Figure 3.

Spectral Evolution and Synchrotron Self-absorption (SSA) Analysis

We use the spectral parameters derived above to calculate the radius (R), magnetic field (B), and energy (U) using the Chevalier (1998) prescription (see Ho et al. 2019, for the relevant equations using p ≈ 3). These calculated values are tabulated in Table 3. For the 1993.87 epoch we consider for demonstrative purposes that the SSA frequency (ν a ) is 1.4 GHz. The average velocities implied by the equipartition (ò e = ò B = 1/3) radii estimated at epochs 27 Although the associated uncertainty is large, the size measurement cannot immediately be dismissed as an upper limit since the source appears to be at least marginally resolved with the EVN. Further observations will be needed to improve the precision on this measurement. We also note that, since the measured size cannot be below zero, any measurement of the size will be biased high in the case of low SNR (like Ricean bias). For both EVN and VLBA, the best fit is about 2σ above zero, so this bias may be significant and there will be a small upward bias in the weighted mean. Quantifying this bias is, however, nontrivial. 28 Fitting a broken power law, which corresponds to a SBPL with smoothness parameter  s inf , we get a peak flux density S p = 3.4 ± 0.3 mJy and a peak frequency 0.28 ± 0.03 GHz. In this case, the parameter estimates given in Section 3.2 change as n µ

1 compared with the SBPL case given in Table 3. Specifically, R increases by a factor of 1.3 compared to the SBPL case. 29 Since the turnover frequency depends on the flux densities of J1419+3940 at 360 MHz and 150 MHz (obtained using different instruments), we verified the spectra of two nearby radio AGNs, FIRST J141849.5+395154 and FIRST J141828.5+393928. Their spectra between 150 MHz, 360 MHz and 1.4 GHz appear perfectly consistent with a single power law with spectral index -0.5, indicating that there are no additional systematic offsets between the LOFAR and VLA P-band data points.

1993.87 and 2019.5 are about 44,000 (t d /1000 d) -1 km s -1 and 7000 (t/26 yr) -1 km s -1 , respectively (the former value is a lower limit since the peak luminosity at 1.4 GHz may be higher than that observed in 1993; t d is the age of the transient at the discovery epoch 1993.87 and t denotes the age around the VLBI observing epoch 2018/19). In comparison, the average velocity as implied by the VLBI radius measurement of ; 1.2 mas is aboutfoot_7 19,000 (t/26 yr) -1 km s -1 . These velocities are reminiscent of Type Ib/c and Type Ic-broad line (BL) SNe and make a Type II SN explanation unlikely.

We also note that the cooling frequency, ν c , is far above our observing band, and therefore irrelevant to this analysis. Specifically, from Sari et al. (1998) we estimate the cooFing frequency at epoch 2019.5 to be ν c ∼ 700 GHz for the ò e = ò B = 1/3 case (Table 3) and even higher for lower values of ò

In Figure 4 we compare the peak luminosity and timescale of J1419+3940 with those of different classes of stellar explosion, including fast blue optical transients (FBOTs). This figure also places the velocities derived above in the context of different radio afterglows, again indicating that a Type II SN explanation is disfavored (see also Bietenholz et al. 2021a).

Energetics and Light Curve Modeling

Given the measured luminosity, spectrum, and radius of J1419+3940, we can place direct constraints on underlying properties of the source. Typically, in modeling synchrotron blastwaves there is a degeneracy between blastwave energy and ambient-medium density. In the following the source-size measurement can be used to break this degeneracy and unambiguously constrain the source energetics.

First, we explicitly define the formalism. We assume that the observed radiation at GHz frequencies is produced by optically thin synchrotron emission from a nonthermal population of electrons for which the momentum, γβ, is distributed as a power law, gb µ -p ( ) . The optically thin synchrotron spectrum is n ~-p 1 2 which, given the observed spectral index α ≈ -1, implies that p ≈ 3 (see Section 3.1). We further consider the standard scenario in which the synchrotron-emitting electrons are accelerated at a nonrelativistic shock (with efficiency ò e ), where magnetic fields are also amplified (with efficiency ò B ). For p = 3, the synchrotron luminosity of such a blastwave can be expressed as

10 erg s 26 yr 4.8 0.5 pc 0.4 ; ISM, 8.9 0.5 pc ; wind,

where we have separated into two cases depending on whether the ambient medium into which the shock expands has a constant number density n (ISM; normalized to n 0 ≡ n/1 cm -3 ) or a wind-like density profile ρ = Ar -2 (wind; which we normalize to A å ≡ A/5 × 10 11 g cm -1 ). We have used the notation q x = (q/10 x ) in the appropriate unit for parameter q, e.g., ò B,-1 ≡ (ò B /0.1).

such that m  1 is a correction factor to the "average" velocity R/t (in the case of a power-law temporal evolution of the shock front, m describes this exponent, i.e., R ∝ t m ). In the context of radio SNe, where the ambient medium is a wind environment, Chevalier (1982b) shows that m = (k -3)/(k -2), where ρ ej ∝ r -k is the outer density profile of the SN ejecta, and k  7 are typical values (implying m  0.8). On the other hand, a blastwave that propagates into a constant-density ISM and that is deep within the Sedov-Taylor regime will be characterized by m = 0.4. Finally, Equation (1) has been derived assuming that the blastwave is in the so-called deep-Newtonian regime discussed by Sironi & Giannios (2013). This regime is relevant if

(assuming p = 3), and is therefore appropriate for J1419+3940 at the current epoch (Equation ( 2)).

Using Equation (1) we can find the ambient density that is required in order to produce the observed spectral luminosity (at 1.5 GHz at epochfoot_8 2019.46) of J1419+3940, νL ν ; 9 × 10 36 erg s -1 . It is

in the ISM case, and

for a wind medium.

Assuming that the blastwave in the ISM case is currently in the Sedov-Taylor regime (otherwise the light curve would be rising rather than declining; e.g., Nakar & Piran 2011), the energy associated with this blastwave is

where ξ ; 1.15 is given by the Sedov-Taylor solution. Note that Equation (5) does not depend on the assumed source age. This energy is reasonable for various astrophysical sources, and in particular for LGRBs.

In the radio-SN case, Chevalier (1982b) showed that the shock radius is

where U c is a parameter governing the outer ejecta density profile, r = --

. This parameter can be related to the ejecta energy that is contained above some velocity coordinate,

) . Using these expressions, we can estimate a lower limit on the ejecta energy that is required for interpreting J1419+3940 within the radio-SN paradigm:

The second line is for k ; 6.3 (and correspondingly m ; 0.77), at which E ej attains a minimum as a function of k (this accounts for the m = (k -3)/(k -2) dependence, and the additional m dependence implied by Equation (4)). Different values of k would imply larger ejecta energies. For example, for the physically motivated value of k ≈ 12 (Matzner & McKee 1999), we find that E ej ∼ 10 52 erg. Furthermore, Equation (6) only accounts for the energy of ejecta material that has velocity v  20, 000 km s -1 (as inferred in Equation (2)) so that the total ejecta energy (including lower-velocity material) is likely to be even higher. This energy constraint disfavors the interpretation of J1419+3940 as a typical 32 radio SN (see also Ofek 2017).

Initial Velocity

Having ruled out the typical radio-SN scenario (and hence wind medium) on energetic grounds, we from here focus on the constant-density (ISM) blastwave scenario. As discussed above, the fact that J1419+3940's light curve is observed to decline implies that the shock is within the Sedov-Taylor regime, t > t dec , i.e., the shock must be decelerating. This implies that the shock velocity at the current epoch (Equation ( 2)) is lower than the initial blastwave velocity, v i . We can roughly constrain this initial velocity by considering the deceleration (or Sedov-Taylor) timescale, t dec (e.g., Hotokezaka & Piran 2015):

The deceleration time depends on the initial blastwave velocity, v i (where

2 is the corresponding Lorentz factor), and we specifically consider also the possibility that the initial blastwave was highly relativistic (bottom case). For a spherical explosion (and neglecting synchrotron selfabsorption), the light curve rises up to t ∼ t dec and subsequently declines. Equation (7) shows that the light curve rise time is ∼decades for outflows whose initial velocity is subrelativistic. In particular, t dec would be  16 yr if the initial velocity, v i , were similar to the current inferred velocity (Equation ( 2)). This is in tension with the data, which instead suggest a much faster light curve rise time. The fact that the observed flux density of J1419 +3940 declines by a factor of ∼2 in the span of ≈1.5 yr (between epochs 1993.87 and 1995.3) since the first detection epoch is suggestive of a rise timescale 2 years. Using Equations (3, 5, 7) we can use this as a constraint on the initial velocity. Requiring that t dec  2 yr implies that the initial velocity must have been

at the very least transrelativistic. Another (albeit not mutually exclusive) possibility is deviation from spherical symmetry, such as an initially collimated explosion pointed off-axis from our line of sight. The short rise time in this case may be attributed to a jet break; however, this scenario would also require an initially relativistic outflow (for the case of nonrelativistic outflows, asymmetry has only a modest effect; e.g., Margalit & Piran 2015).

Synchrotron Self-absorption (SSA)

As a final point, we can estimate the SSA frequency implied by the ISM-case solution for J1419+3940 (Equation ( 3)). We estimate ν SSA by equating the optically thin synchrotron luminosity (Equation ( 1)) to the optically thick luminosity

) and B is the magnetic field. This approximate approach is correct to order-unity correction factors due to the (uncertain) geometry and the electron distribution function. In this manner, we find that (for p = 3) n »

´-

For the inferred ISM density of J1419+3940 (Equation (3)), this implies

which depends almost exclusively on the shock radius R. This result is broadly consistent with the observed spectrum (at this epoch, ∼2019), which exhibits a turnover at a few hundred megahertz, which is compatible with SSA (Figure 3; see also Section 3.1). Note that precise details of the transition between the optically thick and optically thin regime depend on additional geometric effects, which we do not consider here (e.g., Björnsson & Keshavarzi 2017).

Host Galaxy

The host galaxy of the transient, SDSS J141918.80 +394035.9, is a blue compact dwarf and is clearly detected in SDSS survey imaging (see Law et al. 2018). The spectra are shown in Figure 1. We extracted line fluxes by fitting a Gaussian profile to each emission line, and obtained estimates of the gas-phase oxygen abundance using both strong-line methods (all 13 diagnostics enabled within pymcz; Bianco et al. 2015) as well as a T e -based "direct" measurement based on the [O III]λ4363 auroral line (Izotov et al. 2005). Other bulk measurements (including mass and star formation rate (SFR) estimates) were obtained from the NASA-SDSS Atlas (NSA; Blanton et al. 2011). We infer a stellar mass of M * = 3 × 10 7 M e , specific SFR (sSFR) of 2.4 × 10 -9 yr -1 , and an oxygen abundance of 12+log[O/H] = -+ 8.10 0.05 0.06 (T e method; strongline methods give consistent estimates.)

The left panel of Figure 5 shows the spectral energy distribution (SED)-derived stellar mass (M * ) and specific star formation rate (SFR/M * ) of the host galaxy of J1419+3940 in comparison to the SDSS spectroscopic sample and some other populations known to reside within extreme galaxies: SLSNe and GRBs (Perley et al. 2016), Ic-BL SNe (Modjaz et al. 2020a), as well as a 1/V max resampling of the NSA spectroscopic sample. J1419+3940 is at the extreme of all of these groups, but it is not an outlier. Similarly, the right panel of Figure 5 plots two strong-line emission ratios (as in Baldwin et al. 1981); the host lies near one end of the diagram due to its low metallicity but its ionization properties are otherwise consistent with a normal star-forming dwarf.

While it is impossible to come to a secure concPlusion about the nature of the progenitors of J1419+3940-like transients based on a single event, these properties derived above generally support a massive stellar origin and favor a progenitor that is intrinsically more likely to form at low mass, high sSFR, or low metallicity. GRBs, Ic-BL SNe, superluminous SNe, and AT2018cow-like fast transients all seem to have these properties (Japelj et al. 2016;Modjaz et al. 2020b;Schulze et al. 2018;Perley et al. 2021). An exotic origin is otherwise not required, as while the properties of the host are not typical they are far from unprecedented (approximately 1 in 50 SNe explodes in a similarly extreme environment; Taggart & Perley 2021).

Origin of the Broad [O III]λ4959,5007 Line

Here we consider several possibilities for the broad (∼3000 km s -1 ) [O III] line observed in the Keck/LRIS spectrum (Section 2.4).

(a) The broad component could originate from collisional broadening in extremely dense star-forming regions. Broad profiles of similar widths have been seen in some of the most extreme H II regions in the SMC (Testor & Pakull 1985;Kurt et al. 1999), although, to our knowledge, it has not been reported in integrated spectra of entire galaxies. The host galaxy is extremely star forming and moderately metal poor, so this interpretation is quite plausible. A possible challenge to this interpretation is the lack of a similar broad component to the Hα line. (b) It could also originate from a fast outflow driven by star formation. However, the relatively high velocity makes this possibility unlikely (typical galactic winds are of order 100 km s -1 , with even 1000 km s -1 considered to be extreme; Veilleux et al. 2005). (c) An active galactic nucleus (AGN) is another possible explanation. Some dwarf galaxies are known to harbor AGNs (or "proto-AGNs"; III] due to CSM interaction has been observed in decades-old SNe, but in all cases luminous hydrogen emission is also produced (e.g., Milisavljevic et al. 2012). Nebular spectra of hydrogen-poor SLSNe also show hydrogen features (although nebular spectra decades post-explosion have not been reported) (e.g., Yan et al. 2015;Nicholl et al. 2016). Thus SN-CSM interaction cannotfoot_10 explain the spectrum of J1419+3940. Broad nebular emission lacking hydrogen lines years after a SN has been seen in at least one previous case (Milisavljevic et al. 2018), although only on timescales of a few years and not decades, and in that case other oxygen lines were observed that are not apparent here. J1419+3940 occurred prior to the SDSS spectrum shown in Figure 5 being taken, and while (given the lower sensitivity) it is not completely clear whether the line is absent in 2004, it certainly was not any brighter in comparison to 2019, which makes this explanation relatively unlikely.

Given these possibilities, we consider a collisional broadening component to be the most conservative interpretation, but the other possibilities cannot be entirely ruled out. Some discussion of the possible implications if the [O III] line can be attributed to the transient are discussed in Section 6.

Summary and Discussion

Summary of Radio Observations and Derived Spectral Information

We have carried out late-time radio (VLA 1-12 GHz, LOFAR 0.15 GHz, and VLBA 3 GHz) and optical spectroscopic follow-up observations, and reprocessing of archival radio data of the transient FIRST J141918.9+394036 (J1419+3940; Law et al. 2018). For the first time, we unambiguously determine the peak in the afterglow spectrum around 0.3 GHz (at epoch 2019.46, i.e., 26 yr after the first detection of the transient with the VLA at epoch 1993.87). We identify the peak with the SSA frequency and the optically thin part of the spectrum lying between the SSA and cooling frequencies. The optically thin part of the spectrum at epoch 2019.46 satisfies F ν ∝ ν -1 (Figure 3), indicating an electron power-law index of p ; 3. This radio spectral index of α 1-10 GHz ; -1 can be compared with the X-ray-to-radio constraint of α R-X < -0.25 implied by the Swift X-ray flux upper limit (Law et al. 2018). The late-time decline of the 1.4 GHz light curve is consistent with t -2.3 (t -2.7 ) assuming that the first detection occurred 100 days (1000 days) post-explosion. These fits, shown in Figure 3, do not suggest any further steepening of the light curve around epoch 2017/18 (suggested earlier by Law et al. (2018) and Marcote et al. (2019), based on the 3 GHz nondetection in the VLASS Epoch 1 quick-look data and the 1.5 GHz EVN flux density).

Decline of the Late-time Radio Light Curve is Relatively Fast

The inferred light curve decline rate is steeper than that expected ( ∼ t -1.2 ) for a subrelativistic Sedov-Taylor blastwave within the deep-Newtonian regime (Sironi & Giannios 2013). However, this is inferred over a small dynamical range in time (  0.3 dex) and is sensitive to the assumed explosion epoch. We note that a decline rate ∼ t -2.4 , more consistent with the observed light curve, is expected if the minimal Lorentz factor for the population of radio-emitting electrons is ? 1 (e.g., Frail et al. 2000). However, Sironi & Giannios (2013) show that this is only possible if the shock velocity exceeds

, which is not satisfied by J1419+3940 (Equation ( 2)). This velocity threshold could be lowered if only a small fraction ζ e = 1 of electrons that are swept up by the shock participate in diffusive-shock acceleration, as long as the total energy carried by these electrons remains significant (large ò e ). If the decline Figure 5. Left: mass vs. specific SFR for the host galaxy of J1419+3940 and the hosts of a variety of comparison samples: superluminous SNe (blue diamonds), GRBs (purple triangles), B (yellow crosses), FBOTs (cyan pluses) and SDSS galaxies (gray circles). SDSS galaxies have been resampled to simulate a volume-limited survey and the size of the points are scaled by SFR to visually represent a SFR-selected sample. SDSS galaxies with AGN contamination are recolored. Right: BPT line-ratio diagram for the host galaxy of J1419+3940 as compared to the same comparison samples as in the left panel. Solid lines show the AGN separation criterion. The emission properties of the host galaxy of J1419+3940 are typical of low-redshift, star-forming dwarf galaxies.

) and constitute a novel constraint on this parameter (typically implicitly taken to be ζ e = 1). Finally, we note that a steeply declining light curve could point instead toward a drop in the CSM density profile encountered by the shock and/or changing microphysical parameters (time varying ò e , ò B ), though it is unclear whether these scenarios are well motivated in the case of J1419+3940.

Summary of VLBI Results and Parameters Derived from the Radio Analysis

The source radius, measured using VLBA and EVN (data obtained in 2018/19), is measured to be -+ 1.3 0.6 0.3 mas, i.e., = -+ R 0.5 0.2 0.1 pc for an angular-diameter distance of 85 Mpc. This implies an average velocity of 19,000 km s -1 ( ; 0.06c) over the ∼26 yr evolution of J1419+3940 (see Equation ( 2)). This size constraint allows us to constrain the blastwave energy and ambient density of J1419+3940 independently. For a constant-density (e.g., ISM) circumstellar medium, we find n ≈ 40 cm -3 and E ≈ 5 × 10 50 erg (Equations (3), 5). These values also imply a SSA peak at ∼170 MHz. These results are broadly consistent with SSA analysis (Section 3.2, Figure 3), which suggests a blastwave energy 10 49 -10 51 erg, and magnetic field strength ∼10-100 mG.

Summary of Optical Spectroscopic Findings

In the optical we find a broad, ∼3000 km s -1 , [O III] λ4959,5007 line that may be a collisionally excited nebular emission from compact star-forming region(s) in the host galaxy, but we cannot confidently rule out its association with J1419+3940.

Classification of the Transient

No multiwavelength counterparts have been detected for J1419+3940 apart from the emission features described above. Using our radio observations we can estimate the late-time properties of the blastwave, but it is difficult to ascertain the progenitor and whether the transient was initially relativistic. However, we consider the following lines of argument to further investigate the nature of the transient. Given the properties of the host galaxy (Section 4), we primarily consider stellar explosion scenarios. In Figure 6 we compare the light curves of GRBs and SNe with J1419+3940. The peak luminosity of > 2 × 10 29 erg s -1 Hz -1 at 1.4 GHz is unprecedented for regular SN afterglows (Weiler et al. 2002;Bietenholz et al. 2021a), but it is compatible with GRB afterglows (or their associated SNe Ic-BL, e.g., SN 1998bw) and possibly AT2018cow-like events. Acknowledging the caveat that a comparison between the afterglows of optically selected SNe and the radio-selected afterglow of J1419+3940 could be biased, we proceed with the following discussion.

In rare cases dense-CSM interaction, as seen in SN IIn, Ibn and Ic-BL, can produce high radio luminosity. Especially in the case of SN Ic-BL PTF 11qcj (Palliyaguru et al. 2021), the afterglow is currently undergoing rebrightening due to late-time CSM interaction (the 1.4 GHz luminosity is currently about 10 29 erg s -1 Hz -1 at ;3000 days post-explosion; see Figure 6). We cannot immediately rule out CSM interaction based on the properties of the afterglow light curve, but we return to this point below.

In Figures 4 and6 we have shown well-known SN Ic-BL: 98bw, 03lw, 06aj (GRB associated), 09bb, 02ap and 11qcj (not GRB associated). It is evident that J1419+3940 is much more luminous (∼2×-20,000×) and longer lived than normal radioloud Ic-BL (but not rising to late times like the interacting 11qcj, discussed above). Moreover, most SN Ic-BL are not radio detected (e.g., Corsi et al. 2016). Hence, a simple Ic-BL afterglow explanation for J1419+3940 appears to be unlikely (but CSM-interacting or off-axis GRB-associated Ic-BL remains plausible).

We can compare J1419+3940 with AT2018cow-like events (Ho et al. 2019(Ho et al. , 2020;;Margutti et al. 2019;Coppejans et al. 2020), but since only a handful of such events are currently known, their properties remain uncertain and the comparison cannot be conclusive. Although the peak luminosity of J1419 +3940 (at 1.4 GHz) is much larger than that observed for any of the SN2018cow-like events, one such event (CSS161010; Coppejans et al. 2020) has a luminosity approaching 10 29 erg s -1 Hz -1 . However, the radio light curves (around 1.4 GHz) of such FBOTs are generally seen to peak on timescales of few 100 days and decline rapidly (faster than t -3 ; e.g., Ho et al. 2020), unlike the ∼ t -2.5 and decades-long emission seen for J1419+3940.

The sample size of radio-detected SLSNe is even smaller (Eftekhari et al. 2019(Eftekhari et al. , 2021;;Law et al. 2019;Coppejans et al. 2021), and, although their 1.4 GHz radio spectral luminosities are < 10 28 erg s -1 Hz -1 , their association with J1419+3940 cannot be immediately ruled out (but see below).

Only one previous radio-discovered afterglow, SN 1982aa (Mrk 297A; Yin 1994), is known to have a peak spectral luminosity around 10 29 erg s -1 Hz -1 at 1.4 GHz and peak timescale of ∼1000 days. Bietenholz et al. (2021a) suggest, based on the luminosity and timescale, that SN1982aa may be a SN IIn, but the nature of the transient remains uncertain since no optical spectrum is recorded. J1419+3940 differs from SN1982aa in that the late-time decline in the radio light curve (t -2.5 ) and the radio spectrum (ν -1 ) are much steeper than those of SN1982aa (t -1.3 and ν -0.75 ). The properties of SN1982aa are generally in agreement with those measured for radio SNe (Weiler et al. 2002;Yin 1994;Bietenholz et al. 2021a), and, notably, are also similar to SN1998bw, while J1419+3940 appears to be an outlier in this group.

The peak luminosity, timescale and decline of the radio light curve of J1419+3940 are similar to those seen for some GRB afterglows.

LGRBs have peak 1.4 GHz spectral luminosities around 10 30 erg s -1 Hz -1 (short GRBs have lower peak spectral luminosities), peak timescale around ∼few 100 days, and the light curve decline at late times is ∼ t -1t -2 . For offaxis events the peak luminosity and timescale may be longer (as in the case of the neutron star merger GW170817; Dobie et al. 2018;Fong et al. 2019), depending on the observing angle. The late-time evolution of off-axis and on-axis GRB afterglows is expected to be similar (e.g., Granot et al. 2002;Salafia et al. 2016;Kathirgamaraju et al. 2016). Note that if we assume the age of J1419+3940 to be ∼1000 days in 1993, then the light curve of J1419+3940 is strikingly similar to that of GRB030329 (see Figure 6).

Taken together, the light curve of J1419+3940 is dissimilar to the afterglows of optically selected SNe (SN II, Ib/c and Ic-BL; i.e., not associated with jets/central engines) and to the handful of AT2018cow-like events that are currently known, but consistent with GRB afterglows. The connection with SLSN-I cannot be ruled out purely based on light curve arguments, but we return to the case of SLSN-I below.

Velocity

As described in Section 3, an average velocity of at least ∼44,000(t d /1000 day) km s -1 is needed to explain the first radio detection of J1419+3940 in 1993, and if the rise in the light curve is as rapid as the decline, then initial velocity needs to be 0.2 c. The average velocity up to the mean VLBI observing epoch 2019.1 is 0.06c. Such velocities (together with the peak radio luminosity) can be reconciled only by SN Ic-BL or engine-driven explosions/jets and rule out SN II (e.g., Bietenholz et al. 2021a). These velocities are also compatible with the dynamical ejecta from neutron star mergers, as pointed out by Lee et al. (2020), and we return to this point below.

Energetics

In Figure 6, we show blastwave energy as a function of velocity for J1419+3940 and other stellar explosions. The energy and velocities are large compared to most SNe, but comparable to those of GRBs and AT2018cow-like events. Specifically, as also shown in Section 3.3, the energy required to explain the properties of J1419+3940 is very large for typical radio SNe, and makes this scenario unlikely. The energy reservoir of 10 51 erg also makes SN 2009bb-like events (Soderberg et al. 2010;i.e., engine-driven Ic-BLs lacking GRB counterparts) improbable. As noted earlier, a CSM-interacting Ic-BL like 11qcj could still explain the properties of J1419+3940, but given the higher implied energy (and velocity) at the peak of the afterglow light curve (see Figure 7), we disfavor this explanation. We can also consider the case of dynamical ejecta from neutron star mergers (Lee et al. 2020). Ejecta of ∼10 50 erg at speeds ∼0.1-0.5c are expected from the simulations of neutron star mergers, while the outflows from black hole-neutron star mergers could be faster and may reach ∼10 52 erg (e.g., Rosswog et al. 2013). Hence, the dynamical ejecta explanation requires an exceptional circumstance. (e.g., Fong et al. 2015;Troja et al. 2019;Makhathini et al. 2021), while steeper values, p ; 3.0, are more common in mildly/ nonrelativistic blastwaves (with the exception of SN II; e.g., Weiler et al. 2002;Chevalier & Fransson 2006;Ho et al. 2019;Coppejans et al. 2020; but also see Soderberg et al. 2010). The index p ; 3 derived for the late-time afterglow of J1419+3940 may therefore be suggestive of a mildly/nonrelativistic transient, but this does not conclusively rule out an initially relativistic blastwave. For example, the index could change during the relativistic to nonrelativistic transition. Indeed, such a transition has been suggested for the TDE Swift J1644+57, which harbored an initially relativistic jet (Cendes et al. 2021). In the case of GRB 030329, it is found that p ; 2.1-2.5 during the relativistic and nonrelativistic regimes, so the electron power-law index may not have changed appreciably (Frail et al. 2005;van der Horst et al. 2008;Mesler & Pihlström 2013).

Host Galaxy and Local Environment

The host galaxy of the transient is a blue compact dwarf, characterized by a high specific SFR and low metallicity. It is similar to the hosts of LGRBs, SN Ic-BL, SLSNe, and AT2018cow-like FBOTs. Lee et al. (2020) evaluated that ∼1% of neutron star mergers/short GRBs may occur in such hosts. However, the relatively dense CSM (n  10 cm -3 ) needed to explain the afterglow of J1419+3940 demands further fine-tuning for the delay time and makes such an explanation unlikely.

Rates

In order to better understand the ratesfoot_11 of radio transients like J1419+3940, we carried out a flux-density-limited search for transients detected in FIRST (White et al. 1997) and absent in the VLASS Epoch 1.0 (Gordon et al. 2021). We do not find any other radio transient (particularly, having luminosity larger than ∼2 × 10 28 erg s -1 Hz -1 ) that are >4 mJy at 1.4 GHz in FIRST and absent in VLASS (i.e., 3 GHz flux density <1 mJy; in this search we removed candidates that were nuclear and hence very likely to be AGNs or TDEs). Considering this unique event found in the 10,000 deg 2 of FIRST, we calculate the corresponding Poisson 68% confidence interval to be 0.4-2.4 events (Gehrels 1986). This corresponds to > = -

´-4mJy, 1.4 GHz 4 24 10 5 R(

) ( ) deg -2 of the sky. Alternatively, assuming a timescale of ∼10 yr above 4 mJy at 1.4 GHz and a peak luminosity of 10 29 erg s -1 Hz -1 , we can estimatefoot_12 the volumetric rate: 40-240 Gpc -3 yr -1 (68% confidence intervalfoot_13 ; median rate is 100 Gpc -3 yr -1 ). ), is assumed to be 10,000 days post-explosion. Blue circles denote SNe (Weiler et al. 1986;Yin 1994;Kulkarni et al. 1998;Soderberg et al. 2005Soderberg et al. , 2006Soderberg et al. , 2006aSoderberg et al. , 2010;;Salas et al. 2013;Palliyaguru et al. 2021) where the equipartition parameters are calculated at ∼50-1500 days post-explosion (except SN2006aj for which the 5 GHz peak is around 5 days), when the SSA peak lies between 1-5 GHz. Two data points plotted for SN2007bg represent the two peaks observed in its 6 GHz light curve. Grey circles are SNe from the compilation of Bietenholz et al. (2021a), where the SSA peak is taken to be the peak of the light curve (observing frequency between 5-10 GHz). Yellow triangles denote AT2018cow-like events (Ho et al. 2019;Margutti et al. 2019;Ho et al. 2020;Coppejans et al. 2020). GRBs occupy the phase space shown by the gray region. Green squares denote the evolution of the afterglow of GRB030329 (Frail et al. 2005;van der Horst et al. 2008).

This rate corresponds to ∼0.1% of the rate of core-collapse SN (Taylor et al. 2014), ∼1% of SN Ic-BL (Kelly & Kirshner 2012;Graham & Schady 2016;Ho et al. 2020), ∼10% of SN Ic-BL (Graham & Schady 2016;Ho et al. 2020), and is comparable to the rates of SLSN-I (∼30 Gpc -3 yr -1 ), AT2018cow-like events (∼400 Gpc -3 yr -1 ) and estimates for LGRBs (∼60 -3 yr -1 ) (Guetta et al. 2005;Quimby et al. 2013;Goldstein et al. 2016;Ho et al. 2020) in the local universe. The rate is also consistent with that of binary neutron star mergers, but in this case, and similarly for SLSN-I, it requires ∼50%-100% of the mergers/explosions to produce luminous radio emission, which is not consistent with observations (e.g., Fong et al. 2015;Horesh et al. 2016;Eftekhari et al. 2019;Schroeder & Margalit 2020;Makhathini et al. 2021;Abbott et al. 2021). Therefore we can rule out SLSN-I and neutron star merger explanations for J1419+3940.

Putting it All Together

Typical radio SN II/Ib/Ic are ruled out based on energy and velocity. Engine-driven SN Ib/c (2009bb-like) are disfavored on energetic grounds. SN Ic-BL are disfavored based on the peak luminosity and timescale as well as the long-lived radio light curve. Considering the energy and velocity of the blastwave, we believe that SN Ic-BL with CSM interaction (11qcj-like) is unlikely. From the small sample of AT2018cow-like events, this class of transients seems unlikely on account of the shape of the afterglow light curve. Transient rates suggest that SLSN-I and neutron star mergers are unlikely. On the other hand, the afterglow light curve, velocities, blastwave energy, host galaxy properties and transient rates are compatible with the LGRB class. We therefore conclude that, in terms of previously studied transients, the afterglow properties of J1419+3940 are most consistent with those of LGRBs. A further exotic explanation, e.g., involving a magnetar or stellar merger, may not be required (but see below for a short discussion on the magnetar scenario).

Similarity with LGRBs: Inverse Beaming Fraction and Jetopening Angle

We conclude based on the above arguments that an LGRB remains the most likely explanation (Law et al. 2018;Marcote et al. 2019) for J1419+3940. Under this premise, we can calculate the inverse beaming fraction ( q º --

) ) for GRBs, where θ j is the average jet half-opening angle. Using the formalism of Levinson et al. (2002), we calculate this parameter as 37

where N is the number of afterglows detected in our search above the minimum flux-density threshold of 4 mJy,  is the observed rate of GRBs in the local universe (z = 0) beamed toward us, n is the CSM density, R is the measured source radius in pc, t is the age of the transient around epoch 2019.5, and we have used Equation (5). We have used the notation q x = (q/10 x ) for the microphysical parameters, as in Section 3.  degrees (68% confidence), consistent with previous estimates (Frail et al. 2001b;Levinson et al. 2002;Guetta et al. 2005;Gal-Yam et al. 2006;Goldstein et al. 2016).

Predictions for Future Radio Surveys and Future

Evolution of J1419+3940

The rate >4mJy, 1.4 GHz R(

) derived above suggests that a radio survey across a hemisphere with few-milliJansky sensitivity should be able to find ∼1 J1419+3940-like transient. Such surveys are currently being executed with the VLA (the VLASS; Lacy et al. 2020), ASKAP (e.g., McConnell et al. 2020), and an even deeper survey has been proposed for the MeerKAT (Santos et al. 2016). Therefore, we predict that at least a few J1419+3940-like transients will be discovered in the coming years.

Since LGRBs are accompanied by GRB-SNe, it is possible that the late-time radio luminosity of J1419+3940 will enter a second rise phase as the slower-moving SN ejecta collides with the ambient CSM (Barniol Duran & Giannios 2015; Kathirgamaraju et al. 2016;Peters et al. 2019;Margalit & Piran 2020;Eftekhari et al. 2021). The relatively close distance of J1419 +3940 combined with its old age make this source an opportune target for detecting such an emission, which would present unique possibilities for probing the additional physics of the explosion (Margalit & Piran 2020). We therefore recommend continued radio monitoring of J1419+3940.

An Alternative Explanation for J1419+3940 Involving a Magnetar

Finally, we consider an alternative (speculative) possibility that J1419+3940 arose from a SN (that may or may not be associated with a LGRB) that gave birth to a long-lived central engine, such as a millisecond pulsar or magnetar. In particular, a long-lived magnetar could power a nebula of synchrotron radio emission (e.g., Murase et al. 2016;Margalit & Metzger 2018), which would become visible in radio once the SN ejecta shell becomes optically thin to free-free absorption. Although the ejecta in the case of SLSNe may take decades or longer to become optically thin at GHz frequencies (Margalit et al. 2018), inconsistent with the rapid early light curve decay of J1419+3940, this transition could happen sooner for an explosion with a low ejecta mass (e.g., similar to those inferred in FBOTs or ultrastripped SNe). One motivation for this scenario is the speculation that the [O III] emission line observed at late times from J1419+3940 could arise from the nebular phase of an engine-powered SN, as a result of UV/X-ray emission from the engine nebula being reprocessed by the ejecta shell and cooling through the 37 We have ignored the parameter τ i , the time at which the radio source just becomes isotropic. This way, we assign it the same value, 3 yr, considered by Levinson et al. (2002). We note that the dependence on this parameter is weak

emission line. While detailed nebular-phase photoionization calculations are currently challenging, a preliminary examination of spectra produced using CLOUDY (Ferland et al. 2013; see methods in Margalit et al. 2018) indicates that, for parameters typical of SLSN magnetars and their ejecta, strong [O III] emission broadly consistent with J1419+3940 is common at ∼20 yr post-explosion. However, more detailed calculations of the nebular phase of pulsar/magnetar-powered SNe would be required to confirm this possibility and its implications for the ejecta structure.

Conclusions

Based on all the observational data and our analysis of FIRST J141918.9+394036, we arrive at the following conclusions.

1. J1419+3940 is an unprecedented (radio-discovered, luminous, decades-long) transient having a peak radio luminosity > 2.3 × 10 29 erg s -1 Hz -1 at 1.4 GHz and detectable radio emission > 26 yr post-explosion. 2. Average blastwave velocity is >44,000 km s -1 in 1993

(assuming the first radio detection epoch is <1000 days post-explosion). If the rise of the light curve is as rapid as the decline then the initial velocity of the transient is 0.2 c. 3. Average blastwave velocity is ;19,000 km s -1 in 2019 (last observing epoch; assuming ∼26 yr post-explosion). 4. Age-and CSM-density-independent estimate of the blastwave energy is ∼ 5 × 10 50 erg (dependent on the microphysical parameters). 5. Optical spectroscopic observations from 2019 reveal a broad [O III]λ4959,5007 emission line. We find that collisional excitation in compact star-forming region(s) within the host galaxy is the most conservative explanation, but we cannot completely rule out its association with the transient. A transient origin for the broad line could suggest the presence of a magnetar. 6. Host galaxy properties are suggestive of a massive star progenitor that is more likely to form in high-specific-SFR or low-metallicity environments, similar to those observed for LGRBs, SN Ic-BL, SLSN-I, and AT2018cow-like FBOTs. We are able to rule out SLSN-I (and neutron star merger ejecta scenario proposed by Lee et al. 2020) based on rates and peak radio luminosity, and find that SN Ic-BL (not associated with GRBs) is very unlikely based on the energetics. 7. The observed afterglow properties of J1419+3940 are most consistent with those of LGRBs in terms of previously studied transients. The afterglow light curve is especially similar to the late-time evolution of GRB 030329 if we assume that the first radio detection occurred ∼1000 days post-explosion.  degrees (68% confidence). 9. The late-time radio light curve of J1419+3940 may reveal the presence of a GRB-SN and continued radio monitoring of J1419+3940 is therefore recommended. 10. The rates of J1419+3940-like events, which we find to be 4 -24 × 10 -5 deg -2 , or equivalently about 40-240 Gpc -3 yr -1 , suggest that the VLA Sky Survey and surveys with the ASKAP and MeerKAT will find a few such events over the coming years.

Taking the ambient density to be either a constant-density ISM, or an r -2 wind, Combining the above equations, we find that the optically thin synchrotron luminosity is For the case where p = 3, we find quantitatively that where above we have expressed the shock velocity as v = mR/t (Equation ( 2)). This is the same as Equation (1) in the main text, and that is used to infer source properties.