A Mildly Relativistic Outflow Launched Two Years after Disruption in Tidal Disruption Event AT2018hyz

AT2018hyz was observed in targeted observations with the Arcminute Microkelvin Imager Large Array (AMI-LA) 32 days after optical discovery leading to an upper limit of F_ν (15.5 GHz) ≲ 85 μJy (Horesh et al. 2018), and with the Atacama Large Millimeter/submillimeter Array (ALMA) at 45 and 66 days to limits of F_ν (100 GHz) ≲ 38 and 48 μJy, respectively (Gomez et al. 2020). The location of AT2018hyz was also covered in the Australian Square Kilometre Array Pathfinder (ASKAP) Variables and Slow Transients Survey (VAST; Murphy et al. 2021) at 697 days with a limit of F_ν (0.9 GHz) ≲ 0.9 mJy, and in the Very Large Array (VLA) Sky Survey epoch 2.1 (VLASS; Lacy et al. 2020) at 705 days with a limit of F_ν (3 GHz) ≲ 0.45 mJy. We determine the flux density limits in the survey images with the imtool fitsrc command within the pwkit package ¹⁰ (Williams et al. 2017).

As part of a broader study of late-time radio emission from TDEs, we observed AT2018hyz with the Karl G. Jansky Very Large Array (VLA) at 972 days post optical discovery in the C band (Program ID 21A-303, PI: Hajela) and detected a source with a F_ν (5 GHz) = 1.39 ± 0.02 mJy and F_ν (7 GHz) = 1.00 ±0.02 mJy (see Table 1). Following this initial detection we obtained multifrequency observations spanning from the L to the K band (≈1–23 GHz; Programs 21B-357, 21B-360, and 22A-458, PI: Cendes), which resulted in detections across the full frequency range. For all observations we used the primary calibrator 3C147 and the secondary calibrator J1024-0052. We processed the VLA data using standard data reduction procedures in the Common Astronomy Software Application package (CASA; McMullin et al. 2007), using tclean on the calibrated measurement set available in the NRAO archive, using Briggs weighting. Our images were 2048 × 2048 pixels, with a cell size of $\tfrac{1}{3}$ the size of the synthesized beamwidth for a given band and array configuration ¹¹ ; no other bright sources within the region that required additional cleaning. The observations and resulting flux density measurements are summarized in Table 1. Additionally, we use data collected by the commence of the VLA Low-band Ionosphere and Transient Experiment (VLITE; Clarke et al. 2016) at 350 MHz during our multifrequency observations (Table 1). We note that the uncertainties listed in Table 1 are statistical only and do not include a ≈3%–5% systematic uncertainty in the overall flux density calibration; we account for this systematic uncertainty in our subsequent modeling (see Section 4).

We also obtained observations with the MeerKAT radio telescope in the UHF and L band (0.8–2 GHz) on 2022 April 17 (1282 days; DDT-20220414-YC-01, PI: Cendes) and with the Australian Telescope Compact Array (ATCA) on 2022 May 1 (1296 days; Program C3472; PI: Cendes) at 2–20 GHz. For ATCA we reduced the data using the MIRIAD package. Calibrator 1934-638 was used to calibrate absolute flux density and bandpass, while calibrator 1038+064 was used to correct short term gain and phase changes. The invert, mfclean, and restor tasks were used to make deconvolved wideband, natural weighted images in each frequency band. For MeerKAT, we used the flux calibrator 0408-6545 and the gain calibrator 3C237, and used the calibrated images obtained via the South African Radio Astronomy Observatory Science Data Processor (SDP). ¹² We confirmed via the secondary SDP products that the source fluxes in the MeerKAT images were ∼90% of the sources overlapping with the NRAO VLA Sky Survey (NVSS; Condon et al. 1998), and that frequency slices show a steady increase in flux within the MeerKAT frequency range (with a spectral index α ≈−1).

Following the initial VLA radio detection we also observed AT2018hyz with the ALMA at 1141, 1201, and 1253 days (Project 2021.1.01210.T, PI: Alexander). These observations roughly coincide with the three multifrequency VLA observations. The first and third ALMA observations were in band 3 (mean frequency of 97.5 GHz) and band 6 (240 GHz); the second epoch was in band 3 only. For the ALMA observations, we used the standard NRAO pipeline (version: 2021.2.0.128) in CASA (version: 6.2.1.7) to calibrate and image the data. We detect AT2018hyz in all observations with a rising flux density (Table 1). We note that in the first observation, the flux density is at least 30 times higher than the early upper limits.

We obtained a Director's Discretionary Time observation of AT2018hyz on 2022 March 19 (δ t = 1253 days) with ACIS-S on board the Chandra X-ray Observatory with an exposure time of 14.9 ks (Program 23708833, PI: Cendes). We processed the data with chandra_repro within CIAO 4.14 using the latest calibration files. An X-ray source is detected with wavdetect at the position of AT2018hyz with a net count rate of (8.9 ± 2.5) × 10⁻⁴ c s⁻¹ (0.5 − 8 keV), with statistical confidence of 6σ (Gaussian equivalent).

We extracted a spectrum of the source with specextract using a 15 radius source region and a source-free background region of 22'' radius. We fit the spectrum with an absorbed simple power-law model (tbabs ∗ztbabs ∗pow within Xspec). The Galactic neutral hydrogen column density in the direction of AT2018hyz is N_H,MW = 2.67 × 10²⁰ cm⁻² (Kalberla et al. 2005), and we find no evidence of additional intrinsic absorption, to a 3σ upper limit of N_H,int ≲ 2.8 × 10²² cm⁻². The photon index is Γ = 1.5 ± 0.7 (1σ confidence level), and the 0.3 − 10 keV unabsorbed flux is ${F}_{x}={1.78}_{-0.28}^{+0.35}\times {10}^{-14}\,\mathrm{erg}\,{{\rm{s}}}^{-1}{\mathrm{cm}}^{-2}$. The flux density at 1 keV is 4.03 ± 1.20 × 10⁻⁶ mJy.

The X-ray flux is a factor of about 2 lower than the flux measured with XRT in the first 86 days, ${F}_{x}={4.1}_{-0.4}^{+0.6}\times {10}^{-14}\,\mathrm{erg}\,{{\rm{s}}}^{-1}{\mathrm{cm}}^{-2}$ (Gomez et al. 2020). This indicates that AT2018hyz has faded in X-rays over this timescale, and we discuss the implications of this in Section 4.3.

We obtained UV observations of AT 2018hyz from the UV/Optical Telescope (UVOT; Roming et al. 2005) on board the Neil Gehrels Swift observatory (Swift; Gehrels et al. 2004) on 2022 January 7 (1182 days post discovery) to continue tracking the evolution of the UV light curve originally presented in Gomez et al. (2020). We measure an AB magnitude of m_UVW2 = 20.20 ± 0.08 mag using a 5''aperture centered on the position of AT 2018hyz using the HEAsoft uvotsource function (HEARSARC 2014). We correct this magnitude for Galactic extinction using the Barbary (2016) implementation of the Cardelli et al. (1989) extinction law (0.28 mag) and subtract the host contribution of 21.76 mag (Gomez et al. 2020) to obtain a final value of m_UVW2 = 20.14 ± 0.08 mag.

We compare this measurement to the MOSFiT model of AT 2018hyz presented in Gomez et al. (2020), which fit all optical/UV data to about 300 days post discovery. The model predicts a UVW2 magnitude of 23.3 ± 0.2 at the time of our new observation, a factor of 18 times dimmer than observed. This points to excess emission in the UV band.

Similarly, we measure the late-time r-band magnitude of AT 2018hyz at ∼1030 days post discovery by downloading raw Zwicky Transient Facility (ZTF) images from the NASA/IPAC Infrared Science Archive ¹³ and combining all r-band images taken within ±12 days of MJD = 59501. We subtract a template archival pre-explosion ZTF image from the combined image using HOTPANTS (Becker 2015), and perform point-spread function (PSF) photometry on the residual image. For calibration, we estimate the zero-point by measuring the magnitudes of field stars and comparing to photometric AB magnitudes from the PS1/3π catalog. Corrected for galactic extinction, we find r = 21.83 ± 0.36 mag. At this phase, the MOSFiT model predicts a magnitude of r = 24.65 ± 0.30, a factor of 13 times dimmer than measured. We do not detect emission in the g band to a 3σ limit of 21.76 mag. This implies a late-time color of g − r ≳ − 0.1, redder than the latest measurements from Gomez et al. (2020) of g − r =− 0.58 ± 0.23.

The radio/millimeter SEDs, shown in Figure 3, exhibit a power-law shape with a turnover and peak at ≈1.5 GHz through 1251 days. At 1282 days, however, the peak of the SED shifts upwards to ≈3 GHz. A rapid shift to a higher peak frequency is unprecedented in radio observations of TDEs. The power-law shape above the peak is characteristic of synchrotron emission.

Zoom In Zoom Out Reset image size

We fit the SEDs with the model of Granot & Sari (2002), developed for synchrotron emission from gamma-ray burst (GRB) afterglows, and previously applied to the radio emission from other TDEs (e.g., Zauderer et al. 2011; Cendes et al. 2021b), using specifically the regime ¹⁴ of ν_m ≪ ν_a :

$$\begin{eqnarray}\begin{array}{rcl}{F}_{\nu } & = & {F}_{\nu }({\nu }_{m})\left[{\left(\displaystyle \frac{\nu }{{\nu }_{m}}\right)}^{2}{e}^{-{s}_{4}{\left(\nu /{\nu }_{m}\right)}^{2/3}}+{\left(\displaystyle \frac{\nu }{{\nu }_{m}}\right)}^{5/2}\right]\\ & & \times {\left[1+{\left(\displaystyle \frac{\nu }{{\nu }_{a}}\right)}^{{s}_{5}({\beta }_{2}-{\beta }_{3})}\right]}^{-1/{s}_{5}},\end{array}\end{eqnarray} \tag{ 1 }$$

where β₂ = 5/2, β₃ = (1 − p)/2, s₄ = 3.63p − 1.60, and s₅ = 1.25 − 0.18p (Granot & Sari 2002). Here, p is the electron energy distribution power-law index, $N({\gamma }_{e})\propto {\gamma }_{e}^{-p}$ for γ_e ≥ γ_m , ν_m is the frequency corresponding to γ_m , ν_a is the synchrotron self-absorption frequency, and F_ν (ν_m ) is the flux normalization at ν = ν_m .

We determine the best-fit parameters of the model—F_ν (ν_m ), ν_a , and p—using the Python Markov Chain Monte Carlo (MCMC) module emcee (Foreman-Mackey et al. 2013), assuming a Gaussian likelihood where the data have a Gaussian distribution for the parameters F_ν (ν_m ) and ν_a . For p we use a uniform prior of p = 2 − 3.5. We also include in the model a parameter that accounts for additional systematic uncertainty beyond the statistical uncertainty on the individual data points, σ, which is this a fractional error added to each data point. The posterior distributions are sampled using 100 MCMC chains, which were run for 3000 steps, discarding the first 2000 steps to ensure the samples have sufficiently converged by examining the sampler distribution. The resulting SED fits are shown in Figure 3 and provide a good fit to the data.

Our first observation at 972 days only includes 5 and 7 GHz, but clearly points to an optically thin spectrum with a peak at ≲5 GHz. We combine this observation with a VAST detection of the source at 1013 days with a flux of 1.3 ± 0.03 mJy at 0.89 GHz (Horesh et al. 2022), indicating the peak is between these two frequencies. We make the reasonable assumption that the lack of evolution in ν_p between 1126–1251 days indicates no serious change at 972 days as well and fix p = 2.30. We use these values to determine the physical properties of the outflow at 972 days.

From the SED fits we determine the peak frequency and flux density, ν_p and F_ν,p , respectively, which are used as input parameters for the determination of the outflow physical properties. The best-fit values and associated uncertainties are listed in Table 2. We find that ν_p remains essentially constant at ≈1.3–1.6 GHz at 972 to 1251 days, while F_ν,p increases steadily by a factor of 3.7. While the VLITE limits at 350 MHz lie above our SED model fits (Figure 8) we note that they require an SED peak at ≳0.5 GHz as otherwise these limits would be violated. In addition, the single power-law shape above ν_p indicates that the synchrotron cooling frequency is ν_c ≳ 240 GHz at 1126 and 1251 days. For the SED at 1282 days we find that F_ν,p has remained steady, while ν_p increased by a factor of 2. We also note that the spectral index below the peak appears to be shallower than F_ν ∝ ν^5/2.

Given the unusual evolution to higher ν_p in the latest epoch, we have also considered a model with two emission components, peaking at ν_p,1 ≈ 1.5 and ν_p,2 ≈ 3 GHz, in which the lower-frequency component dominates the emission at an early time, and the higher-frequency component rises at later times. The results of this model are presented in the Appendix, but in the main paper we focus on the simpler single-component model.

Using the inferred values of ν_p , F_ν,p , and p from Section 4.1, we can now derive the physical properties of the outflow and ambient medium using an equipartition analysis. In all epochs, we assume a mean value of p = 2.3 in our calculations. We first assume the conservative case of a nonrelativistic spherical outflow using the following expressions for the radius and kinetic energy (Equations (27) and (28) in Barniol Duran et al. 2013):

$$\begin{eqnarray}\begin{array}{rcl}{R}_{\mathrm{eq}} & \approx & (1\times {10}^{17}\,\mathrm{cm})\times {\left(21.8\times {525}^{p-1}\right)}^{\tfrac{1}{13+2p}}{\chi }_{e}^{\tfrac{2-p}{13+2p}}\\ & & \times \,{F}_{p,\mathrm{mJy}}^{\tfrac{6+p}{13+2p}}{d}_{L,28}^{\tfrac{2(p+6)}{13+2p}}{\nu }_{p,10}^{-1}{\left(1+z\right)}^{\tfrac{-(19+3p)}{13+2p}}{f}_{A}^{\tfrac{-(5+p)}{13+2p}}{f}_{V}^{\tfrac{-1}{13+2p}}\\ & & \times {4}^{\tfrac{1}{13+2p}}{\epsilon }^{1/17}{\xi }^{\tfrac{1}{13+2p}}{{\rm{\Gamma }}}^{\tfrac{p+8}{13+2p}}{\left({\rm{\Gamma }}-1\right)}^{\tfrac{2-p}{13+2p}},\end{array}\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}\begin{array}{rcl}{E}_{\mathrm{eq}} & \approx & (1.3\times {10}^{48}\,\mathrm{erg})\times [{21.8}^{\tfrac{-2(p+1)}{13+2p}}]{\left[{525}^{p-1}\times {\chi }_{e}^{2-p}\right]}^{\tfrac{11}{13+2p}}\\ & & \times {F}_{p,\mathrm{mJy}}^{\tfrac{14+3p}{13+2p}}{d}_{L,28}^{\tfrac{2(3p+14)}{13+2p}}{\nu }_{p,10}^{-1}{\left(1+z\right)}^{\tfrac{-(27+5p)}{13+2p}}\\ & & \times \,{f}_{A}^{\tfrac{-(3(p+1)}{13+2p})}{f}_{V}^{\tfrac{2(p+1)}{13+2p}}{4}^{\tfrac{11}{13+2p}}{\xi }^{\tfrac{11}{13+2p}}[(11/17){\epsilon }^{-6/17}\\ & & +(6/17){\epsilon }^{11/17}]\,{{\rm{\Gamma }}}^{\tfrac{-5p+16}{13+2p}}{\left({\rm{\Gamma }}-1\right)}^{\tfrac{-11(p-2)}{13+2p}},\end{array}\end{eqnarray} \tag{ 3 }$$

where d_L = 204 Mpc is the luminosity distance and z = 0.0457 is the redshift. The factors f_A and f_V are the area and volume filling factors, respectively, where in the case of a spherical outflow f_A = 1 and ${f}_{V}=\tfrac{4}{3}\times (1-{0.9}^{3})\approx 0.36$ (i.e., we assume that the emitting region is a shell of thickness 0.1R_eq), while in the case of a jet ¹⁵ ${f}_{A}={\left({\theta }_{j}{\rm{\Gamma }}\right)}^{2}$ and ${f}_{V}=0.27{\left({\theta }_{j}{\rm{\Gamma }}\right)}^{2}$. The factors of 4^1/13+2p and 4^11/13+2p in R_eq and E_eq, respectively, arise from corrections to the isotropic number of radiating electrons (N_e,iso) in the nonrelativistic case. We further assume that the fraction of post-shock energy in relativistic electrons is _e = 0.1, which leads to correction factors of ξ^1/13+2p and ξ^11/13+2p in R_eq and E_eq, respectively, with $\xi =1+{\epsilon }_{e}^{-1}\approx 11$. We parameterize any deviation from equipartition with a correction factor = (11/6)(_B /_e ), where _B is the fraction of post-shock energy in magnetic fields (here we use _B = 0.1; see Section 4.3). Finally, χ_e = (p − 2)/(p − 1)_e (m_p /m_e ) ≈ 42.3, where m_p and m_e are the proton and electron masses, respectively.

Using R_eq we can also determine additional parameters of the outflow and environment (Barniol Duran et al. 2013): the magnetic field strength (B), the number of radiating electrons (N_e ), and the Lorentz factor of electrons radiating at ν_a (γ_a ):

$$\begin{eqnarray}\begin{array}{rcl}B & \approx & (1.3\times {10}^{-2}\,{\rm{G}})\times {F}_{p,\mathrm{mJy}}^{-2}{d}_{L,28}^{-4}{\left(1+z\right)}^{7}\\ & & \times \,{\nu }_{p,10}^{5}\displaystyle \frac{{f}_{A}^{2}{R}_{\mathrm{eq},17}^{4}}{{{\rm{\Gamma }}}^{3}},\end{array}\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}\begin{array}{rcl}{N}_{e} & \approx & (4\times {10}^{54})\times {F}_{p,\mathrm{mJy}}^{3}{d}_{L,28}^{6}{\left(1+z\right)}^{-8}{f}_{A}^{-2}{R}_{\mathrm{eq},17}^{-4}{\nu }_{p,10}^{-5}\\ & & \times \,{\left({\gamma }_{a}/{\gamma }_{m}\right)}^{(p-1)},\end{array}\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&{\gamma }_{a}\approx 525\times {F}_{p,\mathrm{mJy}}{d}_{L,28}^{2}{\nu }_{p,10}^{-2}{\left(1+z\right)}^{-3}\displaystyle \frac{{\rm{\Gamma }}}{{f}_{A}{R}_{\mathrm{eq},17}^{2}}.\end{eqnarray} \tag{ 6 }$$

We note an additional factor of 4 and a correction factor of ${\left({\gamma }_{a}/{\gamma }_{m}\right)}^{p-1}$ are added to N_e for the nonrelativistic regime (Barniol Duran, private communication); here we use ${\gamma }_{m}=\max [{\chi }_{e}({\rm{\Gamma }}-1),2]$. We determine the ambient density assuming a strong shock and an ideal monoatomic gas as n_ext = N_e /4V, where the factor of 4 is due to the shock jump conditions, and V is the volume of the emitting region, $V={f}_{V}\pi {R}_{\mathrm{eq}}^{3}/{{\rm{\Gamma }}}^{4}$, with f_V as defined above.

Using the inferred physical parameters we can also predict the location of the synchrotron cooling frequency, given by Sari et al. (1998):

$$\begin{eqnarray}&&{\nu }_{c}\approx 2.25\times {10}^{14}{B}^{-3}{{\rm{\Gamma }}}^{-1}{t}_{d}^{-2}\,\mathrm{Hz}.\end{eqnarray} \tag{ 7 }$$

As ν_c has a strong dependence on the magnetic field strength, measuring its location directly can determine _B and whether the outflow deviates from equipartition.

In Figure 4 we show our VLA+ALMA+Chandra SED at 1251–1253 days, along with our model SED from Section 4.1, which does not include a cooling break (dashed lines). This model clearly overpredicts the Chandra measurements. The steepening required by the Chandra data is indicative of a cooling break, which we model with an additional multiplicative term to Equation (1) of ${\left[1+{\left(\nu /{\nu }_{c}\right)}^{{s}_{3}({\beta }_{3}-{\beta }_{4})}\right]}^{-1/{s}_{3}}$, where β₄ = − p/2 and we use s₃ = 10 (Granot & Sari 2002; Cendes et al. 2021a). Fitting this model to the data we find ν_c ≈ 10¹³ Hz (solid lines in Figure 4). However, as the X-ray flux measured with Chandra is only a factor of 2 times fainter than the steady early time X-ray flux, and hence can be due to a source other than the radio-emitting outflow, we consider it to be an upper limit on the contribution of the radio-emitting outflow. As a result, our estimate of ν_c is actually an upper limit.

Zoom In Zoom Out Reset image size

With the value of ν_c determined, we adjust the value of _B and solve Equation (7) after repeating the equipartition analysis (Equations (2)–(5)) to account for the deviation from equipartition in those parameters. With this approach, we find that _B ≈ 0.01. Given that the deviation is not significant, and that ν_c is an upper limit (and hence _B ≈ 0.01 is a lower limit) we conservatively assume equipartition (_B = 0.1) in our subsequent analysis and in Table 3. We emphasize that the change would be relatively minor if we adjusted _B - at 1251 days, for the jetted model we find this would correspond with a radius decrease from $\mathrm{log}({R}_{\mathrm{eq}})\approx 18.21$ to ≈18.16, and the energy would increase from $\mathrm{log}({E}_{\mathrm{eq}})\approx 49.80$ to ≈49.98.

Table 3. Equipartition Model Parameters

Geometry	δ t	log(Req)	log(Eeq)	log(B)	log(Ne )	log(next)	Γ	β	γa	νc
	(day)	(cm)	(erg)	(G)		(cm−3)				(GHz)
spherical	972	${17.22}_{-0.06}^{+0.09}$	${49.03}_{-0.09}^{+0.13}$	$-{0.79}_{-0.07}^{+0.05}$	${52.54}_{-0.09}^{+0.13}$	${0.27}_{-0.14}^{+0.10}$	${1.03}_{-0.01}^{+0.01}$	${0.23}_{-0.03}^{+0.03}$	${70}_{-1}^{+1}$	${1061}_{-1212}^{+1039}$
(fA = 1,	1126	${17.41}_{-0.01}^{+0.01}$	${49.47}_{-0.01}^{+0.01}$	$-{0.85}_{-0.01}^{+0.01}$	${52.96}_{-0.01}^{+0.01}$	${0.12}_{-0.03}^{+0.03}$	${1.02}_{-0.01}^{+0.01}$	${0.22}_{-0.01}^{+0.01}$	${63}_{-1}^{+1}$	${579}_{-62}^{+48}$
fV = 0.36,	1199	${17.57}_{-0.03}^{+0.03}$	${49.76}_{-0.03}^{+0.03}$	$-{0.96}_{-0.03}^{+0.03}$	${53.28}_{-0.03}^{+0.03}$	$-{0.04}_{-0.06}^{+0.06}$	${1.03}_{-0.01}^{+0.01}$	${0.25}_{-0.01}^{+0.01}$	${65}_{-1}^{+1}$	${838}_{-209}^{+158}$
B = 0.1)	1251	${17.54}_{-0.02}^{+0.02}$	${49.75}_{-0.02}^{+0.02}$	$-{0.91}_{-0.02}^{+0.02}$	${53.30}_{-0.02}^{+0.02}$	$-{0.02}_{-0.03}^{+0.04}$	${1.02}_{-0.01}^{+0.01}$	${0.22}_{-0.01}^{+0.01}$	${65}_{-1}^{+1}$	${492}_{-77}^{+62}$
	1282	${17.24}_{-0.04}^{+0.04}$	${49.57}_{-0.05}^{+0.04}$	$-{0.56}_{-0.05}^{+0.05}$	${52.92}_{-0.05}^{+0.04}$	${0.54}_{-0.10}^{+0.09}$	${1.01}_{-0.01}^{+0.01}$	${0.12}_{-0.01}^{+0.01}$	${62}_{-1}^{+1}$	${377}_{-122}^{+165}$
10◦ jet	972	${17.89}_{-0.06}^{+0.09}$	${49.09}_{-0.09}^{+0.13}$	$-{1.03}_{-0.07}^{+0.05}$	${52.65}_{-0.09}^{+0.13}$	${0.02}_{-0.14}^{+0.11}$	${1.20}_{-0.04}^{+0.04}$	${0.55}_{-0.04}^{+0.05}$	${96}_{-1}^{+1}$	${5465}_{-1196}^{+1034}$
(fA = θ2Γ2,	1126	${18.08}_{-0.01}^{+0.01}$	${49.52}_{-0.01}^{+0.01}$	$-{1.09}_{-0.01}^{+0.01}$	${53.09}_{-0.01}^{+0.01}$	$-{0.13}_{-0.03}^{+0.03}$	${1.21}_{-0.01}^{+0.01}$	${0.56}_{-0.01}^{+0.01}$	${76}_{-1}^{+1}$	${2965}_{-318}^{+247}$
fV = 0.27fA )	1199	${18.24}_{-0.03}^{+0.03}$	${49.81}_{-0.03}^{+0.03}$	$-{1.20}_{-0.03}^{+0.03}$	${53.36}_{-0.03}^{+0.03}$	$-{0.29}_{-0.06}^{+0.06}$	${1.26}_{-0.02}^{+0.02}$	${0.61}_{-0.02}^{+0.02}$	${78}_{-1}^{+1}$	${4341}_{-1082}^{+819}$
B = 0.1)	1251	${18.21}_{-0.02}^{+0.02}$	${49.80}_{-0.02}^{+0.02}$	$-{1.15}_{-0.02}^{+0.02}$	${53.37}_{-0.02}^{+0.02}$	$-{0.24}_{-0.04}^{+0.04}$	${1.21}_{-0.01}^{+0.01}$	${0.57}_{-0.01}^{+0.01}$	${78}_{-1}^{+1}$	${2523}_{-388}^{+322}$
	1282	${17.93}_{-0.04}^{+0.04}$	${49.60}_{-0.05}^{+0.04}$	$-{0.79}_{-0.05}^{+0.05}$	${53.18}_{-0.05}^{+0.04}$	${0.31}_{-0.10}^{+0.09}$	${1.09}_{-0.01}^{+0.01}$	${0.39}_{-0.01}^{+0.02}$	${79}_{-1}^{+1}$	${188}_{-72}^{+56}$

Note. Values in this table are calculated using an outflow launch time of t₀ = 750 days. For Γ and β we have accounted for the uncertainty in the launch date (${t}_{0,\mathrm{sphere}}={750}_{-127}^{+80}$ days and ${t}_{0,\mathrm{jet}}={750}_{-113}^{+73}$).

Download table as:  ASCIITypeset image

We begin by investigating the properties of the outflow in the conservative case of spherical geometry. We summarize the inferred physical parameters for all epochs in Table 3. We find that the radius increases from $\mathrm{log}({R}_{\mathrm{eq}})\approx 17.22$ to ≈17.57 between 970 and 1251 days, corresponding to a large velocity of β ≈ 0.28 over this time span. However, if we use the time of optical discovery as the outflow launch date, then the inferred velocity in the first epoch (972 days) is β ≈ 0.066 and in the fourth epoch (1251 days) it is β ≈ 0.11. This means that the assumption of a launch date that coincides with the optical discovery is incorrect. Instead, we find that the increase in radius during the first four epochs is roughly linear, and fitting such an evolution with the launch time (i.e., time at which R = 0) as a free parameter, we instead find t₀ ≈ 750 d; see Figure 5. Thus, the physical evolution of the radius during this period confirms our initial argument based on the rapid rise of the radio emission (Section 3) that the outflow was launched with a substantial delay of about 2 yr relative to the optical emission. The inferred outflow velocity using the four epochs is β ≈ 0.24, which is larger than the typical velocities of β ≈ 0.1 inferred for previous nonrelativistic radio-emitting TDEs (Alexander et al. 2016, 2017; Anderson et al. 2019; Cendes et al. 2021a; Goodwin et al. 2022).

Zoom In Zoom Out Reset image size

The outflow kinetic energy increases by a factor of about 5 from E_K ≈ 1.1 × 10⁴⁹ to 5.8 × 10⁴⁹ erg between 970 and 1250 days. The rapid increase in energy indicates that we are observing the deceleration of the outflow. The kinetic energy is larger than that of previous radio-emitting TDEs by a factor of ≈4 (Anderson et al. 2019; Goodwin et al. 2022), although we note that this energy is a lower limit and would be higher if a deviation from equipartition is assumed.

Finally, we find that the inferred ambient density declines from n_ext ≈ 1.9 to ≈1.0 cm⁻³ over the distance scale of 1.7 × 10¹⁷ to 3.7 × 10¹⁷ cm, or ρ(R) ∝ R⁻¹. This is similar to the rate seen in M87* and Sw J1644+57 and Sgr A*, and less steep than the density profiles inferred around previous thermal TDEs (see Figure 7 and citations therein). Combined with the mild decline in density with radius, this indicates that the late turn-on of the radio emission is inconsistent with a density jump.

To conclude, even in the conservative spherical scenario we find that radio emission requires a delayed, mildly relativistic outflow with a higher velocity and energy than in previous radio-emitting TDEs.

In light of the large outflow velocity and kinetic energy in the spherical case, we also consider the results for a collimated outflow. In particular, we choose an outflow opening angle of 10°, typical of GRB jets and the jet in Sw J1644+57 (Berger et al. 2012; Zauderer et al. 2013). For a collimated outflow the resulting radius (and hence velocity) are larger than in the spherical case, so we need to also solve for the Lorentz factor, Γ, which impacts the values of R_eq (as well as the other physical parameters). We begin by using the launch time inferred in the spherical case, and then iteratively recalculate Γ, R_eq, and the launch date. In the process we also include the relevant modifications due to Γ in f_A and f_V , which also impact the value of R_eq (Equation (2)). We note that these corrections are relatively small but we account for them for completeness.

With this approach, we find that the launch date for a 10° jet is t₀ ≈ 750 days post optical discovery. The resulting radii are a factor of 4.7 times larger than in the spherical case, leading to a mean velocity to 1251 days of β ≈ 0.57, or Γ ≈ 1.23. The kinetic energy is E_K ≈ 6.3 × 10⁴⁹ erg at 1251 days. Finally, we find the density declines from n_ext ≈ 1.0 to ≈0.5 cm⁻³ over the distance scale of 7.8 × 10¹⁷ to 1.7 × 10¹⁸ cm. This is less dense than the spherical case but follows a similar profile of ρ(R) ∝ R⁻¹ (Section 5.1).

We can also consider the possibility that the late-time emission from AT2018hyz is caused by a relativistic jet with an initial off-axis viewing orientation. The radio emission from an off-axis relativistic jet will be suppressed at early times by relativistic beaming but will eventually rise rapidly (as steep as t³) when the jet decelerates and spreads. The time and radius at which the radio emission will peak are given by (e.g., Nakar & Piran 2011):

$$\begin{eqnarray}&&{t}_{\mathrm{dec}}\approx 30\,{\rm{d}}\,{E}_{\mathrm{eq},49}^{1/3}{n}_{\mathrm{ext}}^{-1/3}{\beta }_{0}^{-5/3},\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&{R}_{\mathrm{dec}}\approx {10}^{17}\,\mathrm{cm}\,{E}_{\mathrm{eq},49}^{1/3}{n}_{\mathrm{ext}}^{-1/3}{\beta }_{0}^{-2/3},\end{eqnarray} \tag{ 9 }$$

where β₀ is the initial velocity, which for an off-axis relativistic jet is β₀ = 1. Using the equipartition parameters in the spherical case, we find t_dec ≈ 25 days and R_dec ≈ 8.3 × 10¹⁶ cm. The value of t_dec is substantially smaller than the observed delay of ∼10³ days (as is R_dec). This agrees with our general argument in Section 3 that the lower radio luminosity of AT2018hyz and its much later appearance compared to the radio emission of Sw J1644+57 (which became nonrelativistic at ≈700 days) argues against an off-axis jet launched at the time of disruption. Moreover, the radio emission from an off-axis jet is expected to rise no steeper than t³ (Nakar & Piran 2011), whereas in AT2018hyz the emission rises as t⁵ if the outflow was launched at the time of disruption.

We can consider the possibility of an off-axis jet expanding initially into a low-density cavity, followed by a denser region (thus delaying t_dec to a longer timescale as observed). However, we can rule out this model because the observed rise in radio emission spans ≈300 days, while Equation (8) indicates deceleration over a timescale of t_dec ≈ 30 days even if the time at which deceleration starts is itself delayed. Similarly, if there was initially a higher density environment which we did not capture in our observations, this would only correspond to a faster t_dec.

We also consider the hypothesis that the rise in radio emission is caused by unbound material from the initial disruption colliding with a surrounding dense circumstellar material (CSM). Theoretical modeling of such unbound debris indicates the fastest speeds reached are v ≈ 0.03c (Guillochon et al. 2016; Yalinewich et al. 2019), which is significantly smaller than what we infer in both the spherical and jetted cases for AT2018hyz. We thus conclude this scenario is unlikely.

Finally, we considered a model in which the change in the SED properties during the latest epoch could be due to a combination of two outflows, with one dominating at a lower frequency and fading and the second dominating at higher frequencies and rising (see Appendix). Such a model may be expected if internal shocks within the outflow are leading to dissipation at more than one radius. However, we find that such a model requires the initial emission component to decline very rapidly (effectively turn off), while the later emission component has to rise by more than an order of magnitude in only 30 days between about 1250 and 1280 days. Such rapid evolution does not seem feasible, even in the context of AT2018hyz.

In Figure 6 we plot the kinetic energy and velocity (Γβ) of the delayed outflow in AT2018hyz in comparison to previous TDEs for which a similar analysis has been carried out, using the highest energy inferred in those sources (Zauderer et al. 2011; Alexander et al. 2016, 2017; Anderson et al. 2019; Cendes et al. 2021a; Cendes et al. 2021b; Stein et al. 2021; Goodwin et al. 2022). We find that at its peak, in the spherical case the energy is ≈4 times larger and the velocity is ≈2 times faster than in previous nonrelativistic TDEs. If we compare to ASASSN-15oi specifically, using the observation with the highest peak frequency and peak flux (i.e., 182 days post optical discovery), with _B = 0.1 and p = 2.4 (which best fits the SED; see Horesh et al. 2021a), we obtain an energy for ASASSN-15oi that is ≈23 times lower than in AT2018hyz. To infer the velocity in ASASSN-15oi, we subtract 90 days from the date of the observation (the last date of nondetection; see Horesh et al. 2021a) to find β ≈ 0.13, which is lower than β ≈ 0.25 for AT2018hyz.

Zoom In Zoom Out Reset image size

If the outflow in AT2018hyz is collimated, then its velocity is significantly higher than in previous nonrelativistic TDEs, placing it in an intermediate regime with the powerful jetted TDEs such as Sw J1644+57.

In Figure 7 we plot the ambient density as a function of the radius (scaled by the Schwarzschild radius) for AT2018hyz and previous radio-emitting TDEs. Here we use M_BH ≈ 5.2 ×10⁶ M_⊙ for AT2018hyz, as inferred by Gomez et al. (2020). We find that the density decreases with the radius and is consistent with the densities and circumnuclear density profiles of previous TDEs, including Sw J1644+57. Crucially, we do not infer an unusually high density, which might be expected if the radio emission was delayed due to rapid shift from low to high density.

Zoom In Zoom Out Reset image size

Two previous TDEs with delayed radio emission have been published recently. The radio emission from iPTF16fnl is detected on a much earlier timescale than in AT2018hyz and rises more gradually to a peak luminosity of only ≈10³⁷ erg s⁻¹ (Horesh et al. 2021b). This is almost a factor of 100 less luminous than AT2018hyz. We do not consider the radio emission from iPTF16fnl to be similar in nature to that from AT2018hyz.

On the other hand, ASASSN-15oi exhibits two episodes of rapid brightening, at ≈200 and ≈1400 days (Horesh et al. 2021a). While the second brightening is not well characterized temporally or spectrally, it has a comparable rise rate and luminosity to AT2018hyz. We speculate that it may be due to a delayed outflow with similar properties to that of AT2018hyz, including a delay of several hundred days, which would make it distinct from the first peak in the ASASSN-15oi light curve.

There are at least two broad possibilities for the origin of the delayed mildly relativistic outflow, both of which connect to its assumed origin in a fast disk wind or jet (hereafter, collectively referred to as jet) from the innermost regions of the black hole accretion flow.

One possibility is that the jet was weak or inactive at early times after the disruption and then suddenly became activated at ≈750 days. Such sudden activation could result from a state change in the SMBH accretion disk, such as a thin disk that transitioned to a hot accretion flow. This is predicted to occur—in analogy with models developed to explain state changes in X-ray binaries—when the mass accretion rate falls below a few percent of the Eddington rate (e.g., Tchekhovskoy et al. 2014). The optical light curve of AT2018hyz extends to ≈800 days (Gomez et al. 2020; Short et al. 2020; Hammerstein et al. 2022), overlapping the time at which we estimate the outflow was launched. When examining the data from ZTF at later times, we find there is an excess (Section 2). Based on analytical modeling of the optical data using the MOSFiT code (Gomez et al. 2020) we estimate that the mass accretion rate at the time the radio outflow was launched is ≈$0.05{\dot{M}}_{\mathrm{Edd}}$, consistent with the possibility of a state change.

An alternative explanation is that powering a relativistic jet via the Blandford–Znajek process requires a strong magnetic flux threading the black hole horizon. The original magnetic field of the disrupted star in a TDE is not expected to contain a strong enough magnetic field to power a relativistic jet (Giannios & Metzger 2011), which requires an alternative origin. The first possibility is that the magnetic flux could be generated through a dynamo by the accretion disk itself; Liska et al. (2020) found that it may take only ∼10 days to generate poloidal flux from the toroidal field through a dynamo effect once the disk is sufficiently thick, thereby connecting jet production to the disk getting thinner as the accretion rate drops. Alternatively, Tchekhovskoy et al. (2014) and Kelley et al. (2014) suggest that the required magnetic flux may originate from a preexisting AGN disk, which is "lassoed" in by the infalling fallback debris; because the matter falling back at later and later times in a TDE reaches larger and larger apocenter radii, depending on the radial profile of the magnetic flux in the preexisting AGN disk, this could delay the jet production.

Another possibility is that the delayed radio emission is due to the timescales for debris circularization and viscous accretion (Hayasaki & Jonker 2021). In this scenario, the first stream–stream collisions produces the optical/UV emission, creating a debris-circularized ring. The ring evolves viscous-diffusively and reaches the innermost stable circular orbit on a timescale of months to years after the initial flare, with a disk wind velocity of ∼0.4c, consistent with the timescale and velocities we infer for AT2018hyz.

It is also possible that instead the jet has been present for the entire duration of the TDE. However, due to the combination of the high density of the large cloud of circularizing TDE debris (e.g., Bonnerot et al. 2022) and the potential for jet precession (e.g., due to misalignment of the disk angular momentum relative to the black hole spin axis; e.g., Stone & Loeb 2012), the jet is initially choked. At later times, as the accretion rate and gas density surrounding the black hole drop, eventually the jet is able to propagate through the debris cloud and escape.

Planned additional monitoring of AT2018hyz may more clearly elucidate the mechanism responsible for the delayed launch.