AT2019azh: an unusually long-lived, radio-bright thermal tidal disruption event

Tidal disruption events (TDEs) occur when a star is destroyed by a supermassive black hole at the center of a galaxy, temporarily increasing the accretion rate onto the black hole and producing a bright flare across the electromagnetic spectrum. Radio observations of TDEs trace outflows and jets that may be produced. Radio detections of the outflows from TDEs are uncommon, with only about one third of TDEs discovered to date having published radio detections. Here we present over two years of comprehensive, multi-radio frequency monitoring observations of the tidal disruption event AT2019azh taken with the Very Large Array (VLA) and MeerKAT radio telescopes from approximately 10 days pre-optical peak to 810 days post-optical peak. AT2019azh shows unusual radio emission for a thermal TDE, as it brightened very slowly over two years, and showed fluctuations in the synchrotron energy index of the optically thin synchrotron emission from 450 days post-disruption. Based on the radio properties, we deduce that the outflow in this event is likely non-relativistic and could be explained by a spherical outflow arising from self-stream intersections, or a mildly collimated outflow from accretion onto the supermassive black hole. This data-set provides a significant contribution to the observational database of outflows from TDEs, including the earliest radio detection of a non-relativistic TDE to date, relative to the optical discovery.


INTRODUCTION
The central mass in a galaxy influences the dynamics and spatial distributions of stars in the inner regions.In some rare cases a supermassive black hole (SMBH) can capture and destroy a star if it passes within the radius at which the tidal shear forces on the star from the black hole exceed the star's self-gravity (Rees 1988).Such tidal disruption events (TDEs) produce bright flares across the electromagnetic spectrum that are usually visible for timescales of 1-2 years, with approximately half of the debris remaining in orbits bound to the black hole, and other parts flung out on hyperbolic orbits with large velocities (e.g.Lacy et al. 1982;Rees 1988;Evans & Kochanek 1989;Lodato et al. 2009).The bound stellar debris may circularise and form an accretion disk (e.g.Shiokawa et al. 2015;Bonnerot & Lu 2020;Bonnerot et al. 2016;Liptai et al. 2019;Hayasaki et al. 2016; Mummery & Balbus 2020), producing X-ray emission from the accretion stream as the material falls back towards the black hole, and optical emission possibly from the re-processing of the X-rays in the disk (e.g.Auchettl et al. 2017;Cannizzo et al. 1990;van Velzen et al. 2020;Gezari 2021;Strubbe & Quataert 2009;Metzger & Stone 2016;Roth et al. 2016).The disk circularisation time, and thus the time taken for disk formation and accretion onto the SMBH to begin, is a subject of debate (Bonnerot & Stone 2021).Some models predict varying circularisation timescales depending on the physical properties of the star, SMBH, and system (e.g.Lu & Bonnerot 2020;Liptai et al. 2019;Hayasaki et al. 2016;Bonnerot & Stone 2021).X-ray observations of TDEs that trace accretion onto the SMBH have shown variable behaviour, with some TDEs showing bright X-ray emission early on (e.g.Miller et al. 2015), and others showing delayed X-ray flares (e.g.Hinkle et al. 2021).In some cases radio emission is also observed from outflowing material.
Radio observations of TDEs trace the outflows produced by the debris from the destroyed star that may be ejected from the black hole (see Alexander et al. 2020, for a review), enabling detailed insight into the launching of outflows from SMBHs, the circumnuclear density, and how such outflows might provide feedback and influence the evolution of the host galaxy.The radio properties of the TDEs observed to date are diverse, with some exhibiting high luminosity emission (νLν > 10 40 erg/s) that is well-described by a relativistic jet (e.g.Swift J164449.3+573451(Sw J1644+57), Levan et al. 2011;Burrows et al. 2011;Zauderer et al. 2011;Bloom et al. 2011) and others exhibiting lower-luminosity emission (νLν < 10 40 erg/s) that could be described by synchrotron emission from a non-relativistic spherical or mildly collimated outflow (e.g.ASASSN-14li, van Velzen et al. 2016;Alexander et al. 2016).The time relative to the optical/X-ray flare at which the radio emission is observed could provide insight into the mechanism launching the outflow, the star that was destroyed, and the nature of its orbit around the black hole.
Recently it has been suggested that delayed radio emission may be common in TDEs (Horesh et al. 2021a), based on latetime radio flares observed for the thermal TDEs ASASSN-15oi and iPTF16fnl.However, radio observations of TDEs at early times are uncommon, and the apparent lack of detected early-time radio emission could naturally result from a lack of early-time observations, as we demonstrate through earlytime radio observations of the thermal TDE AT2019azh in this work.
Due to the diverse properties of the TDEs that have been observed in the radio to date, there is no consistent explanation for the type of outflows that may be produced in any single event.In some rare cases, as in Sw J1644+57, a relativistic jet is produced with energy ∼ 10 51 erg (Burrows et al. 2011;Brown et al. 2017;Pasham et al. 2015;Cenko et al. 2012;Zauderer et al. 2011).These arise in TDEs that present a non-thermal X-ray spectrum, and the radio emission can be well described by a relativistic jet model in which synchrotron emission is produced as the jet shocks the circumnuclear medium (CNM) and slows, producing emission similar to a gamma-ray burst (e.g.Metzger et al. 2012;Kumar et al. 2013).In other cases, where the TDEs exhibit a thermal X-ray spectrum (as for example in ASASSN-14li and AT2019dsg), a less energetic (E ∼ 10 46 − 10 50 erg), nonrelativistic outflow is produced.
There are a few possible scenarios for producing the observed properties of these non-relativistic outflows, some of which are difficult to rule out with the current set of observations.The prevalent models include a disk wind model in which the outflow is produced early on by accretion onto the SMBH, and emits synchrotron radiation as it moves through the interstellar medium around the black hole (e.g.Alexander et al. 2016).Alternatively, in a similar scenario, a mildly collimated, non-relativistic jet could be produced by the accretion onto the SMBH (e.g.van Velzen et al. 2016).In this case, the jet, emitting radio emission by an internal emission mechanism, could switch on at later times and a constant injection of energy could be observed, with the energy increasing with time (e.g.Falcke & Biermann 1995).Another possibility is a collision-induced outflow, in which the debris from the destroyed star undergoes stream-stream collisions as it circularises into an accretion disk, with a significant amount of gas becoming unbound and ejected in an approxi-mately spherical outflow (Lu & Bonnerot 2020).Finally, the radio emission could also be produced by the unbound tidal debris stream, which would be ejected from the system with escape velocities ∼ 10 4 km/s in a concentrated cone close to the orbital plane (Krolik et al. 2016).
Recently there has been an increase in the number of TDEs with radio detections (e.g.Alexander et al. 2016;van Velzen et al. 2016;Alexander et al. 2017;Cendes et al. 2021;Horesh et al. 2021b), with the next few years expected to bring a large number of new radio observations of these unique events due to targeted radio campaigns to follow-up optical and Xray detected events.These new observations will be crucial in characterising the mechanism behind the radio-emitting outflows that can be produced, and identifying if there is a single mechanism behind all radio outflows, or if the type of outflow differs between individual systems.
In this work we present over two years of radio monitoring observations of the thermal TDE AT2019azh.AT2019azh was first discovered on 2019 February 22 by the All-Sky Automated Survey for Supernovae (ASASSN) (Brimacombe et al. 2019) and named .It was detected by the Zwicky Transient Facility (ZTF) on 2019 February 12 and denoted ZTFaaazdba (van Velzen et al. 2019).The source was coincident with the nucleus of the E+A galaxy KUG 0810+227, with a redshift of z = 0.022 (luminosity distance of 96 Mpc).Spectra obtained by the Nordic Optical Telescope Unbiased Transient Survey (NUTS) on 2019 February 22 (Heikkila et al. 2019), ePESSTO on 2019 February 25 (Barbarino et al. 2019), and the Spectral Energy Distribution Machine mounted on the Palomar 60-in telescope on 2019 February 24 and March 10 (van Velzen et al. 2019) all revealed a blue, featureless spectrum with narrow emission and Balmer absorption features associated with the host galaxy.The event was also detected in X-ray and UV with the Neils Gehrels Swift Observatory and was found to have a high Xray blackbody temperature of kT = 0.06 eV assuming a thermal spectrum (van Velzen et al. 2019).The combination of optical spectral properties, high blackbody temperature, location at center of host galaxy, and lack of spectroscopic AGN or supernova features led van Velzen et al. (2019) to classify the source as a tidal disruption event.The first reported radio detection of the event was by Perez-Torres et al. (2019) with the electronic Multi-Element Remotely Linked Interferometer Network (e-MERLIN) at 5 GHz on 2019 May 21 and 2019 June 11.In this work we present an earlier radio detection, on 2019 March 9. Hinkle et al. (2021) analysed the optical, UV, and X-ray observations of AT2019azh from 30 d before to ∼300 d after the optical peak.During the first 200 d there was very low-level X-ray emission with a harder spectral index, and strong optical/UV emission that evolved as expected for a thermal TDE.The X-rays brightened by a factor of 30-100 approximately 250 d after discovery, and the X-ray spectrum became softer.The optical/UV flare was observed to begin on MJD 58528 (2019 Feb 14) and peaked on approximately MJD 58560 (2019Mar 18, Hinkle et al. 2021).From a powerlaw fit of the optical rise observed by ASASSN, Hinkle et al. (2021) inferred that the time of first light for the event was MJD 58522 (2019 Feb 8).Liu et al. (2019) found X-ray flaring episodes during the early times, which were temporally uncorrelated with the optical/UV emission.They deduced that the optical and X-ray data could be explained by a two-process scenario in which the early emission in UV/optical is explained by emission from debris stream-stream collisions as the bound debris is becoming circularised, and the low-level early X-ray emission is due to a low-mass accretion disk forming during this time.Liu et al. (2019) explain the late X-ray brightening as due to the major body of the disk forming after circularisation of the bound stellar debris.However, Hinkle et al. (2021) concluded that the late-time X-ray brightening was a consequence of an increase in the area of the X-ray emitting region via the blackbody radius, while the short term X-ray variability was due to changes in the X-ray temperature.In this scenario, the late-time X-ray brightening is not due to delayed accretion disk formation, but rather an expansion of the X-ray emitting region.However, Mummery (2021b) recently showed that the X-ray blackbody radius is not a good measure of length scales in a TDE system, implying that the change in blackbody radius may not be caused by a change in the disk radius.The nature of AT2019azh is thus a subject of debate, with evidence both for and against a delayed accretion scenario to explain the multiwavelength emission from the event.
In this work we present 13 epochs of radio observations of AT2019azh taken with the Very Large Array (VLA) and MeerKAT beginning 2019 March 9 (before the optical peak) and spanning until 2021 June 5.These radio observations enable further insight into the nature of the TDE, including the disk formation and launching of the outflow.The paper is outlined as follows: in Section 2 we describe the radio observations and data processing.In Section 3 we present the radio spectral observations and synchrotron emission fits.In Section 4 we describe the modelling of the radio emission to predict physical properties of the outflow.In Section 5 we present a more detailed accretion-disk model of the multiwavelength observations.In Section 6 we discuss the implications of these results, the possible nature of the outflow in AT2019azh, and relate the outflow properties to those of other TDEs.Finally in Section 7 we provide a summary of our results and concluding remarks.

VLA observations
We obtained radio observations of AT2019azh with the NRAO's Karl G. Janksy Very Large Array (VLA) spanning from 2019 March 9 to 2021 February 26 across 300 MHz-24 GHz (P-to .In our first observation, on 2019 March 9, we observed the optical position of the source (RA, Dec) 08:13:16.945,+22:38:54.03and detected faint radio emission at 10 GHz with a flux density of 150 ± 12 µJy.The position of this radio emission was (RA, Dec) 8:13:16.95,+22.38.54.02 with a positional accuracy of 1 arcsecond, coincident with the optical position.We subsequently triggered follow-up observations over a broader frequency range, and continued to monitor the source evolution over the following two and a half years, taking 13 epochs of observations in total.The observations are summarised in Table 2.2 (flux densities and frequencies are available in the online machine-readable format of the table ).
All data were reduced in the Common Astronomy Software Application package (CASA 5.6, McMullin et al. 2007) using standard procedures.Where possible, we calibrated the data using the VLA calibration pipeline available in CASA.In all observations 3C 147 was used as the flux density calibrator.For phase calibration we used ICRF J082324.7+222303 for 12-26 GHz (K-and Ku-band), ICRF J083216.0+183212 for 4-12 GHz (X-and C-band); ICRF J084205.0+183540 for 1-4 GHz (S-and L-band); and PKS J0801+1414 for 0.23-0.47GHz (P-band).The P-band data were reduced manually using standard procedures in CASA, including phase and amplitude self-calibration.Images of the target field of view were created using the CASA tasks CLEAN or TCLEAN (for epochs post April 2019) for all bands except P-band, where we used the WSCLEAN (w-stacking CLEAN) imager (Offringa et al. 2014;Offringa & Smirnov 2017).The source flux density was measured in the image plane, by fitting an elliptical Gaussian fixed to the size of the synthesized beam using the CASA task IMFIT.The errors associated with the measured flux densities include a statistical uncertainty and a systematic one due to the uncertainty in the flux-density bootstrapping, estimated at 5%.Where enough bandwidth was available, we split the L-, C-, and S-band data into four sub-bands when imaging, and the X-band data into two sub-bands.The source was detected at a 4-σ confidence level in the P-band observations, so we did not split these data into sub-bands.

MeerKAT observations
We also observed AT2019azh with MeerKAT on four occasions between 2019 November 29 and 2020 November 14.We used the 4K (4096-channel) wideband continuum mode and observed with bandwidth of 856 MHz around a central frequency of 1.28 GHz.Each observation was about 2 h long in total, except for that on 2020 Nov. 14, which was only 1 h long.
The data were reduced using the OxKAT scripts (Heywood 2020).We used PKS J0408−6544 (QSO B0408−65) to set the flux density scale and calibrate the bandpass, and ICRF J084205.0+183540 as a secondary calibrator.The final images were made using the WSCLEAN imager (Offringa et al. 2014;Offringa & Smirnov 2017), and resolved into 8 layers in frequency.WSCLEAN deconvolves the 8 frequency layers together by fitting a polynomial in frequency to the brightness in the 8 frequency-layers.Our flux densities include both the statistical uncertainty and a systematic one due to the uncertainty in the flux-density bootstrapping, estimated at 5%.
To ensure no systematic offset between epochs and instruments, in Appendix A we present an analysis of flux density measurements of three background sources for 9 epochs of VLA data and 4 epochs of MeerKAT data.We found no significant systematic offset between the two instruments, and found flux densities between VLA epochs were consistent to within ∼ 10%.The flux scale obtained through calibration of the VLA data is consistent across epochs to within a few percent, indicating that the flux density fluctuations we infer between epochs are larger than that expected through calibration differences alone.However, there is no systematic frequency dependence for these inter-epoch flux density variations, and these differences between epochs could be due to intrinsic variability of the background sources, which are expected to be variable at some level.

Multiwavelength observations
We obtained forced point-spread function fitting (PSF) photometry of AT2019azh from the public ZTF MSIP data through the ZTF forced-photometry service (Masci et al. 2019).We filtered the resulting optical light curves for observations impacted by bad pixels, and required thresholds for the signal-to-noise of the observations, seeing, the sigmaper-pixel in the input science image, and several parameters relating to the photometric and astrometric calibrators.
The majority of the Swift UVOT observations were published in van Velzen et al. (2021).Here, we include new observations taken after the publication of that work.We used the uvotsource package to analyze the Swift UVOT photometry and the resulting UV data have been host galaxy subtracted.We also include NICER and XMM-Newton observations reported in Hinkle et al. (2021).

RADIO LIGHTCURVE AND SPECTRA
The 2.25, 5, and 9 GHz VLA lightcurves for AT2019azh are plotted in Figure 1, as well as a comparison of the 5 GHz lightcurve with other thermal TDE lightcurves.The radio emission from AT2019azh rose relatively slowly at all radio wavelengths until approximately 625 d post optical discovery, at which time the higher frequency (> 4 GHz) emission started to decrease while the 2 GHz emission remained relatively constant.Such a slow rise in the radio relative to the optical peak, which occurred around the time of our first radio detection, places AT2019azh in the slow-rising thermal TDE population (Figure 1).In contrast, some thermal TDEs have been observed to begin fading in the radio soon after the optical peak (e.g.Alexander et al. 2016;Horesh et al. 2021b).
The 5 GHz luminosity of AT2019azh increases approximately linearly with time, similar to that of the relativistic event Sw J1644+57.However, Sw J1644+57 rose to a peak within ∼100 d (Eftekhari et al. 2018), whereas AT2019azh took over ∼600 d.We note that AT2019azh was detected in the radio significantly earlier relative to the optical peak than the other thermal TDEs, and a similar slow rise cannot be ruled out for ASSASN-14li, CNSS J0019+00, or XMMSL1 J0740-85.The rise observed for AT2019azh is significantly different than those of ASASSN-15oi, which had early radio non-detections (Horesh et al. 2021b)), and AT2019dsg, which rose to a peak over < 350 d with L ∝ t 2.5 (Stein et al. 2021).

Spectral fitting
We fit the observed radio spectrum at each epoch using the synchrotron spectrum model described in Granot & Sari (2002).The flux density of the synchrotron emission spectrum, assuming νm < νa < νc (where νm is the synchrotron minimum frequency, νa is the synchrotron self-absorption frequency, and νc is the synchrotron cooling frequency), is described by where ν is the frequency, Fν,ext is the normalisation, s1 = 3.63p − 1.60, s2 = 1.25 − 0.18p, β1 = 5 2 , β2 = 1−p 2 , and p is the synchrotron energy index.
Archival Faint Images of the Radio Sky at Twenty-cm (FIRST) survey observations from January 1996 of the host galaxy of AT2019azh show no detection at 1.4 GHz, with a 3 σ upper limit of 0.41 mJy (Becker et al. 1995).These observations place a strong upper bound on the host galaxy contribution to the radio emission and make recent AGN activity in the galaxy unlikely.Whilst Hinkle et al. (2021) were also able to rule out strong AGN activity based on optical properties of the host galaxy, they found that optical and X-ray observations cannot rule out the presence of a low-luminosity AGN in the host galaxy, KUG 0810+227.In order to account for 58600 58800 59000 59200 59400 Time (MJD)  the possibility of some low-level contribution from the host galaxy to the observed radio flux densities of the outflow, we add a host component to the spectral fitting that is described by where F0 is the flux density measured at 1.4 GHz (F0 < 0.41 mJy) and α0 is the spectral index of the host galaxy.
The total observed flux is then given by We fit the spectra for all epochs using a Python implementation of Markov Chain Monte Carlo (MCMC), emcee (Foreman-Mackey et al. 2013).We use a Gaussian likelihood function where the variance is underestimated by some fractional amount f .We assume flat prior distributions for all parameters, and allow p to fall in the range 2.5 − 4.0.Whilst it is common practice to fix p between epochs (e.g.Alexander et al. 2016;Cendes et al. 2021), there is evidence in the observations that the spectral slope is not constant for this event (see Figure 2), so we do not fix p in our modelling.To constrain the host flux density and spectral index (Equation 2), we first ran a MCMC fit for epoch 11 (2020 Dec 12) only, where the synchrotron spectrum is very well constrained by the observations, ensuring that F0 < 0.41 mJy and −2 < α0 < 2. We found the best solution for F0 = 0.175 mJy and α0 = −0.84.Next we fit the total flux density using the determined host contribution as a function of frequency for the energy index, p, the flux normalisation, Fν,ext, the minimum frequency, νm, and the self absorption frequency, νa for all epochs using Equation 3.For the firs three epochs the peak frequencies and flux densities are not well-constrained due to the lack of low-frequency radio coverage.The MCMC spectral fitting results for the synchrotron self-absorption break and peak flux density for these epochs are dependent on the choice of prior.Thus, for these epochs we provide upper and lower limits respectively for νa and F peak .
The spectral fits for each epoch are plotted in Figure 2, and the best fit peak flux densities and frequencies from the spectral fits are reported in Table 2 and plotted in Figure 3.For the epochs where the peak of the synchrotron spectrum is well constrained, we find the peak frequency, ν peak , remained relatively constant at ν peak = 1.1 ± 0.3 GHz, with a slight downwards trend over the 800 d spanned by our observations, whilst the peak flux density, F peak , increased approximately linearly with time and only showed signs of decreasing in the final epoch, 820 d post-disruption.The index of the electron energy distribution, p, remains roughly constant at p ≈ 2.7 ± 0.2, similar to that of other thermal events (e.g.Alexander et al. 2016;Cendes et al. 2021;Stein et al. 2021;Horesh et al. 2021b) excepting three epochs at 254, 459, and 749 d post disruption at which the energy index shows a significant steepening to p ≈3, 3.7, and 3.3 respectively.We note that the slight spectral steepening at 255 d is likely not real based on an analysis of the flux density of background sources in the field for the different epochs presented in Appendix A and that the significance of the change in p is 3-σ and 2-σ respectively for the other two epochs.We found no evidence in the background source measurements for inconsistent calibration with frequency for the two other epochs where we observed a steepening, although there are small systematic flux density offsets between most epochs (see Appendix A), which we account for with the added systematic uncertainty of 5% to the flux densities.We thus conclude that the steepenings we observed are real for May 2020 (epoch 10) and February 2021 (epoch 12), and are not artefacts of inconsistent calibration over the frequency ranges for those epochs.We note that the 10% flux density offsets between epochs on the background sources suggests a possible systematic uncertainty of 10% in the flux density calibration.Such a systematic un-certainty would affect our peak flux density values, however, the uncertainty of the peak flux density is dominated by the peak not being well constrained in many of the epochs due to the paucity of the data at low frequencies.

MODELLING OF THE RADIO EMISSION
We model the radio emission from the outflow using the standard synchrotron emission model outlined in Barniol Duran et al. (2013), in which the ambient electrons are accelerated into a power-law distribution by the blastwave from the outflow, N (γ) ∝ γ −p , where γ is the electron Lorentz factor (γ ≥ γm, where γm is the minimum Lorentz factor) and p is the synchrotron energy index.We assume equipartition between the electron and magnetic field energy densities in order to derive the equipartition energy and radius.Although the emitting region may not be in equipartition, we can derive estimates of the physical system parameters by parameterising the deviation from equipartition and accounting for its effect.This approach allows us to estimate key physical quantities such as the ambient electron density, magnetic field strength, mass of the emitting region, and velocity of the ejecta.We assume the fraction of the total energy in the magnetic field is 0.1%, B = 10 −3 , based on observations of other TDEs and supernovae (e.g.Eftekhari et al. 2018;Horesh et al. 2013).We assume that 10% of the total energy is carried by the electrons, e = 0.1, given that electrons are typically accelerated much less efficiently than protons in astrophysical accelerators (e.g.Morlino & Caprioli 2012).We find no evidence for a relativistic outflow (see Section 6.2), and assume the outflow is non-relativistic (i.e. the bulk Lorentz factor, Γ = 1), and that the peak of the radio spectrum is associated with the synchrotron self-absorption frequency (i.e.νa = νp).We model two different geometries, one where the emitting region is approximately spherical (with geometric factors1 fA = 1 and fV = 4/3), and one where the emitting region is conical (with geometric factors fA = 0.13 and fV = 1.15), corresponding to a mildly collimated outflow with a half-opening angle of 30 degrees.
In the Newtonian regime, the equipartition energy, corresponding to the minimum total energy in the observed region, assuming νa > νm, is given by (Barniol Duran et al. 2013) where d is the distance from the observer, z is the redshift, χe = p−2 p−1 e mp me (me is the electron mass and mp is the proton mass), or χe = 2 if Γ = 1 (Newtonian case), and The equipartition radius is given by (5) To infer the physical system parameters we first correct the energy and radius for the system being out of equipartition using the following assumptions where = B e 11 6 , i.e. the total energy is minimised with respect to R at Req when the energy in the magnetic field is 6/11 the energy in the electrons, so the deviation from equipartition is parameterised by (Barniol Duran et al. 2013).
Using the corrected radius, R and corrected energy, E, we then calculate the total number of electrons in the observed region (Ne), ambient electron density (ne), magnetic field (B), mass of the emitting region (Mej), and outflow velocity (β = v/c), using Equations 8-13 (Barniol Duran et al. 2013), as where γm γa 1−p is a correction factor to account for the regime where νm < νa and the extra factor of 4 arises due to the Newtonian correction; γm = 2, and γa is given by The ambient electron density is then inferred via We determine the velocity of the outflow, β, by rearranging Equation 22 from Barniol Duran et al. (2013), where the observer time, t is given by t where we set the time t relative to the approximate launch date of the outflow, MJD 58522; inferred based on a linear fit to the predicted radius and estimated optical time of first light (see below).The magnetic field is given by ( and the approximate mass in the emitting region of the ejecta is given by noting that this is a lower limit on the mass in the emitting region of the outflow due to the energy estimate from equipartition also being a lower limit on the energy.
The physical outflow properties as predicted by these equations are listed in Table 2 and plotted in Figure 4.All uncertainties reported correspond to the 1σ uncertainty obtained through propagating the uncertainty of ν peak , F peak and p obtained through the MCMC modelling of the observed spectra.
The derived radius of the outflow increases with time, following the relation R ∝ t 0.65 (reduced χ 2 =1.79), or R ∝ t (reduced χ 2 =2.67).We thus deduce that the outflow is roughly undergoing free expansion, with the velocity remain- ing approximately constant at β ≈ 0.2 ± 0.1 (conical) or β ≈ 0.1 ± 0.06 (spherical) until the final epoch in which the velocity shows a slight decrease.We note that the powerlaw fit to R, indicating a decelerating outflow, is statistically preferred to the constant velocity case.However, given the underlying assumptions of the radius calculation, that the synchrotron peak flux density and frequency were not resolved by the observations for the first few epochs, and that outflows from other thermal TDEs have all been observed to undergo approximately free-expansion at early times (e.g.Alexander et al. 2016;Cendes et al. 2021), in the sections that follow we operate under the assumption that the outflow was freely expanding with ∼constant velocity until at least the radio lightcurve peak (t ≈ 650 d).
The derived energy of the outflow increases approximately linearly with time for all observations, but also seems to show statistically significant non-uniform fluctuations with time.The magnetic field shows a slight decrease with time and the inferred mass in the emitting region of the outflow also increases with time (based on the energy prediction).
A simple linear fit (assuming constant velocity) to the predicted radius gives an outflow launch date of MJD 58435±10 (2018 Nov 13, concial) or MJD 58432±10 (2018 Nov 10, spherical), approximately 120 d before the first radio detection on MJD 58551 (2019 Mar 09).The optical/UV flare was observed to begin on MJD 58528 (2019 Feb 14), ∼ 90 d after the estimated outflow launch date, and peaked on approximately MJD 58560 (2019 Mar 18 Liu et al. 2019;Hinkle et al. 2021).From a power-law fit of the optical rise observed by ASASSN, Hinkle et al. (2021) inferred that the time of first light for the event was MJD 58522 (2019 Feb 8), indicating the radio outflow was likely launched later than MJD 58433 with an initial velocity higher than 0.1 c.Thus in this work we assume an outflow launch date of MJD 58522; coincident with the optical time of first light.

Expected future evolution of the outflow in the
Sedov-Taylor decay phase The predicted evolution with time of the outflow's velocity (Figure 4) is consistent with either an outflow expanding with constant velocity until the last three epochs (post-radio peak), or a gradually decelerating outflow.In the case that the outflow had approximately constant velocity of ≈ 0.1 c until the last three epochs, we suggest that this could be indicative of an outflow that was "coasting" until the peak radio flux was reached (t ≈ 650 d), and is now decelerating as the flux decays.However, the data could be equally well explained by a model in which the outflow consistently decelerates over the course of the observations.Under the assumption that the outflow did exhibit a coasting phase, the outflow sweeps up material from the circumnuclear medium (CNM), increasing the energy released in the emitting region as the outflow impacts the CNM (e.g.Generozov et al. 2017).The onset of deceleration could indicate the outflow is entering the Sedov-Taylor phase of its evolution, at which time the outflow has reached peak energy/flux emission in the free expansion phase by sweeping up mass from the CNM, and begins to decelerate with constant energy as the blastwave approaches spherical symmetry (e.g.Sironi & Giannios 2013).In Figure 4 there is some evidence for deceleration of the outflow in the last three epochs of observations from 660 d post-disruption, corresponding to the epochs post-radio luminosity peak, with the deceleration most evident in the final epoch.We note, however, that an alternate velocity evolution consisting of a power-law evolution in radius and velocity with time fits the radius measurements equally well.
Under the assumption that the outflow velocity only began decelerating after the lightcurve peaked, the Sedov-Taylor phase enables predictions about the rate of decay of the emission for observing frequencies much above the self-absorption frequency and much below the cooling frequency for a stratified CNM density profile and the synchrotron spectrum power law index p (Sironi & Giannios 2013).In Figure 5 we show the  Sedov-Taylor solution for different CNM density stratifications with p = 3 (suitable for the final epoch of observations) for the 5.5 GHz light curve of AT2019azh.The predicted flux evolution is calculated assuming late-time radio emission in the Sedov-Taylor phase, in which a spherical shock runs into a stratified medium with density profile n ∝ r −k and the electrons in the shock are accelerated into a synchrotron powerlaw distribution (Sironi & Giannios 2013).The best-fit solution is for a steep CNM density gradient, with k = 2.5, similar to the density gradients observed for ASASSN-14li and AT2019dsg (Figure 7).The current decay-rate of the radio emission at 5.5 GHz is well-fit by the Sedov-Taylor approximation for a CNM density n ∝ −2.5 and p = 3.

MULTIWAVELENGTH MODELLING: ACCRETION DISK EVOLUTION
To assess the possibility that the radio outflow from AT2019azh was produced by accretion onto the SMBH, in this section we model the accretion disk emission based on the optical, UV, and X-ray properties of the event.The evolving disk density profiles are calculated using the method developed in Mummery & Balbus (2020) and Mummery (2021a), to which the reader should refer for detailed information.We assume that the black hole environment is initially completely devoid of material, before feeding disk material into a ring according to the following relationship: where δ is the Dirac delta function and ∆t ensures the feeding rate is finite as t → 0, and was taken equal to one code time step.The feeding radius was taken to be r0 = 50Rg, appropriate for a TDE around a low-mass black hole like AT2019azh, and the disk evolution was started at a time corresponding to the first observed optical emission of the source.This model assumes that the matter is fed into the disk at the rate at which disrupted stellar material returns to pericentre (the so-called 'fall-back' rate Ṁfb , Rees 1988).The material fed into the disk then evolves according to the equations of disk angular momentum conservation and disk mass conservation.Energy conservation then allows the disk temperature to be calculated at each radius and time.We assume that each disk radius emits like a colour-corrected blackbody (using the Done et al. ( 2012) colour-correction model), and ray-trace the resulting disk emission profile.We include all relevant relativistic effects, such as Doppler and gravitational redshift, and gravitational lensing.The evolution of the X-ray and UV light curves of a thermal TDE are therefore determined by three fitting parameters: the black hole mass M , the total accreted mass Macc (a normalisation on the source term, Eq. 14), and the viscous timescale of the evolving disk tvisc.We compute Macc in the following manner where Ṁ (rI , t) is the inner-most stable circular orbit (ISCO) mass accretion rate.
Given the simplifications applied in the modelling (the black hole spin is fixed to zero, and the observer inclination angle is fixed to θ obs = 30 • ) the best fit parameter values and their associated uncertainties should be treated with some caution, although Mummery & Balbus (2020) found only a moderate change (factor of ∼ 1.5) in best fit black hole mass over a wide range of black hole spins.
Finally, we determine the best fit system parameters by simultaneously minimizing the sum of the squared differences between the model and the different UV and X-ray light curves of AT2019azh.As in previous works, we anticipate large formal values of the reduced χ 2 .These large values result as a consequence of short time-scale fluctuations present in the well-sampled TDE light curves, and are to be expected in any theoretical model using a smooth functional form for the time-dependence of the turbulent stress tensor.Our standard approach implicitly averages over rapid turbulent variations.Short time-scale fluctuations are likely to be highly correlated so that accurately assessing the statistical significance of the fit is not straightforward.We have therefore used χ 2 minimisation as a sensible guide towards finding a best fit.
The observed and modelled optical, UV, and X-ray lightcurves for AT2019azh are plotted in Figure 6 for two scenarios: early accretion and delayed accretion.The observed optical and UV data are well fit by both the early and delayed accretion models, with little difference between the two.The observed X-ray lightcurve is significantly better fit by the delayed accretion model at early times, whilst the late-time X-ray lightcurve is well-fit by either model.We note that in the case of significant X-ray obscuration at early times the early accretion model may also be viable.

Different disc evolution scenarios
The X-ray evolution of AT2019azh is, as discussed in the Introduction, somewhat atypical for a TDE.Thus, we present The ZTF r-and g-band observations are shown in red and green respectively and the Swift UVM2 filter and 2-10 keV X-ray observations are shown in magenta and black respectively.We model two scenarios: a delayed disk-formation (dashed lines) and an early disk-formation (solid lines).The optical and UV observations are well-fit by either model, while the X-ray observations are better-fit by the delayed accretion scenario, unless in the case of significant X-ray obscuration at early times.The X-ray observations are from Hinkle et al. (2021) and only the 0.3-10 keV data from each telescope was used in the fitting.
two models for the light curve evolution of AT2019azh: 'prompt' and 'delayed' accretion.There are two ways in which the peak X-ray luminosity of a TDE may be delayed.The X-ray luminosity of a TDE is primarily a function of the hottest temperature in the TDEs disc (Mummery & Balbus 2020), as the TDEs peak X-ray luminosity is only reached once the TDEs inner disc density reaches its maximum.As typical TDE feeding radii are tens to hundreds of gravitational radii, a large viscous timescale (which delays the build up of the inner disc density by increasing the length of time it takes for the disc material to propagate inwards) can suppress early TDE X-ray emission.Alternatively, if there is substantial obscuration of the inner regions of the accretion disc at early times, which then clears at larger times, this may also lead to a late-time rise in X-ray luminosity.We model both scenarios in this section.In the delayed accretion model we fit to the entire X-ray light curve, finding as expected a large viscous timescale tvisc = 220 ± 20 d.The accreted mass, Macc = 0.1 ± 0.02M , is consistent with a star of stellar mass M 0.2M (i.e., no missing energy problem, as in Mummery 2021b).The best fitting black hole mass MBH = 3.2 +0.15 −0.1 × 10 6 M , is consistent with the galaxy scaling measurement MBH < 4 × 10 6 M (van Velzen et al. 2016).The peak Eddington ratio of this model is l peak = 0.22 (L peak = 8.5 × 10 43 erg/s), and the total radiated energy was E rad 8 × 10 51 erg.In the prompt accretion model, which requires substantial early time obscuration, we find a viscous timescale tvisc = 65 ± 10 d.The accreted mass, Macc = 0.07 ± 0.01M , and best fitting black hole mass MBH = 2.1 +0.15 −0.1 × 10 6 M are again consistent with their expected values.The peak Eddington ratio of this model is l peak = 0.75 (L peak = 2.1 × 10 44 erg/s), and the total radiated energy was E rad 5 × 10 51 erg.The difference between the two accretion evolution scenarios can be seen in Figure 6, which shows the accreted mass as a function of time (plotted in the same units as Figure 4).

DISCUSSION
Our findings reveal a likely non-relativistic outflow with constant (or gradually decreasing) velocity and continuous kinetic energy increase during the radio rise (up to 666 d), and constant energy post-radio peak (from 666-849 d).We infer that the outflow ranges from radii of ∼ 3×10 16 cm-2×10 17 cm with energies of ∼ 3 × 10 47 erg-1 × 10 51 erg.These energy and radii correspond to a magnetic field of ∼ 0.05 G and ambient electron density of ∼ 50-3000 cm −3 .The observed energy and radius increased with time until the peak radio luminosity was reached.
The optical and X-ray observations, particularly the latetime X-ray rise, have been explained with either a delayedaccretion disk formation scenario (Liu et al. 2019) or due to an increase in the X-ray emitting region (Hinkle et al. 2021).We modelled two disk emission scenarios in Section 5, and found that the observed UV/optical behaviour of the event was well-fit by either model (Figure 6).The X-ray emission is better fit by a delayed accretion scenario except in the case of significant X-ray obscuration at early times.

The unusual late time steepening of the synchrotron spectrum
We observed statistically significant fluctuations (at a 3-σ level) in the synchrotron energy index, p, of the opticallythin part of the synchrotron spectrum from ∼ 450 d postdisruption.We found that the energy index steepened to p = 3.7 ± 0.2 at t = 460 d, reducing back to p ≈ 2.6 at t = 660 d, steepening again to p = 3.3 ± 0.1 at t = 750 d, and finally reducing back to p ≈ 3 in the final epoch at t = 850 d post-disruption.The mean value of p for the epochs without spectral steepenings is 2.7 ± 0.2 (1-σ error, excluding the first 4 epochs where the spectra were not well constrained), indicating the spectral steepening at t = 460 d is significant to 3-σ, and the steepening at t = 750 d is significant to 2-σ.After detailed investigation (see Appendix A), we found that these fluctuations are not due to calibration issues with the data or inconsistent flux density scaling between epochs as there is no such systematic difference in three background sources that we examined in the field of view for each epoch (excepting the first minor steepening at 200 d, which we conclude is not statistically significant, see Appendix A).Fluctuations in the energy index have not been observed in the radio emission from thermal TDEs to date, and are difficult to explain in the current (single-zone) synchrotron emission model.Usually the steepening of a synchrotron spectrum can be attributed to adiabatic cooling of the electrons, and would indicate the detection of a cooling break in the spectrum (Granot & Sari 2002).However, the adiabatic cooling timescales are too long to explain the fluctuations on timescales of ∼months that we observed in this event, and we find no evidence for the presence of a cooling break in any of the spectra (see Appendix B).
We propose that the energy index fluctuations could be due to a spherical or collimated outflow encountering an inhomogenous CNM, or fluctuations in the energy injection rate of a collimated jet-like outflow.In the spherical outflow scenario, different populations of electrons from different regions of the outflow might encounter inhomogenous clumps of the CNM, changing the emitting properties of different popula-tions of electrons, each with their own synchrotron spectrum.Smaller populations will fade quickly, contributing less to the total radio emission, and allowing the flux to fall at higher frequencies.The synchrotron spectra we observe at each epoch is the sum of the emission from these different populations of electrons.If the fluctuations in p are due to the changes in the shock acceleration efficiency, we expect more fluctuations in the radio lightcurve at 9 GHz than at 5 GHz, as indeed is seen in Figure 1.

The outflow mechanism
With the addition of the radio observations to the multiwavelength data we are able to obtain a more robust picture of the event and how the different types of emission may have been produced.The empirical properties of the outflow in AT2019azh obtained partly from the radio observations are key to modelling and understanding the mechanism that produced the outflow.What is crucial to understanding the event is that the radio outflow was first observed early, around the time of disruption, and well before the X-ray emission brightened, in contrast to the suggestion that delayed radio emission is common in TDEs (Horesh et al. 2021a).Furthermore, the energy index fluctuations observed in the radio emission place strong constraints on the geometry of the outflow; a spherical homogeneous outflow cannot produce the observed fluctuations.Below we discuss different scenarios that could explain the observed properties of the outflow in AT2019azh.

Accretion-driven wind outflow
Accretion onto an SMBH can produce winds and outflows that would be observable in the radio as they travel through the CNM at velocities ∼ 0.01 − 0.1 c (e.g.Mohan et al. 2021;Strubbe & Quataert 2009;Tchekhovskoy et al. 2014).A popular model to explain non-relativistic outflows from TDEs is a spherical wind driven by accretion onto the SMBH (e.g.Alexander et al. 2016;Cendes et al. 2021).In this scenario, the radio outflow should appear at approximately the time that high accretion luminosities are observed at X-ray wavelengths, provided there is no obscuration of the X-ray emission.A wind outflow would have approximately spherical geometry (e.g.Mohan et al. 2021), and produce radio emission from a forward shock from the non-relativistic outflow expanding into the CNM driven by the gas accretion onto the SMBH.
The early radio emission for AT2019azh is difficult to explain as an accretion wind-induced outflow in the delayed accretion scenario, due to the lack of bright X-ray emission indicative of significant accretion (requiring FX ∼ 10 11 erg/s/cm 2 ) and the lack of any X-ray/optical correlation, unless there was strong obscuration of the X-ray emission.In Section 5, we found that the optical, UV, and X-ray observations are well-fit by either delayed accretion or early accretion with significant X-ray obscuration (Figure 6).In the delayed accretion scenario, the radio outflow could not be produced by an accretion-driven wind due to the lack of significant amounts of material reaching the black hole to be ejected into the outflow at early times, when the first radio emission was observed.

Sub-relativistic jet
In the scenario of a mildly relativistic or sub-relativistic jet, a collimated outflow is produced by accretion onto the SMBH and the radio emission could be produced by either a forward shock that the jet drives into the surrounding medium, or internally through shocks inside the jet (e.g.van Velzen et al. 2016), both of which would produce synchrotron emission.A sub-relativistic or mildly relativistic jet was proposed initially for ASASSN-14li (van Velzen et al. 2016) and AT2019dsg (Stein et al. 2021).The main argument against a jet-like outflow relies on the geometric factors, and the level of collimation required to obtain a self-consistent solution for the outflow properties.
We deduce that a relativistic jet explanation for the radio properties of AT2019azh is not supported by the data.Similar to the argument against a relativistic jet provided in Alexander et al. ( 2016), if we introduce an additional parameter, the bulk Lorentz factor (Γ) to the synchrotron equipartition model outlined in Section 4 (Barniol Duran et al. 2013), in order to obtain a self-consistent result where Γ 2 (i.e. a relativistic outflow) requires fA 0.01, i.e. a jet with opening angle of 0.1 deg.Such a small opening angle is not possible for SMBH outflows (e.g.Jorstad et al. 2005) and rules out the possibility of a relativistic jet for the outflow from AT2019azh.
The late-time evolution of a sub-relativistic jet and a mildly relativistic spherical outflow appear very similar at radio frequencies (Nakar & Piran 2011), however, early on an initially on-axis relativistic jet that decelerates to non-relativistic velocities would appear much more energetic (with energies comparable to Sw J1644+57 ∼ 10 52 erg).The luminosity we observed for AT2019azh (L ∼ 10 38 erg/s) disfavours the possibility of an initially relativistic on-axis jet for the early radio emission.With observations of this event spanning the peak of the radio lightcurve, we can also deduce that the outflow is likely non-relativistic due to the observed behaviour of the lightcurve as the outflow transitioned from freely expanding to decelerating.The Doppler factors are no longer important when the outflow begins decelerating (e.g.Sironi & Giannios 2013) so the radius constraint obtained is the true radius of the outflow.If the radio emission was produced by an offaxis relativistic jet, at the time of deceleration it would be emitting isotropically and we would expect to see a flux increase, which is not observed in the radio lightcurve.Under the assumption that the outflow transitioned into the subrelativistic Sedov-Taylor decay phase after the peak of the radio lightcurve (Figure 5), the inferred radius at the time of transition yields an average speed of the outflow that is significantly less than the speed of light.This would further confirm the sub-relativistic nature of this event, in contrast to the assumed relativistic event Arp 299-B AT1, which was found to initially move at relativistic speeds for the first ∼ 760 d (Mattila et al. 2018).Alternatively, in the scenario in which the outflow was constantly decelerating over the course of the radio observations, we cannot rule out initially relativistic speeds of the outflow prior to the first radio detection.
A sub-relativistic jet, with Γ ≈ 1, would not require such extreme collimation of the emission.A sub-relativistic jet may present similarly to our conical geometry model (Table 2 and Figure 4).Such a jet would have slightly larger radii, higher velocity, increased energy and require a lower CNM density to self-consistently explain the observed properties of the emission, than would a spherical outflow.A collimated outflow may also explain the energy index fluctuations in the case of an inhomogenous CNM.As a sub-relativistic jet travels through the CNM, it would sweep up material, slowing down and causing the jet to spread laterally (e.g.van Velzen et al. 2016).In this scenario the emission from the jet is more isotropic due to the Doppler factor close to unity, and a narrow viewing angle is not necessary in order to observe the radio emitting region.A sub-relativistic jet driven by accretion may have continued energy injection until the central engine switches off (e.g.Mohan et al. 2021), which could explain the continuous energy increase observed for AT2019azh (Figure 4).Furthermore, the material ejected from close to the black hole could easily reach velocities and energies as high as we predict for AT2019azh due to energy conservation and the angular momentum available at the inner orbits of the SMBH.However, through our disk modelling (Section 5), we infer that the accretion rate never exceeds 0.2 times the Eddington rate in the delayed accretion scenario, and 0.7 in the early accretion scenario.For the observed radio emission to be explained by an accretion-driven outflow it would require inefficient accretion onto the SMBH, which has not been confirmed in observations of SMBHs.Thus we deduce that a mildly collimated sub-relativistic jet may explain the observed properties of the outflow in AT2019azh, under the condition that the accretion disk formation was not delayed and there was significant X-ray obscuration at early times.

Collision induced outflow
Lu & Bonnerot (2020) modelled the self-intersection of tidal debris streams during a TDE and deduced that during stream-stream collisions of the tidal debris, a significant amount of gas will become unbound and ejected as a collision induced outflow (CIO).This kind of outflow would have kinetic energies between 10 50 -10 52 erg and velocities between 0.01-0.1 c, similar to the properties we infer for the outflow of AT2019azh (Figure 4).
Due to the lack of evidence at early times of significant accretion that could produce an outflow from inefficient accretion onto the SMBH (unless in the case of significant X-ray obscuration), the CIO model is quite promising to explain the radio emission from AT2019azh.A CIO will be launched when the streams intersect, an event that precedes the start of accretion onto the black hole, which could well explain the early radio detection.Liu et al. (2019) proposed that a collision induced outflow model could explain the UV/optical and X-ray emission from AT2019azh at early times as the debris is becoming circularised; our radio detection pre-optical peak supports this theory.
In the case of a CIO, the outflow would be produced by a prompt injection of energy during the circularisation, and the outflow emission would evolve over time as the spherical ejecta is slowly shocked by the CNM and sweeps up material.In order to reach velocities as high as 0.1 c, the pericenter of the destroyed star would need to be within 10-15 Rg of the SMBH (depending on the black hole spin), otherwise the CIO would not be strong enough to reach the observed velocities and energies of AT2019azh.A CIO outflow is well-described by a "coasting", free-expansion and a Sedov-Taylor decay phase once the peak luminosity is reached (Lu & Bonnerot 2020), in contrast to a jet-like outflow which could begin with an expansion phase during which the outflow is powered by the jet-engine, then a Sedov-Taylor decay phase when the jet switches off.The deceleration radius and transition from the free-expansion to Sedov-Taylor phase corresponds to the time at which the lightcurve peaks, and is characterised by E k = (1/2)N (r dec mpv 2 0 (e.g.Lu & Bonnerot 2020), i.e. the deceleration is caused by the outflow interacting with the CNM.

Unbound debris stream
When a star is destroyed by a SMBH, approximately half of the stellar debris will be captured by the gravitational well of the black hole to be accreted, whilst the remaining half of the star is unbound, and may be ejected from the system with high escape velocities (v > 10 4 km s −1 ) (Rees 1988).This unbound debris will interact with the circumnuclear medium, emitting synchrotron radiation in the bow shock that forms along the leading edge of the debris stream (e.g.Krolik et al. 2016).The earliest emitting region will correspond to the fastest unbound debris, expanding at velocities of ∼ 0.05 c (Krolik et al. 2016).Over time, the bulk of the unbound debris will decelerate and eventually become visible, adding to the emitting region of the outflow.The unbound material would be confined to a very small solid angle (Guillochon et al. 2014;Kochanek 1994;Coughlin et al. 2016), which is often used as a justification to rule-out radio emitting outflows being produced by the unbound debris stream (e.g.Alexander et al. 2016, for ASASSN-14li).
The predicted mass in the outflow for the radio-emitting region of AT2019azh is significantly less than expected for the entire unbound debris (Figure 4, assuming a ∼ 1 M star was destroyed and approximately half of the stellar debris is ejected in the unbound debris stream (e.g.Rees 1988)).However, only the debris with the fastest escape velocities would be visible early on, corresponding to a very small fraction of the unbound debris.If the outflow in AT2019azh were produced by the unbound debris we may expect to see the radio evolving at later times as the slower debris catches up.However, in the unbound debris model, the outflow should continue to expand at a constant velocity without slowing down until the ejecta has swept up a mass comparable to its own (Krolik et al. 2016).In Figure 4 we find weak evidence of fluctuations in the velocity of the outflow, and some downwards trend between t = 250 − 850 d.The velocity of the outflow we infer even in the spherical geometry model ( 0.1 c) is higher than expected in models of the unbound debris (≈ 0.05 c) and the inferred energy of the outflow ( 10 48 erg) is also greater than expected for the unbound debris stream (∼ 10 47 erg) (e.g.Krolik et al. 2016).An unbound debris stream outflow would require slightly higher collimation of the emission than we consider in the conical model (fA ≈ 0.2) in Section 4, which would only increase the predicted energy and velocities.
We thus conclude that the multiwavelength emission of the outflow from AT2019azh could be explained self-consistently by either a collision induced outflow or, less likely, an accretion-driven wind or sub-relativistic jet.The equipartition analysis in Section 6.2 provides a robust estimate of the size of the emitting region, and thus its velocity, but it does not enable strong discrimination between the energy source of the emission, and thus the driving source of the outflow.The early radio emission combined with the low X-ray emission and lack of early optical/X-ray correlation points towards the radio outflow not being produced through accretion onto the SMBH, unless there was significant X-ray obscuration.Our disk modelling of the multiwavelength observations in Section 5 indicates that the optical/UV and X-ray emission is well-fit by either an immediate accretion scenario or a delayed accretion scenario with X-ray obscuration.The observed fluctuations of energy and p could be explained by an inhomogenous CNM, in which different populations of electrons in the outflow are emitting and being observed at different times.
A more energetic outflow is possible in a disk-driven outflow than a CIO as the ejected material can be ejected with larger energy, translating to a faster, more energetic outflow, but both scenarios could reach the energies and velocities we deduce for the outflow in AT2019azh.Further, detailed modelling of the multiwavelength properties of this event, building on the detailed optical and X-ray analysis presented in Hinkle et al. (2021), are necessary to truly gain a deep understanding of how to reconcile the X-ray, optical, and radio observations in a self-consistent way.

AT2019azh in the context of other TDEs
On comparison of the outflow properties of AT2019azh with other TDEs with spectrally-resolved radio observations, we find AT2019azh fits well into the population of thermal TDEs (Figure 7).There is only one other thermal TDE (AT2019dsg) with multi-frequency radio coverage at early times (t < 100 d post-distribution).AT2019dsg showed an order of magnitude increase in energy early on, in contrast to the slow rise in radio that we observed for AT2019azh (Figure 1).
We find the ambient density is approximately proportional to R −2.5 for most TDEs, whilst AT2019azh shows significant variation and could be described by n ∝ R −1 to R −2.5 (Figure 7). Figure 7 indicates that for higher ambient densities the radio emission is brighter and peaks earlier (Lu & Bonnerot 2020), as indeed is the case for the light curve of AT2019dsg compared to AT2019azh (Figure 1).
Inhomogeneity of the CNM density for AT2019azh compared to that of other TDEs could possibly explain the fluctuations in energy index we observed that was not seen in the other early-time TDE observations, however other studies of thermal TDE radio spectra do not comprehensively assess the possibility of variations in p.In Figure 7 there is some evidence that the ambient density around AT2019azh is more inhomogenous than the other thermal TDEs based on the inferred spectral properties.The kinetic energy and velocity we infer for AT2019azh are somewhat larger than those for other thermal events, but does not reach the energies (or velocities) observed for relativistic events.
The late-time X-ray brightening of AT2019azh resembles the behaviour observed in the TDEs ASASSN-15oi and OGLE16aaa (Horesh et al. 2021b;Kajava et al. 2020).However, in the case of ASASSN-15oi, there was no radio emission detected early on, in stark contrast to the pre-optical peak radio detection of AT2019azh.Thus, the argument for delayed accretion, while likely for ASASSN-15oi, is difficult to justify for AT2019azh if the outflow was accretion disk-driven.The thermal TDE ASASSN-15oi also exhibited a change in energy index of the optically-thin part of the synchrotron spectrum at late times (Horesh et al. 2021b).ASASSN-15oi exhibited some lower-level radio activity, which faded with time, ∼ 100 d post-disruption, and a drastic increase in radio flux at ∼ 600 d post-disruption (Horesh et al. 2021b).Interestingly, the energy index of the optically thin synchrotron emission became much flatter in the late-time flare compared to the earlier emission.The early emission exhibited a standard energy index of p = 2 − 3, whereas for the later flare it was much flatter at p = 0.2 (Horesh et al. 2021b).Horesh et al. (2021b) deduce that the standard spherical outflow model from a super-Eddington wind cannot explain both the delayed onset of radio emission and late-time flare from this event, which latter would require an outflow to be launched at late-times, and possibly into an inhomogenous circumnuclear medium.They also propose that the emission could be explained by a transition in accretion state onto the SMBH at late-times.However, for AT2019azh we observed the opposite behaviour in p; a steepening of the energy index at later times, followed by fluctuations over the next months.
The late-time radio emission of AT2019azh behaves similarly to that of the TDE AT2019dsg (Stein et al. 2021), which is explained by Cendes et al. (2021) to likely be driven by a spherical outflow from a super-Eddington wind.However, the outflow properties Cendes et al. (2021) determined for AT2019dsg could equally-well be explained by the CIO model.Stein et al. (2021) and Mohan et al. (2021) suggest that the initially increasing energy and flux density observed from AT2019dsg (as in AT2019azh) is produced by a con-stant energy injection at the source of the outflow, which later switches off and causes the radio emission to fade.At this point the outflow naturally decelerates due to the shutting off of the central engine, rather than due to interactions with the CNM.This could be the case for AT2019azh in the sub-relativistic jet outflow scenario, however, the increasing energy can also be explained by the outflow sweeping up CNM material at a constant velocity (Sironi & Giannios 2013).

SUMMARY
We followed the radio evolution of the tidal disruption event AT2019azh for 850 days post-disruption.The radio emission rose slowly to a peak over 650 d, and is now decaying following the expected decay rate for the Sedov-Taylor solution.We modelled the radio emission as a spherical (or conical), non-relativistic outflow and infer energies of 3.5 × 10 47 -2.7 × 10 50 erg (9 × 10 47 -1 × 10 51 erg), radii of 2 × 10 16 -2 × 10 17 cm (5 × 10 16 -5 × 10 17 cm), and circumnuclear density of 70-1100 cm −3 (15-1000 cm −3 ), for a velocity of ≈ 0.1 c.We detected radio emission approximately 10 d before the optical peak, which is the first radio detection of a thermal TDE at such early times, and well before the X-ray emission indicated any signs of significant accretion.This early time radio detection is in contrast to the suggestion that late-time radio flares or delayed radio emission is common in TDEs (Horesh et al. 2021a).Such an early radio detection could rule out the accretion-driven wind outflow model at early times, but through detailed disk modelling we found that both a sce-nario with delayed accretion or one with early accretion but significant X-ray obscuration could explain the optical, UV, and X-ray properties of the event.
Interestingly, we observed the electron energy index determined from the optically thin part of the synchrotron emission to fluctuate between p ≈ 2.6 to p = 3 − 3.5 from 400 d post-disruption on timescales of months.We deduce that these fluctuations could be due to either inhomogenous emitting regions in a spherical or conical outflow or fluctuations in the energy injection rate of a conical outflow.We rule out the possibility of a relativistic jet producing the outflow in this event since it would require an opening angle of the jet that is smaller than expected for a jet from an SMBH.We found that the mean speed of the outflow to t = 183 d was ∼ 0.1c, thus ruling out a long-lived relativistic outflow.We also found a possible increase in deceleration, suggesting a transition to a Sedov-like evolution after t ∼ 600 d.
We deduce that the unbound debris stream is unlikely to explain the radio properties of this event due to the high energy inferred through an equipartition analysis of the synchrotron emission.We propose the outflow could have originated from an accretion-driven wind or sub-relativistic jet, or a collision induced outflow from stream-stream intersections of the tidal debris.Further detailed modelling of the multiwavelength properties of this event is required to differentiate between outflow models and explain some of the unique properties observed.Future observations of the radio decay of AT2019azh will enable more robust constraints on the CNM density and further insight into the nature of the outflow.

APPENDIX A: CONSISTENCY OF RADIO OBSERVATIONS ACROSS EPOCHS
Due to the use of both the MeerKAT and VLA telescopes, and different configurations at the latter, we performed a consistency check of the data to ensure there were no instrumental or systematic offsets between epochs and between the VLA and MeerKAT L-band observations.We measured the flux densities of 3 background sources in the field of view The data indicate that no systematic offset is present between epochs or between instruments, but that fluctuations on the order of ∼ 10% are present between epochs of VLA observations.for 9 of the VLA epochs and four of the MeerKAT epochs.We applied a primary beam correction to the images before extracting the flux densities.The flux densities for these 3 sources in all bands are shown in Figure A1, demonstrating that there is no significant systematic offset between the two instruments.
In Figure A1 there is evidence for fluctuations between VLA epochs on the order of ∼ 10%.We examined the flux density scale obtained during calibrations of the secondary calibrator for the VLA observations.We found fluctuations of only a few percent between epochs, consistent with the expected flux-density calibration errors.Thus, the fluctuations between epochs are larger than expected due to flux calibration alone, and may be attributed to either intrinsic fluctuations in the sources we examined, primary beam corrections, resolution changes from different configurations, or other calibration inconsistencies.
In the spectral observations of AT2019azh we identify steepenings of the synchrotron energy index.Here we rule out the possibility of these being artificial steepenings due to inconsistencies in the data calibration between epochs.In Figure A1 there is no evidence for inconsistent calibration with frequency for single epochs, except the October 2019 epoch.If the steepenings were artificial, we would expect to see the epochs where the steepenings were observed to show a trend of being lower at the higher frequencies and higher at the lower frequencies than the epochs in which no steepenings were observed.The only epoch where evidence of this frequency trend is present is the October 2019 epoch, and thus we deduce that the slight steepening observed in this epoch (from p ≈ 2.7 to p ≈ 3) is likely not real and merely an artefact of the data calibration.For the other epochs, May 2020 and Feb 2021, there is no evidence in Figure A1 of any such trend with frequency, and we deduce that the spectral steepenings we observed are real.
The flux density offsets between epochs present in Figure A1 would only affect the peak flux densities that we infer from the spectra, and not the spectral slope, since the offsets are not frequency dependent within epochs.Since the peak flux density is often not well constrained in our observations, the uncertainty in the peak flux density is dominated by the spectral fit uncertainty, and the fluctuations of order ∼ 10% between epochs are outweighed.

APPENDIX B: SYNCHROTRON EMISSION FITS INCLUDING A COOLING BREAK
Recently Cendes et al. (2021) detected evidence for a cooling break in radio observations of the synchrotron emission from AT2019dsg.Here we analyse whether a cooling break is detected in our radio observations of AT2019azh.In the case of a cooling break, the data would indicate an additional steepening at a frequency νc > νa, and equation 3 is multiplied by where β3 = -p/2 and s2 is a softening factor (Granot & Sari 2002).
To determine if there is a clear preference for a model with or without a cooling break, we compute the Akaike Information Criterion (AIC) and the Bayesian Information Criterion (BIC) of a spectral fit with a cooling break and one without.In both cases, the best model is indicated by the one with the smallest AIC or BIC.The AIC selects the best predictive model among a number of possibly misspecified models and is given by where M k is the model under consideration, l k is the loglikelihood of the model given its parameters θ k , and p k is the number of parameters estimated by the model M k .
The BIC selects the true model, using a minimal number of parameters and sets a large penalty for models with a larger set of parameters to prevent over-fitting.The BIC is given by where n is the total number of data points that the model is being fit to.
Here we define one model being significantly better than other if the preferred model has an AIC or BIC score of at least 2 units lower.A score of 0-2 means minimal confidence, a score of 2-6 means positive confidence, and a score >6 indicates strong confidence that the model is preferred.
We carried out the same MCMC fitting as described in Section 2 but with the inclusion of the cooling break term in the model for epochs t > 321 d.When the additional cooling break was included in the fits, we found that the uncertainties on the peak frequency, peak flux, and p increased, and the cooling break was identified to be around 4 GHz.In Table B1 we show the calculated AIC and BIC for the original model (model 1) and the model including a cooling break (model 2) for the epochs fit.
In Table B1 it is clear that the AIC and BIC for model 1 is always lower than for model 2, with the model without a cooling break clearly preferred for the epoch at 850 d, and marginally preferred for the other epochs.In this analysis we do not conclusively detect the presence of a cooling break.
This paper has been typeset from a T E X/L A T E X file prepared by the author.

FrequencyFigure 2 .
Figure2.MCMC spectral fits (solid lines) of 12 epochs of radio observations (scatter points) of AT2019azh using the combined VLA and MeerKAT data.random samples from the MCMC chains are plotted after discarding the first 1000 steps for burn-in.Note that the peak flux density and frequency of the first four epochs are not well-constrained due to the lack of low-frequency coverage.

Figure 3 .
Figure3.Spectral fit properties of the radio emission inferred using Equation 3 and shown in Figure2.Upper and lower limits (triangles) are given for the epochs where the peak flux and frequency are not well-constrained by the radio observations and all error-bars represent the 1σ confidence intervals from the MCMC fitting.The peak flux density increases approximately linearly with time, excepting the final epoch, which showed a significant drop in peak flux density.The peak frequency is approximately constant for all well-constrained epochs, with a slight downwards trend after the radio lightcurve peak (t > 600 d).The energy index, p, shows significant fluctuations after 400 d.

)Figure 4 .Figure 5 .
Figure4.Physical properties of the outflow produced in the TDE AT2019azh inferred from an equipartition analysis of the peak radio flux and frequency assuming a spherical, non-relativistic outflow (black) and a conical, non-relativistic outflow (red).Upper limits (triangles) are given for the epochs where the peak flux and frequency are not well-constrained by the radio observations.The energy and radius increase approximately linearly with time until the final epoch.The velocity and magnetic field remain approximately constant over the ∼ 800 d spanned by our observations.

Figure 6 .
Figure 6.The observed values, shown as points with error bars, and the modelled values, shown by solid or dashed lines of the optical/UV (top left) or X-ray (top right) properties of the disk and modelled accreted mass (mass past the ISCO radius, bottom) in the TDE AT2019azh.The ZTF r-and g-band observations are shown in red and green respectively and the Swift UVM2 filter and 2-10 keV X-ray observations are shown in magenta and black respectively.We model two scenarios: a delayed disk-formation (dashed lines) and an early disk-formation (solid lines).The optical and UV observations are well-fit by either model, while the X-ray observations are better-fit by the delayed accretion scenario, unless in the case of significant X-ray obscuration at early times.The X-ray observations are fromHinkle et al. (2021) and only the 0.3-10 keV data from each telescope was used in the fitting.

Figure A1 .
Figure A1.Flux density measurements of 3 background sources for 8 epochs of observations including VLA and MeerKAT data.The data indicate that no systematic offset is present between epochs or between instruments, but that fluctuations on the order of ∼ 10% are present between epochs of VLA observations.

Table 1 .
Radio observations of AT2019azh taken with the VLA and MeerKAT.

Table 2 .
Predicted outflow properties of AT2019azh based on MCMC fitting of the observed radio spectrum and a synchrotron equipartition analysis.Note:All times are reported with reference to t 0 , MJD 58522.The peak of the synchrotron spectrum is not well-constrained by the radio observations. *

Table B1 .
AIC and BIC values for the observed flux spectra for each epoch fit with the original model (model 1) and the model including a cooling break (model 2).