Late-time Hubble Space Telescope Observations of AT 2018cow. II. Evolution of a UV bright Underlying Source 2--4 Yr Post-discovery

In this second of a two-paper series, we present a detailed analysis of three HST observations taken $\sim$2--4 years post-discovery, examining the evolution of a UV-bright underlying source at the precise position of AT 2018cow. While observations at $\sim$2--3 years post-discovery revealed an exceptionally blue ($L_\nu\propto \nu^{1.99}$) underlying source with relatively stable optical brightness, fading in the NUV was observed at year 4, indicating flattening in the spectrum (to $L_\nu\propto \nu^{1.64}$). The resulting spectral energy distributions can be described by an extremely hot but small blackbody, and the fading may be intrinsic (cooling) or extrinsic (increased absorption). Considering possible scenarios and explanations, we disfavor significant contributions from stellar sources and dust formation based on the observed color and brightness. By comparing the expected power and the observed luminosity, we rule out interaction with the known radio-producing circumstellar material as well as magnetar spin down with $B\sim10^{15}\,\mathrm{G}$ as possible power sources, though we cannot rule out the possible existence of a denser CSM component (e.g., previously ejected hydrogen envelope) or a magnetar with $B\lesssim10^{14}\,\mathrm{G}$. Finally, we find that a highly-inclined precessing accretion disk can reasonably explain the color, brightness, and evolution of the underlying source. However, a major uncertainty in this scenario is the mass of the central black hole (BH), as both stellar-mass and intermediate-mass BHs face notable challenges that cannot be explained by our simple disk model, and further observations and theoretical works are needed to fully constrain the nature of this underlying source.

Following the discovery of AT 2018cow, several analogs (sometimes referred to as "Cow-like transients"), CSS161010 , AT 2018lug ("The Koala"; Ho et al. 2020), AT 2020xnd ("The Camel"; Perley et al. 2021;Bright et al. 2022;Ho et al. 2022), and AT 2020mrf (Yao et al. 2022), have also been discovered with defining characteristics being superluminous optical brightness, rapid timescale, and mildly relativistic outflow accompanied by bright multi-wavelength emissions 1 . It has been suggested that these luminous FBOTs likely make up an entirely new population of transients distinct from established SNe and typical FBOTs (Ho et al. 2023). Unfortunately, all analogs were discovered either at much larger distances (z ∼ 0.1−0.3) or after the fact, and thus, the most well-studied case at the moment is still AT 2018cow.
For AT 2018cow, a Tidal Disruption Event (TDE) by an intermediate-mass BH (IMBH, with M BH ≲ 10 5 M ⊙ ; 1 Bright radio emissions were observed for all analogs over the first few hundred days. Bright X-ray emissions were typically observed as well, with the exception of AT 2018lug which did not have any follow-up X-ray observation (i.e., no confirmed X-ray emission). Also, the X-ray detection of CSS161010 was weaker and relatively less luminous. see review by Greene et al. 2020) or a supermassive BH (SMBH, with M BH ≳ 10 5 − 10 6 M ⊙ ) was also proposed initially as a possible explanation (Perley et al. 2019;Kuin et al. 2019). However, the TDE hypothesis gradually became less favored given that it is difficult to prove the existence of such a BH at the outskirt of the galaxy where the gas velocity is smoothly varying without any signs of a coincident massive host system (Lyman et al. 2020), and that a mass limit of M BH < 850 M ⊙ was derived from NICER X-ray quasiperiodic oscillations (QPO; Pasham et al. 2021). The dense CSM around AT 2018cow needed to explain the bright radio emission and non-detection of radio linear polarization would also be difficult to explain unless the BH was already embedded in a gas-rich environment (Margutti et al. 2019;Huang et al. 2019). The exact nature of AT 2018cow is still open to debate and additional constraints are needed to distinguish the viable models. One way of obtaining additional constraints is through late-time observations, useful for probing the immediate surrounding as well as any fading transient emission. For AT 2018cow, the Hubble Space Telescsope (HST) was used to acquire six late-time observations, with the first three tracking the fading prompt emission over ∼50-60 days post-discovery, and the latest three monitoring the field over ∼2, 3, and 4 years post-discovery (in 2020, 2021, and 2022). The HST images taken in 2020 and 2021 were initially examined by Sun et al. (2022), which led to the discovery of an UV-bright underlying source at the precise location of AT 2018cow years postdiscovery. Sun et al. (2022) found this source to be bluer and brighter than any known stars with a stable optical brightness for over a year, which initially led to the suggestion that it may be a young stellar cluster. We independently discovered this underlying source and requested an additional HST observation to monitor any evolution over time. This resulted in the most recent HST observation taken in 2022. The photometry from this 2022 HST observation was briefly reported (but not used) in an environmental study of AT 2018cow by Sun et al. (2023), which showed fading in the NUV, suggesting the existence of a transient component undergoing spectral evolution.
Such transient component that persists over multiple years post-discovery is an interesting discovery that had never previously been observed for an FBOT. In the case of AT 2018cow, compared to the rapid fading of ≳ 8 mag in the NUV over the first ∼60 days, the underlying source faded extremely slowly in the NUV by only ∼0.4-0.5 mag between 703-1453 days post-discovery. This (apparent) drastic difference in timescale may imply a transition to a new evolutionary stage, possibly connected to the commonly-proposed power sources, i.e., ejecta-CSM interaction or central engines. For example, ejecta-CSM interaction is well-known to be capable of producing long-lasting emission over multiple years given a sufficient amount of CSM (e.g., Smith et al. 2009). A remnant central engine can also evolve over long timescales, manifesting at late times as pulsar wind nebulae (e.g., Metzger et al. 2014) or transient accretion disks (e.g., Strubbe & Quataert 2009). Therefore, by connecting the underlying source to a specific remnant power source, significant constraints could be placed on the true nature of AT 2018cow.
In this study, we present a detailed analysis of the UV-bright underlying source at the precise position of AT 2018cow in the three latest HST observations taken ∼2-4 years post-discovery. We examine the highly unusual properties and evolution of this source and use them to place constraints on possible origins and power sources. This is paper II of a two-paper series, while paper I (Chen et al. 2023) focused on the first three HST observations of AT 2018cow that tracked the fading prompt emission (50-60 days post-discovery).
Note that throughout this paper, we adopt the term underlying source to refer to the newly-discovered point source spatially coincident with AT 2018cow, and as with paper I, we use the term prompt emission to refer to the initial evolution of AT 2018cow over the first two months. We opted to use a neutral term for the spatially-coincident object, rather than immediately calling it AT 2018cow, because of the current ambiguity in (i) the exact classification of this object (e.g., partially or entirely transient) and (ii) the exact physical processes producing the emission and if they are related to the initial evolution of AT 2018cow (e.g., entirely new processes or the same processes in a new environment). Therefore, we do not use a term that would suggest a particular classification, physical process, or hypothetical phase of AT 2018cow. This choice of terminology also aligns with our strategy of examining the underlying source as an individual object and associating the hypothetical scenarios with the prompt emission to establish possible links with AT 2018cow.
The HST observations as well as a new late-time Xray constraint from Swift are described in 2. In Section 3, we examine the properties of the underlying source and model the SEDs using simple models. In Section 4, we analyze these properties in the context of different physical scenarios and place constraints on the origin of the underlying source. Finally, we summarize the results and overall implications on the nature of AT 2018cow and the new class of luminous FBOTs in Section 5.
For our analyses, we adopt a luminosity distance of d L = 60 Mpc for AT 2018cow (Perley et al. 2019;Margutti et al. 2019). We assume an R V = 3.1 Milky Way extinction curve with E B−V = 0.07 mag (Schlafly & Finkbeiner 2011) and no internal extinction in the host galaxy of AT 2018cow. Throughout this study, we refer to t as the rest-frame time after the first discovery date of AT 2018cow, MJD 58285.441 Prentice et al. 2018 The latest three HST observations of AT 2018cow monitored the field using WFC3 between ∼2-4 years post-discovery. The first observation (PI Levan) was made during 2020-05-29 (MJD 58998 or t ≃ 703 days), designed to study the local host environment, with images were taken in the F225W, F336W, F555W, F657N, F665N, and F814W bands spanning λ ≃ 2358 − 8029Å. The second observation (PI Filippenko) was made as part of a Snaphot program during 2021-07-25 (MJD 59420 or t ≃ 1119 days) in the F555W (λ = 5308Å) and F814W (λ = 8029Å) bands.
We requested the most recent third observation (PI Chen) after independently discovering an underlying source spatially coincident with AT 2018cow in the first and second observations (subsequently reported in Sun et al. 2022). This observation was made during 2022-06-29 (MJD 59759 or t ≃ 1453) in the F225W, F336W, F555W, and F814W bands. A key goal of this observation was to determine whether the underlying source is transient or stable, which was difficult to discern due to the lack of UV spectral coverage in the 2021 Snapshot program and the higher background level present in the optical bands. A composite image from this observation is shown in Figure 1. All HST data can be found on the Mikulski Archive for Space Telescopes (MAST): 10.17909/fmz6-9b21.
Data reduction and PSF photometry was performed on all three late-time epochs via the same procedure described in paper I (Chen et al. 2023). We identified sources using the drizzled WFC3/UVIS F336W frame from t ≃ 50.3 days as a reference and performed forced photometry at the location of AT 2018cow, allowing slight recentering in each image. We note that there is some astrometric uncertainty resulting from both our centroid on AT 2018cow at t ≃ 50.3 days (≈ 0.002 ′′ based on the signal-to-noise of this detection and the PSF size in F336W) and frame-to-frame alignment (0.005-0.049 ′′ ). Regarding the spatial coincidence between AT 2018cow and the underlying source, we also verified by comparing their relative positions from the nucleus of the host galaxy and found an average offset of ≲ 0.02 ′′ , consistent with the alignment uncertainty.
Final HST photometry of the underlying source are listed in Table 1 and shown in Figure 4. For epochs in which we did not detect any significant (≥ 3σ) emission at the site of AT 2018cow, we injected artificial stars using built-in methods in dolphot. We injected 50,000 stars at magnitudes varying from 20-28 AB mag and estimated the threshold where dolphot could detect 99.7% of all stars at ≥ 3σ significance, which we report as the 3σ limiting magnitude in Table 1.
We found that the underlying source was still present in all bands in our 2022 observations. However, we observed fading of ∼0.4-0.5 mag in the UV bands (F225W Figure 1. Composite HST image of the host galaxy of AT 2018cow from the F225W (blue), F555W (green), and F814W (red) images taken at t ≃ 1453 days (around 4 years post-discovery). A crosshair is shown at the position of AT 2018cow, marking the existence of a spatially-coincident underlying source. The underlying source is much bluer even visually than most other sources in the galaxy. Note that a faint tidal tail can be seen south of the galaxy. and F336W), thus confirming the transient nature of the underlying source (see Section 3). Note-1σ errors are given inside the brackets. The < symbol indicates a 3σ upper limit.
2.1.1. Comparison to Sun et al. (2022Sun et al. ( , 2023 Although our HST photometry generally agree with values presented in Sun et al. (2022Sun et al. ( , 2023 to within about 1σ, there are three discrepancies that we outline in this Section. First, our photometry are almost consistently fainter with larger uncertainties, e.g., our magnitudes in the 2021 Snapshot epoch are m F555W = 25.84± 0.17 AB mag and m F814W = 26.53 ± 0.24 AB mag compared with the previously reported m F555W = 25.60 ± 0.08 AB mag and m F814W = 26.39 ± 0.24 AB mag (note that the latter measurements were originally given in Vega mag in Sun et al. 2022). This discrepancy is likely due to the fact that there is extended background emission close to the site of AT 2018cow. If we instead perform aperture photometry with a radius of 0.12 ′′ , we find that the F555W and F814W measurements are 1 and 0.6 mag brighter than our PSF photometry, respectively. Moreover, while dolphot classifies the F814W source as a point-like object, its crowding is 0.31 mag, which can be understood as the difference between a PSF and aperture magnitude. We conclude that there is some other source of background emission that is compara-  ble in brightness to the underlying source that dolphot deblends from the transient. Differences in dolphot parameters between our work and Sun et al. (2022) may therefore contribute to our fainter measurement. The second discrepancy from Sun et al. (2022) is the F665N observation from 2020, which probes Hα emission at the redshift of AT 2018cow. We report a 3σ upper limit of 23.44 mag, while Sun et al. (2022) reported a 4.6σ detection of 24.47±0.24 mag. While there does ap-pear to be some flux inside a 0.12 ′′ aperture centered at the site of AT 2018cow (Figure 3), it is nominally offset from the transient position, leading to a non-detection in our analysis (which utilized forced photometry based on the transient position as described above). This offset is similar to the frame-to-frame alignment uncertainty, which for F665N is on average 0.028 ′′ or 8 pc at the assumed distance to AT 2018cow. In addition, if we perform aperture photometry with a 0.12 ′′ radius, we only find excess emission at a ≈ 2.8σ level in the F665N band. Therefore, it is not clear if the excess emission in the F665N band is (i) associated with the underlying source and (ii) significantly above the diffuse background. We discuss the potential excess narrow band emission further in our interpretations in Section 4.
The third discrepancy is the F225W measurement at the latest HST epoch (t ≃ 1453 days), for which we report 24.93±0.11 AB mag while Sun et al. (2023) reported a fainter 25.21 ± 0.07 AB mag (again from Vega mag). This difference is unexpected, especially since the background level is the lowest in this band, and our other measurements at this epoch are generally consistent. We also found that aperture photometry leads to an F225W  HST photometry are shown as stars, and earlier measurements at similar wavelengths taken from paper I (Chen et al. 2023) are also shown. Offsets are added for visual clarity. Dash-dotted lines corresponding to power law declines of flux density are also added for reference (note these are not best-fit lines). Fading is observed in the F225W and F336W bands for the underlying source, but at a much slower rate than the initial decline of AT 2018cow. measurement 0.37 mag brighter than the PSF photometry, meaning that the issue is unlikely related to excess background. The exact cause of this discrepancy is unclear, and we chose to use our PSF photometry for our analyses. Note that this discrepancy does not significantly impact the results presented in this study because the fading would be even more prominent according to the measurement from Sun et al. (2023).

Late-time Swift XRT Observation (t ≃ 1360 days)
Finally, to constrain any high energy emission associated with the underlying source, we analyzed two X-ray observations (totaling 8.5 ks) obtained by the Swift X-Ray Telescope (XRT; Burrows et al. 2005) on 2022-03-25 and 03-27 (or t ≃ 1360 days). This corresponds to approximately 93 days prior to the latest HST observations taken on 2022-06-29. We downloaded the cleaned event files and the exposure images from the Swift archive and extracted the X-ray images using xselect. We then followed the mosaic routine in ximage to combine the images from the two observations.
In the combined image, no X-ray emission was detected at the site of AT 2018cow. Using the sosta function of ximage, we inferred a 3σ upper limit count rate of 0.00157 cts/s at the location of AT 2018cow. Assuming a power law with a photon index of Γ = 2 and a galactic neutral hydrogen column density of N H = 5 × 10 20 cm −2 (Margutti et al. 2019), the count rate corresponds to an unabsorbed flux limit of F 0.3−10keV < 6.27 × 10 −14 erg s −1 cm −2 (or L 0.3−10keV < 2.70×10 40 erg s −1 at the distance of AT 2018cow). However, we emphasize that since the origin of the underlying source is unknown, there is significant uncertainty in the X-ray spectral shape. Varying the spectral shape can change this flux limit by more than an order of magnitude. This X-ray non-detection is further assessed in the context of various transient models in Section 4.

EVOLUTION OF THE UNDERLYING SOURCE
Cutout HST images of the underlying source are shown in Figure 2 (wide bands) and Figure 3 (narrow bands). The underlying source is detected as a point Figure 5. The HST SEDs of the underlying source (dereddened from Galactic extinction) and best-fit power laws (dotted lines) and blackbodies (dashed lines), as well as the 3σ Swift-XRT upper limits (downward triangles) assuming a power law and fiducial photon indices Γ. Shaded regions represent 1σ uncertainties for the best-fit blackbodies. In the zoomed-in plot, the HST SED of AT 2018cow at t ≃ 60.1 days with brightness lowered by a factor of 20 is also shown for comparison. Although the optical emission seems relatively stable, the fading UV emission confirms the transient nature of the underlying source. source in all three epochs, most clearly in the UV bands (F225W and F336W). Compared to its immediate surrounding, the underlying source is the only bright point source in the UV (also see Figure 1) but is hardly discernable visually from the local diffuse emission in the red optical bands. The light curves of the underlying source are shown in Figure 4, revealing fading in the NUV. In this section, we describe the basic properties and evolution of the underlying source. Figure 5 shows the SEDs of the underlying source derived from the three HST observations as well as the Swift-XRT upper limits at t ≃ 1360 days (plotted assuming a power-law spectrum with various photon indices Γ for illustrative purposes). The underlying source is detected well above the background in the F225W, F336W, and F555W bands, with signal-to-noise ratios of SNR ≃ 9 − 14. The detection is weaker in the F814W band, with SNR = 5.4 and 3.3 at t ≃ 703 days and 1453 days, respectively. In the F657N and F665N narrow bands (t ≃ 703 days), we did not detect emissions above 3σ (see discussion about F665N band in Section 2.1.1). Here we highlight several important properties of the underlying emission.

Basic Properties of the UV-Optical SEDs
SED Shape: The exceptionally blue color was the most remarkable property of the underlying source, quite distinct from most other sources in the galaxy (visually recognizable in Figure 1). The SED peak was not directly constrained in any of the HST epochs, suggesting λ SED,peak ≲ 2358Å, while the lack of an (0.3 − 10 keV) X-ray detection shortly before the latest HST epoch implies that the peak must be further in the UV at t ≳ 1360 days. At t ≃ 703 days and 1453 days, we found a dereddened color of F336W − F555W = −1.3 mag and −1.0 mag, respectively, which are much bluer than even AT 2018cow at peak. Finally, there is tentative evidence for a change in the SED slope between the UV and optical bands, with a shallower slope between F225W and F336W and a steeper slope between F336W and F555W relative to a smooth power law (see Figure 5).
Brightness: The NUV brightness of the underlying source is also surprisingly bright, with an absolute (dereddened) magnitude of -10.1 AB mag in the F225W filter at t ≃ 703 days. This is only about 20 times fainter than AT 2018cow at t ≃ 60 days, despite being observed almost two years later ( Figure 5). In contrast, we found that the NUV brightness of AT 2018cow decreased by more than 6000 times over the initial 60 days. In terms of a power law decline, we found F NUV ∝ t k with k ∼ 0.7 − 1.2 over t ≃ 60 − 1453 days, much slower compared to the rates of t −2 − t −3 over the first 60 days (see Figure 4). We note that the power law decline of flux density in the NUV does not directly correspond to the decline of bolometric luminosity due to the significant change in color/spectrum. Without assuming the spectral shape, we can calculate the minimum UVoptical luminosity L UVO,min of the underlying source by directly integrating the SEDs. The values of L UVO,min at t ≃ 703 days and 1453 days are given in Table 2, which are on the order of 10 39 erg s −1 .

Transient Nature of the Underlying Source
Previous studies (e.g., Sun et al. 2022;Metzger 2022) suggested that the relative stability of the optical emission of the underlying source between t ≃ 703 days and 1119 days could be an indication of a stable stellar source. However, from the most recent HST observation, we found significant fading in the UV bands between t ≃ 703 days and 1453 days. Specifically, we found the fading to be 0.54 ± 0.15 mag (3.5σ) and 0.43 ± 0.11 mag (3.9σ) in the F225W band and F336W band, respectively. In contrast, we found no significant fading (only ∼1σ) in the optical bands, indicating ongoing spectral evolution of the underlying source.
Note that we have performed additional checks to verify the robustness of the observed fading in the UV. First, we checked the photometry of stable sources in the HST images and found consistent brightnesses between different epochs, meaning that the measured fading of the underlying source was not due to any systematic calibration issue. Second, we performed aperture photometry on the underlying source (using apertures show in Figure 2) and recovered similar fading in the UV bands, confirming that the fading was not caused by, e.g., dolphot deblending parameters. Together, these confirm the transient nature of the underlying source, which could be intrinsic (i.e., remnant transient of AT 2018cow) and/or extrinsic (i.e., increased absorption along the line of sight).

Modeling of the UV-Optical SEDs
We performed a set of simple fits to the HST SEDs to further constrain the properties of the underlying source, which we use as a basis for our discussion in Section 4 when considering specific physical scenarios. Specifically, we characterized (i) the SEDs through power law and blackbody models and (ii) the observed fading through an extinction law.

Power Law and Blackbody Models
We fit two models to the HST SEDs: a power law in the form L ν ∝ ν α and a blackbody. We performed forward modeling using the Markov Chain Monte Carlo (MCMC) sampler in the Python package emcee (Foreman-Mackey et al. 2013). The best-fit parameters and uncertainties were derived from the 50th, 15.9th, and 84.1th percentile of the resulting samples. Note that we did not fit the SED at t ≃ 1119 days as photometry was only available in two bands.
The best-fit spectral index α, blackbody temperature T BB , blackbody radius R BB , as well as the blackbody luminosity L BB are given in Table 2. The resulting fits are plotted in Figure 5. It is worth noting that the SED at t ≃ 703 days was so blue that the spectral index was α ≃ 2, the expected value for the Rayleigh-Jeans tail (i.e., L ν ∝ ν 2 ). At this epoch, if the emission was blackbody, the blue color would imply a very high temperature (T BB > 10 5 K), but the brightness would suggest an incredibly small size (R BB ≲ 15 R ⊙ ). The flattening in the SED at t ≃ 1453 days can be seen in the decreasing α, which could suggest cooling (lower T BB ) and expansion (larger R BB ) of a blackbody.
These properties of the underlying source are all quite extreme in the context of late-time observations of transients. In particular, if the blackbody characterization is accurate, the derived temperatures are even higher than those of AT 2018cow, while the radii are no different from the sizes of stars. The derived blackbody luminosities are on the order of L BB ∼ 10 40 − 10 43 erg s −1 , similar to those of AT 2018cow at the second month postdiscovery. Note that because the temperatures are so high and the peak wavelengths are further in the UV, the observed luminosities L UVO,min are only ∼0.03-2% of L BB , meaning that only a tiny fraction of the radiation was actually observed in the NUV-optical bands. Overall, the high temperature and luminosity could be an indication of additional energy injection from a remnant power source of AT 2018cow. In this case, constraining the nature of the underlying source and the associated power source could be crucial in revealing the exact identity of AT 2018cow. We further examine possible power sources through additional modeling and discuss the implications in Section 4.

Extinction Law
We also considered a different scenario where the transient phenomenon was not in the radiation but rather in the extinction that led to preferential dimming at shorter wavelengths. In this case, newly-formed dust grains over t ≃ 703 − 1453 days increased the extinction along the line of sight and caused the apparent fading in the UV. To test this case, we assumed the source spectrum to be the best-fit power law at t ≃ 703 days and reddened the spectrum according to the Cardelli extinc-  tion law (Cardelli et al. 1989) with R V = 3.1 to fit the observed SED at t ≃ 1453 days through synthetic photometry. Figure 6 shows the resulting reddened power law with a best-fit color excess of E B−V ≃ 0.072, which matches the observed SED at t ≃ 1453 days. Therefore, dust extinction could be a possible explanation for the observed fading of the underlying source, especially given that dust has been proposed to explain the excess IR of AT 2018cow (Metzger & Perley 2023). We further discuss the implications of dust formation in Section 4.2.

CONSTRAINTS ON THE ORIGIN OF THE UNDERLYING SOURCE (t ∼ 703 − 1453 days)
Here, we briefly summarize the properties of the underlying source outlined in Section 3 before discussing implications for its origin. The underlying source was quite bright, with an integrated (minimum) luminosity of L UVO,min ∼ 10 39 erg s −1 over λ ∼ 2358 − 8029Å. HST photometry of the underlying source showed an exceptionally blue continuum (F336W − F555W = −1.3) without constraining the peak of the SED (λ peak ≲ 2358Å). A Swift-XRT non-detection at t ≃ 1360 days suggests that the peak was in the UV during the latest HST observation ( Figure 5). The NUV-optical continuum at t ≃ 703 days matches L ν ∝ ν 2 , the Rayleigh-Jeans tail, and can be described by a blackbody with a high temperature T BB ≳ 10 5 K and a small radius R BB ≲ 20 R ⊙ . Transient nature was also confirmed by the 3.5σ and 3.9σ fading in the two NUV bands (F225W and F336W) between t ≃ 703 days and 1453 days, which flattened the spectrum and could indicate cooling and expansion of a blackbody (see Table 2).
The transition between rapid fading over the first two months to slow fading over a timescale of years is likely associated with a sustained power source. However, the exceptionally blue color immediately rules out synchrotron emission and warm dust emission with T dust < 2000 K while the brightness also rules out radioactive decay. If we assume that L UVO,min ∼ 10 39 erg s −1 at t ≃ 1453 days was powered by radioactive decay, following the equations in Afsariardchi et al. (2021) and assuming full γ-ray trapping, we find a completely unphysical M Ni ∼ 74 M ⊙ .
In this section, focusing on the HST SEDs and the Swift-XRT non-detection, we discuss the constraints on the possible power sources of the underlying source. We consider five possible origins of the underlying emission and the observed fading: significant stellar contribution (Section 4.1), dust extinction (Section 4.2), ejecta-CSM interaction (Section 4.3), a magnetar (Section 4.4), and an accreting BH (Section 4.5). Note that we assume the fading to be smooth over t ≃ 703 − 1453 days and do not consider any sporadic activities (e.g., multiple flares) that could explain the observations because we do not have sufficient temporal coverage to distinguish such a possibility.

Significant Stellar Contribution
Although the fading of the underlying source over ∼2-4 years post-discovery is likely associated with some transient phenomenon, there is a possibility that an underlying stellar source contributed to the emission. In particular, the slow fading in the UV and relative stability in the optical might be easier to explain if stellar emission contributed significantly to the HST SED at t ≃ 1453 days. In this case, the actual transient emission would have faded faster, but could not be observed once it became dimmer than the stellar emission. To understand the implications of this scenario, we consider the limiting case where the entire HST SED observed at t ≃ 1453 days was stellar emission. We consider implications for both the stellar and transient emission.

Implications for the Stellar Population
First, we note that any stellar emission is unlikely associated with an isolated "usual" massive star (M ≲ 100 M ⊙ ). This was disfavored by Sun et al. (2022) based on brightness and color. Hypothetically, Very Massive Stars (VMS) with M ≳ 100 − 200 M ⊙ , L ≳ 10 39 erg s −1 , and T eff ≳ 10 4.5 K (e.g., Sabhahit et al. 2022) could perhaps explain the color and brightness of the underlying source. However, this would require a surviving single VMS star to be coincident within ≲6 pc of the location of AT 2018cow, which is unlikely unless AT 2018cow also came from a VMS progenitor. The relative isolation of two VMS at the outskirt of the galaxy would be unusual, and although failed explosion is a proposed scenario for FBOTs (e.g., Kashiyama & Quataert 2015b), a VMS progenitor for AT 2018cow is still less favored for various reasons, such as the small ejecta mass inferred from the prompt transient emission (Margutti et al. 2019) and the possibility that AT 2018cow comes from an older stellar population that resides in the foreground relative to the nearby star-forming regions (Sun et al. 2023).
Thus, a more natural scenario may be that any stellar contribution comes from an underlying star cluster, possibly the host of the progenitor of AT 2018cow. To investigate the properties of this hypothetical cluster, we compared the HST SED at t ≃ 1453 days with stel The left panel of Figure 7 shows the HST SED at t ≃ 1453 days and some selected clusters from BPASS and LEGUS. Compared to the observed SED, the youngest BPASS clusters (1 Myr old) have the best-matching color but overpredict the luminosity by several orders of magnitude. The normalization of the BPASS clusters is based on a cluster with M clus = 10 6 M ⊙ , meaning that a lower cluster mass (M clus ∼ 10 3 − 10 4 M ⊙ ) is required to explain the brightness of the underlying source.
On the other hand, we were able to identify a number of LEGUS clusters with brightness and color similar to the underlying source. However, these are extremely rare cases (we found less than a dozen searching through the cluster catalogs), highlighting the peculiarity of the brightness and color. In the images, these clusters are seen in blue crowded regions of the galaxy that appear to be associated with young stellar populations and star-formation activity 2 . In the catalogs, they are classified as asymmetric and multi-peaked sources (Class = 2 and 3) with a best-fit age ∼few Myr and best-fit M clus ∼ few × 10 3 M ⊙ , consistent with our expectation from the BPASS clusters (however, note that based on the given Q probability, the SED fits are quite poor). Overall, these comparisons suggest that it is possible to explain the HST SED at t ≃ 1453 days as a star cluster (with real examples). However, it would imply a rare young (age ∼ Myr) cluster with a relatively small mass (M clus ∼ 10 3 M ⊙ ) and suggest that AT 2018cow came from a massive progenitor.
We briefly note that we have not considered cases involving long-term evolution of rare stellar merger/interaction events that leave behind compact or stripped stars with some energy source powering a low-mass envelope. An example is the case considered by Cohen & Soker (2023) where, after the common envelope jets supernova imposter event that produced AT 2018cow (Soker 2022), a NS is inside the red supergiant during a second common envelope phase and influencing the evolution through launching jets. Although these objects could in theory be blue and luminous and become redder over time, their exact spectral evolution is not clear at the moment, and it is uncertain if they can fully explain the underlying source.

Implications for the Transient Emission
Next, we consider the implications on the transient emission if the entire HST SED at t ≃ 1453 days was from a star cluster. In this case, we can subtract the cluster SED at t ≃ 1453 days from the observed SED at t ≃ 703 days to obtain a transient SED at t ≃ 703 days. The right panel of Figure 7 shows the transient SED and the best-fit power law. The transient emission in this case is still bright, with a L UVO,min = (5.3 ± 1.0) × 10 38 erg s −1 . The transient SED is much bluer now with a best-fit spectral index of α = 2.5. This extremely hard transient spectrum would be inconsistent with a blackbody and may be significantly challenging to explain with other emission mechanisms.
In addition, if we consider a scenario where only the optical emission at t ≃ 1453 days is dominated by a stellar population and the UV is dominated by the transient (e.g., an older stellar population), then the inferred spectral index of the transient emission at t ≃ 703 days would be even steeper. Thus, while a star cluster could have contributed significantly to the HST SED at t ≃ 1453 days, the implied blue color (spectral index α ≳ 2.5) of the transient emission at may be difficult to explain, posing a significant challenge to this scenario.

Dust Extinction
As we have shown in Section 3.3 and Figure 6, an increase in extinction could theoretically account for the preferential fading in the UV of the underlying source. Assuming the Cardelli extinction law (Cardelli et al. 1989) with R V = 3.1, we found that the fading between t ≃ 703 days and 1453 days would correspond to a color excess of E B−V ≃ 0.072. The increase in extinction would imply dust formation over these epochs, likely associated with AT 2018cow and its surrounding CSM. Here, we examine the properties of a hypothetical dust cloud that could produce the observed fading and discuss the implications of dust formation on the nature of the underlying source.
We can make an order-of-magnitude estimate of the required dust mass column density Σ d from the inferred extinction. Assuming spherical dust grains with radius a and density ρ and uniform size and composition, the optical depth through the dust cloud τ λ (a) can be written as where Q ext λ is the extinction coefficient. The extinction coefficient is given by the sum of the absorption coefficient and scattering coefficient: Q ext λ = Q abs λ + Q sca λ . The optical depth is related to the extinction in magnitude through A λ = 2.5 log 10 (e)τ λ while the visual extinction in magnitudes is calculated by Therefore, we can derive a dust mass column density as: We considered two types of dust, graphite with ρ = 2.26 g cm −3 or silicate with ρ = 3.30 g cm −3 , and two grain sizes, a = 0.1 µ m and 1.0 µ m. We use Q abs λ and Q sca λ derived in Draine & Lee (1984) and Laor & Draine (1993) 3 and interpolated the coefficients to λ V = 5500Å to obtain Q ext V . From these values and assuming R V = 3.1 and E B−V = 0.072, we estimated Σ d (a = 0.1 − 1.0 µ m) ∼ 10 −5.7 − 10 −4.6 g cm −2 ∼ 10 −4.9 − 10 −4.4 g cm −2 for graphite (top) and silicate (bottom).
To convert the column density to a volume density or a mass, we have to consider the size of the dust cloud. For the underlying source, dust formation over t ≃ 703 − 1453 days likely occurred in the region where pre-existing dust grains were initially destroyed by AT 2018cow, but the material cooled over time and started forming new dust grains at these late epochs. The UV radiation of AT 2018cow would have destroyed the dust grains inside the sublimation radius, where the temperature is higher than the sublimation temperature of the dust grains (∼1100-1500 K for silicate and ∼2000 K for graphite). If we assume a luminosity of L ∼ 4 × 10 44 erg s −1 (peak of AT 2018cow), a grain size of a ∼ 0.1 − 1.0 µm, and a sublimation temperature of T s ∼ 1100 − 2000 K, then we can follow the equations in Metzger & Perley (2023) and obtain a sublimation radius r s ∼ 10 16.2 − 10 16.4 cm. The ejecta of AT 2018cow also could have destroyed the dust, which would have reached r ∼ 10 17.6 − 10 17.9 cm by t = 1453 days assuming a velocity of v ∼ 0.1c − 0.2c. Therefore, we can assume that the newly-formed dust cloud has a radius on the order of R d ∼ 10 16.2 − 10 17.9 cm.
Note that if the underlying source was at the center of the dust cloud (i.e., associated with AT 2018cow), its UV emission at t ≃ 1453 days (L ∼ 10 40 erg s −1 ) would create a dust-free cavity with a radius of r s ∼ 10 15.8 − 10 16.0 cm. This size is slightly less significant than the order-of-magnitude size we consider here, so we ignore this cavity. Detailed models should consider dustforming and dust-destroying regions inside an evolving radiation field, which is beyond the scope of this study.
From the column density and physical size, we estimated the dust mass from M d ∼ πR 2 d Σ d , which gave assuming Σ d ∼ 10 −5.7 − 10 −4.4 g cm −2 . Note that as shown in paper I (Chen et al. 2023), the dust mass derived from the excess IR emission of AT 2018cow was on the order of M d ∼ 10 −6 − 10 −4 M ⊙ , suggesting that the dust mass needed to explain the fading of the underlying source is not unreasonable. We can also infer a gas density by assuming a dust-togas mass density ratio ρ d /ρ g = X d ∼ 0.1. We estimated the gas density from ρ g ∼ Σ d /R d X d , which gave ρ g ∼ 10 −20.9 − 10 −19.6 g cm −3 (R d ∼ 10 16.2 cm) ∼ 10 −22.6 − 10 −21.3 g cm −3 (R d ∼ 10 17.9 cm) assuming Σ d ∼ 10 −5.7 − 10 −4.4 g cm −2 . Interestingly, this gas density is roughly consistent with the radioproducing CSM density at R ∼ 10 16 − 10 17 cm (see top panel of Figure 8). However, we note that the exact value of X d is very uncertain and can depend on properties such as the composition. In particular, if X d is orders of magnitude lower than 0.1 (e.g., hydrogen-rich gas), the inferred gas could be much denser than the radio-producing CSM. In that case, a possible explanation is that the dust was associated with the dense hydrogen-rich envelope ejected by the progenitor of AT 2018cow located at R ≳ 10 17.3 cm, which is outside the range covered by published radio observations.
Overall, the properties of the hypothetical dust cloud appear reasonable in the context of AT 2018cow and stellar explosions. However, there are challenges in explaining the actual underlying emission in this scenario. If the fading in UV over t ≃ 703 days − 1453 days was purely due to dust extinction (with E B−V ≃ 0.072), then the underlying emission may have been stable with a spectrum similar to or likely bluer than the observed HST SED at t ≃ 703 days (i.e., spectral index α ≳ 2). Although the stability would suggest a stellar source, if we follow the findings in Section 4.1, there are no known stellar sources that can explain such a blue color. Thus, in this scenario, we are left with the possibility of a remnant transient emission that is stable over several years with an extremely blue spectrum (α ≳ 2), which is still fairly challenging to explain. Therefore, if dust was relevant to the underlying source, it is more likely that dust extinction only contributed partially to the fading. This scenario could be very complex, possibly involving some combination of dust, stellar source, and transient source. For this case, the dust and gas densities derived above assuming E B−V ≃ 0.072 should be considered as upper limits. Constraining such a scenario would require additional observations monitoring the evolution of the underlying source and detailed modeling. Deep images through the JWST may also be useful in constraining potential IR echo produced by the hypothetical dust cloud.

Ejecta-CSM Interaction
Ejecta-CSM interaction is known to be capable of producing prolonged emission in interacting SNe such as SN 2005ip (Smith et al. 2009) and SN2010jl (Fransson et al. 2014Jencson et al. 2016), with optical luminosities of L ≳ 10 41 erg s −1 over hundreds of days. Although no known interacting SN matches the description of both AT 2018cow and the underlying source, the existence of CSM around AT 2018cow may nonetheless contribute significantly to the underlying emission. Here, we model the interaction between the fast ejecta (v ∼ 0.1c) of AT 2018cow and the known radio-producing CSM to check if this interaction years after the explosion can power the underlying emission.
Our ejecta-CSM interaction model is similar to the method described in Section 5.2 of Chandra et al. (2015) for arbitrary ejecta and CSM density distributions. Under the thin shell approximation (Chevalier 1982), we solve the equation of motion to advance the shock velocity v s using the Runge-Kutta method. We advance the shock radius R s and calculate the ejecta and CSM densities, ρ ej and ρ CSM . At the new time step, for the forward shock (fs) and reverse shock (rs), we calculate the velocities v fs = v s and v rs = R s /t − v s , temperatures T fs and T rs assuming a strong shock and ion-electron equipartition, and cooling times t c,fs and t c,rs (following Nymark et al. 2006). Finally, we calculate the kinetic luminosities L kin,fs = 2πR 2 s ρ CSM v 3 fs and L kin,rs = 2πR 2 s ρ ej v 3 rs and derive radiated luminosities in the form of L rad = ηL kin , where we define a radiation efficiency factor η = t/(t + t c ).
We assumed power law distributions for the ejecta: ρ ∝ r −δ with δ = 0 for the inner region and ρ ∝ r −n with n = 12 for the outer region (Chevalier & Fransson 1994). From the radio observations of AT 2018cow (i.e., Ho et al. 2019;Nayana & Chandra 2021), we adopted a two-component CSM distribution withṀ = 4 × 10 −4 M ⊙ yr −1 inside a radius of R ≤ 1.7 × 10 16 cm andṀ = 4 × 10 −6 M ⊙ yr −1 outside a radius of R ≥ 6 × 10 16 cm, assuming ρ =Ṁ /(4πR 2 v w ) with v w = 1000 km s −1 for both distributions. A steeper power law was used to connect the two distributions between R = 1.7 × 10 16 cm and 6 × 10 16 cm. For the explosion, we assumed an ejecta mass of M ej = 0.5 M ⊙ and a total energy of E = 3 × 10 51 erg, motivated by the short rise time and high peak optical luminosity of AT 2018cow (Margutti et al. 2019). We also assumed solar composition with a mean molecular weight of µ = 0.61. Figure 8 shows the model CSM density distribution, shock temperatures, cooling times, and kinetic and radiated luminosities. Dashed lines are shown to mark the three HST epochs. Properties for the forward shock and reverse shock are shown in blue and red, respectively. The cooling times indicate that both the forward and reverse shocks are adiabatic (i.e., t cool ≫ t) from early times because of the high velocity and low ejecta mass. The radiation efficiency η is thus very low, and the radiated luminosity is orders of magnitude lower than the kinetic luminosity.
From our model, the radiated forward shock luminosity and reverse shock luminosity are L rad,fs < 10 33 erg s −1 and L rad,rs < 10 36 erg s −1 , respectively, which are orders of magnitude lower than the minimum UV-optical luminosity of the underlying source, L UVO,min ∼ 10 39 erg s −1 (bottom panel of Figure 8). The kinetic luminosity of the forward shock at L kin,fs ∼ 5 × 10 40 erg s −1 could explain the luminosity of the underlying source if the radiation efficiency is high enough (η ≳ 0.02). However, given the adiabatic nature of the shock, most of the radiation should be in the form of free-free emission in the X-ray and would have difficulty accounting for the observed UV-optical emission. Such a high radiation efficiency would also imply that the interaction should have dominated the fading prompt emission of AT 2018cow, which was not observed.
Therefore, we conclude that the fast ejecta of AT 2018cow interacting with the known radio-producing CSM could not account for the observed UV-optical emission of the underlying source. Note that our model does not include the proposed dense equatorial CSM with slow interaction for AT 2018cow (Margutti et al. 2019). Including such a component may result in an appreciable change but is outside the scope of this study because it would first require a detailed aspherical ejecta-CSM models for AT 2018cow to constrain the correct CSM profile.
We do not rule out the possibility of ejecta-CSM interaction involving a new CSM component, i.e., the previously ejected hydrogen envelope of the progenitor of AT 2018cow. Such a phenomenon has been previously observed in cases such as SN 2014C, where a hydrogenpoor SN transitioned to an interacting hydrogen-rich SN (e.g., Milisavljevic et al. 2015;Margutti et al. 2017;Mauerhan et al. 2018;Brethauer et al. 2022). However, from the interaction with hydrogen-rich CSM, an expected signature is strong Hα emission (e.g., Fransson et al. 2014). For the underlying source, the detection of Hα emission is uncertain (see Section 2.1.1). Even assuming an Hα luminosity of L Hα ∼ 4 × 10 36 erg s −1 (from Sun et al. 2022), this would still be many orders of magnitude fainter than those typically observed in  . CSM density profile used in our ejecta-CSM interaction model (top panels) and the forward and reverse shock properties derived from the model (bottom three panels). See text for details regarding our model. The three HST epochs are marked with gray dashed lines. We also show LUVO,min (triangles), and a key finding here is that LUVO,min is orders of magnitude above the L rad . interacting SNe, e.g., L Hα ∼ 10 39 erg s −1 for SN 2014C (Mauerhan et al. 2018). Detailed modeling as well as additional observations could help constrain this hypothesis. In particular, in addition to deeper Hα images, radio and X-ray follow-ups may also provide significant constraints on any possible rebrightening from the new interaction, and deep IR observations could also probe warm dust associated with the CSM.

Magnetar Spin Down
Studies have suggested that a millisecond magnetar could power the fast-rising luminous optical peak of AT 2018cow (e.g., Prentice et al. 2018;Margutti et al. 2019;Ho et al. 2019;Mohan et al. 2020;Xiang et al. 2021). Following this hypothesis, the underlying source could be remnant emission powered by such a magnetar.
In the magnetar scenario, a fraction of the spin-down luminosity L sd is radiated as the observed underlying source, depending on the efficiency. L sd depends on the rotational energy E sd and spin-down timescale t sd , which are related to the surface magnetic field strength B and the initial spin period P 0 of the magnetar. Therefore, we can place constraints on the magnetar configuration (i.e., B and P 0 ) that would be required to explain the underlying source by comparing L sd with the observed L UVO,min .
We calculated E sd , t sd , and L sd as a function of B and P 0 following the formulas in Section 4.3 of Kasen (2017). We assumed a neutron star with a mass of M ns = 1.4 M ⊙ and a radius of R ns = 10 km and assumed a moment of inertia in the form of I ns = 2M ns R 2 ns /5. Note that for timescales relevant for the underlying source (t ≃ 703 − −1453 days), t ≫ t sd for large B and L sd ≈ E sd t sd /t 2 ∝ B −2 , i.e., L sd is actually independent of P 0 . In Figure  9, we show the magnetar parameter space (B vs. P 0 ) with lines of constant luminosity for L sd = L UVO,min , 10L UVO,min , 100L UVO,min .
We found that L sd = L UVO,min for B ≃ 10 15 G, meaning that if the magnetar has an extreme field strength, by these late times, all of its spin-down luminosity is required to power just the observed UV-optical luminosity over λ ∼ 2300 − 8000Å. However, since the spectral peak of the underlying source could be further in the UV (λ peak ≲ 2300Å), the observed UV-optical luminosity is likely only a small fraction of the total luminosity. Therefore, a magnetar with B ≳ 10 15 G cannot power the emission of the underlying source.
Since L sd ∝ B −2 at these late times, a larger L sd can be obtained from a smaller B. We found that at B ≃ 10 14 G, L sd ≃ 100L UVO,min , the spin-down luminosity is more than two orders of magnitude larger than the observed luminosity. Therefore, a magnetar with B ≲ 10 14 G could hypothetical have enough energy output to power the emission of the underlying source.
We also show a shaded purple region in Figure 9 to indicate the parameters B and P 0 required to produce the energy and timescale associated with the optical peak of AT 2018cow: E sd ∼ 10 50.5 − 10 51.5 erg and t sd ∼ 10 3 − 10 5 s (ranges from Margutti et al. 2019, also see Prentice et al. 2018 andXiang et al. 2021 regarding magnetar parameters). Although a magnetar with B ∼ 10 15 G could power the optical peak of AT 2018cow, we argued above that it does not have sufficient energy output to power the underlying source. On the other hand, a magnetar with a lower field strength B ∼ 2 − 4 × 10 14 G may be able to power both emissions if the observed underlying emission L UVO,min is a at least few percent of the total luminosity. However, the observed underlying emission could be much less than 1% of the total luminosity, which would be the case if the emission at t ≃ 703 days was blackbody (see Section 3.3). In that case, a NS could either power the peak of the AT 2018cow (with a magnetar) or the underlying emission (with a less magnetized pulsar), but not both. A third possibility is that AT 2018cow did not involve a NS at all. Detailed magnetar models and additional observations of the underlying source are required to further constrain these possibilities.

Black Hole Accretion Disk
Many models of AT 2018cow involve various kinds of winds and jets that could be driven by newborn central engines (e.g., Piro & Lu 2020;Gottlieb et al. 2022;Metzger 2022). Some studies have also suggested that AT 2018cow could be a TDE involving an IMBH or a SMBH (Kuin et al. 2019;Perley et al. 2019). Following these scenarios, a remnant accretion disk is a natural hy- HST (t = 703d) HST (t = 1119d) HST (t = 1453d) Figure 10. Disk blackbody fits to the latest HST SED of the underlying source at t ≃ 1453 days. All three (dereddened) HST SEDs are also shown for reference. The order-ofmagnitude estimates of the upper limit Lν derived from the Swift-XRT non-detection at t ≃ 1360 days are also shown. The horizontal span of the limit represents the range of 0.3 − 10 keV, while the vertical span represents the range of Lν derived from the range of MBH.
pothesis that could explain the slow-evolving, extremely hot and small underlying source. Here, we examine the remnant accretion disk scenario and place constraints on the disk configuration and the central mass by modeling the HST SEDs and the Swift-XRT upper limit. Note that in the following sections, although our simple disk model for the underlying source is agnostic to the exact nature of the central object, we refer to the it as a BH. While the central object in principle could be a NS under certain configurations, we prefer a BH in the context of our analysis because (1) much of the mass range we examine (up to 10 6 M ⊙ ) are BHs and (2) for masses that can be NSs, the extreme super-Eddington accretion required in our model (up to about 1 M ⊙ yr −1 ; see Section 4.5.3 and 4.5.4) over multiple years would likely provide enough mass to exceed the NS mass limit.

Multi-Temperature Disk Blackbody Fit
We fit the latest HST SED of the underlying source (t ≃ 1453 days) using a multi-temperature disk blackbody model (Mitsuda et al. 1984;Makishima et al. 1986), which describes a standard geometrically thin accretion disk (Shakura & Sunyaev 1973;Pringle 1981). This model assumes a disk with near-Keplerian orbits that is radiatively efficient (cool disk with negligible radial pressure gradient and supersonic orbits) and geometrically thin (small aspect ratio). The disk is also assumed to be optically thick, with each annulus having a local temperature and radiating as a blackbody.
We use this simple model as an initial approach to broadly explore the parameter space for the accretion scenario in the context of the underlying source. We note that while the assumptions may be reasonable for sub-Eddington accretion, the thin disk approximation does break down for super-Eddington accretion where the disk is radiatively inefficient and geometrically thick (e.g., Paczyńsky & Wiita 1980;Abramowicz et al. 1988; also see reviews by Frank et al. 2002;Abramowicz & Fragile 2013). This caveat becomes relevant for lower BH masses when modeling the underlying source (see Section 4.5.3 and 4.5.4).
We assumed a temperature profile as a function of radius of T (r) = (3GṀ M BH /8πσr 3 ) 1/4 , where G is the gravitational constant, σ is the Stefan-Boltzmann constant,Ṁ is the mass accretion rate, and M BH is the central BH mass. The model disk spectrum is a superposition of the multi-temperature blackbody components, where superposition produces L ν ∝ ν 1/3 with a high energy turnover to the Wien's tail L ν ∝ e −hν/kT associated with the hotter inner edge of the disk and a low energy turnover to the Rayleigh-Jeans tail L ν ∝ ν 2 associated with the colder outer edge of the disk.
Since the HST SED at t ≃ 1453 days was close to the Rayleigh-Jeans tail, the emission would mostly be from the outer edge of the disk, and the fit cannot constrain the inner edge radius (and thus the BH mass). Therefore, we considered a range of BH masses M BH = 1 − 10 6 M ⊙ and fiducial values of the inclination angle θ = 0 • , 40 • , 80 • . From this, we derived a range ofṀ and T out = T (R out ) from the fits, where R out is the outer edge radius. Note that we do not fit the SED at t ≃ 703 days because it is essentially at the Rayleigh-Jeans tail, meaning that even the outer edge radius cannot be well-constrained.
We also utilized the Swift-XRT upper limit from t ≃ 1360 days to place additional constraints on the inner region of the disk which, depending on the configuration, could be hot enough to produce bright X-ray emission. These constraints can have important implications for the BH mass because brighter X-ray emissions are expected from a smaller central mass with a smaller disk inner edge radius (and vice versa). With the upper limit count rate (0.00157 cts/s), we used the xspec fakeit function and the built-in tbabs×diskbb model to derive the corresponding upper limitṀ and L 0.3−10keV considering the same range of M BH and fiducial values of θ. For fakeit, we used the standard XRT response files: swxs6 20010101v001.arf and swxpc0to12s6 20130101v014.rmf. For tbabs, we assumed a N H = 0.05 × 10 22 cm −2 , similar to AT 2018cow (Margutti et al. 2019). For diskbb, we assumed the inner edge radius to be the innermost stable circular orbit for a non-rotating BH and derived the apparent radius following Kubota et al. (1998).
In Figure 10, we show example disk blackbody fits to the latest HST SED assuming θ = 0 • with M BH spanning six orders of magnitude. We also show a shaded region indicating the upper limits derived from the Swift-(a) Extrapolations with θ as a fixed parameter (b) Extrapolations with θ as a free parameter Figure 11. Extrapolated disk blackbody models at t ≃ 703 days (dashed lines) and t ≃ 1119 days (dotted lines) based on the best-fit models at t ≃ 1453 days (solid lines). The observed (dereddened) HST SEDs are also shown for reference. Top and bottom panels show models for MBH = 10 M⊙ and MBH = 10 6 M⊙, respectively. The left panels assume θ = 0 • for the best-fit models with θ fixed for the extrapolations. The right panels assume θ = 80 • for the best-fit models with θ being a free parameter for the extrapolations and θ fit being the best-fit value that match the observed amplitude.
XRT non-detection. Note that the limits shown on this plot are for illustrative purposes, calculated assuming L ν = L 0.3−10keV /ν with hν = 5.15 keV. The disk blackbody can fit the latest HST SED reasonably well for the entire range of M BH . As expected, models with different M BH are indistinguishable in the optical but predict entirely different X-ray brightness. In particular, the best-fit models predict that the inner regions of the disk around a stellar-mass BH should produce extremely bright X-ray emission: L ν ∼ 10 25 − 10 26 erg s −1 Hz −1 over 0.3 − 10 keV or νL ν ∼ 10 43 − 10 44 erg s −1 at 5.15 keV. This prediction contradicts the Swift-XRT non-detection, which sets an upper limit at L ν ≲ 10 22 erg s −1 Hz −1 over 0.3 − 10 keV. The Swift-XRT non-detection is more consistent with best-fit models that involve IMBHs or SMBHs with much larger inner edge radii that predict much fainter X-ray emission. We discuss this further in Section 4.5.3.

Predicted Evolution of a Remnant Accretion Disk
With the best-fit models, we checked to see if the predicted evolution of an accretion disk could explain the evolution of the observed SEDs. Specifically, we followed the predictions of R out ∝ t 2/3 andṀ ∝ t −4/3 from selfsimilar solutions (e.g., Metzger et al. 2008) and extrapolated the best-fit R out andṀ to the earlier HST epochs, t = 703 days and 1119 days. We built the extrapolated disk blackbody models based on these parameters.
We show some example comparisons between extrapolated disk blackbody models and observed SEDs in Figure 11. The best-fit models at t = 1453 days are shown as solid lines while the extrapolated models are shown as dotted lines at t = 1119 days and dashed lines at t = 703 days. We show two cases at the two ends of the mass range, M BH = 10 M ⊙ and M BH = 10 6 M ⊙ , for comparison at the top and bottom panels, respectively. We perform two sets of fits: one where the inclination angle of the disk is fixed between epochs (examples shown in the left panels with θ = 0 • ) and one where we allow the angle to vary between epochs which changes the model normalization (examples shown in the right panels; starting from θ = 80 • ), a situation that might correspond to a precessing disk. Fitting the inclination angle yields the required precession to explain the SED amplitude at each epoch. When θ is fixed, the extrapolation predicts a lower NUV-optical brightness at earlier epochs (left panels of Figure 11) because the smaller R out shifts the turnover (i.e., the outer edge blackbody) to a shorter wavelength, which has a greater effect than the largerṀ . Although we only show two cases in Figure 11, we found this behavior to hold regardless of M BH and θ. This prediction contradicts the observed HST SEDs, which show a brighter NUV-optical emission at earlier epochs. Therefore, pure extrapolation of R out andṀ cannot explain the evolution of the HST SEDs of the underlying source.
On the other hand, if θ can be varied, the extrapolated models can actually match the observed SEDs quite well if the disk was less inclined (by ∼2 • -8 • ) at the earlier epochs for the case of θ = 80 • at t = 1453 days (right panels of Figure 11). In particular, the color of the extrapolated models are very consistent with the observed color at t ≃ 703 days. Note that the scenario with varying θ only works for high inclination angles because in this case, a small change in θ can cause an appreciable change in cos θ to change the model amplitude. A high inclination angle is an interesting point in the context of AT 2018cow because this was also proposed to explain properties such as the receding photosphere and the asymmetric line profiles (Margutti et al. 2019).

Comparison of Fit Parameters
In Figure 12, we show the parameter spaceṀ vs. M BH with allṀ values derived in our analyses: from the best-fit models to the latest HST SED (blue lines) and the Swift-XRT upper limit (black lines) considering θ = 0 • , 40 • , 80 • (solid, dashed, dotted). Additionally, we show the Eddington accretion rateṀ Edd = L Edd /ϵc 2 (orange lines) considering fiducial values for the radiative efficiency ϵ = 0. 01, 0.1, 0.4 (solid, dashed, dotted). For comparison, a vertical dashed line is also shown to indicate the limit M BH < 850 M ⊙ derived from the NICER QPO (Pasham et al. 2021).
This parameter space highlights the discrepancy between the optical and X-ray in the disk blackbody model for stellar-mass BHs. For M BH ≲ 100 M ⊙ , extremely high (super-Eddington)Ṁ ∼ 10 −2 − 10 M ⊙ yr −1 is required to explain the optical emission, while the Swift-XRT non-detection places a limit atṀ ≲ 10 −4 − 10 −6 M ⊙ yr −1 . The orders of magnitude differences and the required super-Eddington accretion rates are significant challenges to the disk blackbody model if AT 2018cow involved a stellar-mass BH. This discrepancy disappears at M BH ≳ 10 4 M ⊙ , the IMBH and SMBH range, where (sub-Eddington)Ṁ ∼ 10 −4 − 10 −6 M ⊙ yr −1 is required to explain the optical emission. Although this finding is consistent with the scenario that AT 2018cow was a TDE involving an IMBH or a SMBH, this interpretation contradicts the limit M BH < 850 M ⊙ derived by (Pasham et al. 2021) from NICER QPOs.
Additional absorption may be another way of explaining the lack of X-ray. We tested this case by multiplying the original xspec model by phabs and found the required N H such that the best-fitṀ from fitting the HST SED can produce the upper limit Swift-XRT count rate (0.00157 cts/s). For M BH ∼ 10 − 10 3 M ⊙ , we found that significant absorption with N H ≳ 10 24 − 10 25 cm −2 is required explain the lack of X-ray. Such high absorption has been observed for some Galactic X-ray binaries (e.g., Matt & Guainazzi 2003) and Compton-thick active galactic nuclei (e.g., Ricci et al. 2015;Marchesi et al. 2018). We do note that our estimated N H is highly uncertain since the X-ray spectral shape and emitting mechanism are essentially unconstrained.

Summary of Implications For Remnant Disk
Here we summarize our main findings about the viability of the accretion disk model for the underlying source and implications for various transient models.
• An accretion disk can naturally explain the high temperature and small size of the underlying source, and a disk blackbody can reasonably fit the observed HST SED at t = 1453 days.
• The predicted evolution of an accretion disk (i.e., R out ∝ t 2/3 andṀ ∝ t −4/3 ) is consistent with the color evolution of the underlying source, but inconsistent with the brightness/amplitude evolution. We found that if the inclination angle is large (close to edge-on) at t = 1453 days, then a precessing disk with an increasing inclination angle can explain the brightness/amplitude evolution.
• For stellar-mass BHs (M BH ≲ 100 M ⊙ ), a super-Eddington accretion rate is required to explain the optical emission but cannot explain the lack of Xray emission. For IMBHs and SMBHs (M BH ≳ 10 4 M ⊙ ), both the optical emission and the lack of X-ray emission can be explained by sub-Eddington accretion rate, but the central mass would violate the limit derived from the NICER QPO, M BH < 850 M ⊙ (Pasham et al. 2021).
In the context of transient models, precessing accretion disks are often considered in TDE frameworks where the disk angular momentum vector is misaligned with the BH spin axis (e.g., Shen & Matzner 2014;Bonnerot et al. 2016;Hayasaki et al. 2016;Liska et al. 2018), and the precession is due to the Lense-Thirring effect. Intriguingly, the same effect was proposed to explain the zero time lag between the optical and X-ray QPO in AT 2018cow (Li 2022). Recent TDE simulations have also shown that high inclination angles can lead to significant X-ray reprocessing and bright UV-optical emission (Dai et al. 2018;Thomsen et al. 2022). Latetime UV-optical emission have also been observed for TDEs (van Velzen et al. 2019), interpreted as emission from viscously spreading accretion disks. In this context, the precessing highly-inclined accretion disk scenario appears to be compatible with the observational constraints of the underlying source.
In terms of the precessing timescale, as an example, if we consider the Bardeen-Petterson effect (Bardeen & Petterson 1975;Hatchett et al. 1981;Nelson & Papaloizou 2000;Fragile et al. 2001) and derive the alignment timescale (Scheuer & Feiler 1996), we find t align ≈ 147 α where α 0.1 = α/0.1 with α being the viscosity parameter, δ 0.5 = δ/0.5 with δ being the disk aspect ratio, and a * ,0.1 = a * /0.1 with a * being the dimensionless specific angular momentum of the BH. The rate of precession is given by 2π/t align , which following the equation above can be a few degrees per year for M BH ∼ 10 M ⊙ anḋ M ∼ 10 −2 M ⊙ yr −1 . However, note that the assumed values for the dimensionless parameters are uncertain. Perhaps the largest uncertainty in the accretion disk scenario for the underlying source is the mass of the BH involved, which has significant implications on the progenitor system of AT 2018cow. A stellar-mass BH would likely imply a stellar progenitor for AT 2018cow (core collapse of single star or binary merger with a BH) but would have to maintain super-Eddington accretion for many years to explain the observed underlying UVbright emission and not produce bright X-ray emission. We note that if we follow the wind-reprocessed framework for prompt emission of AT 2018cow as discussed in paper I (Chen et al. 2023) and assume that the accretion rate is on the order ofṀ ∼ 10 −1 − 10 1 M ⊙ yr −1 over the first two months post-discovery, then a constant decline ofṀ ∝ t −4/3 would predict an accretion rate on the order ofṀ ∼ 0.5 − 1.5 × 10 −3 M ⊙ yr −1 at t = 1453 days. This is inconsistent with the accretion rates derived from our analysis for stellar-mass BHs using a disk blackbody model (Ṁ ∼ 0.1 − 1 M ⊙ yr −1 ; see Figure 12). We also note that while the stellar-mass BH scenario may produce significant precession over a few years (Equation 3), it is unclear if core collapse can naturally lead to a misalignment between the remnant disk and the spin of the compact object in the first place to cause the precession.
However, since the thin disk approximation does break down for super-Eddington accretion, the disk blackbody model is likely not appropriate for accretion around stellar-mass BHs and our derived parameters may not be entirely accurate. Properly accounting for the super-Eddington accretion may also help explain the lack of X-ray emission through mechanisms such as the obscuration by a geometrically thick disk or dense outflow (Done et al. 2007). For example, Metzger (2022) constructed a model to explain AT 2018cow involving super-Eddington accretion from a delayed binary merger between a Wolf-Rayet star and a BH or NS. They predicted the late-time evolution of the accretion disk (Section 2.3.3 in Metzger 2022) with the accretion rate decreasing with radius as a power law due to outflow carrying away angular momentum, which may lower the inner accretion rate by a few orders of magnitude and explain the lack of X-ray. We note that Metzger (2022) also predicted late-time thermal emission from the outer accretion disk, and while their model radius and temperature appear similar to those observed for the underlying source, their model luminosity (being the Eddington luminosity L Edd ∼ 10 38 M BH erg s −1 ) would require M BH ≳ 10 4 M ⊙ (i.e., an IMBH) to match the derived blackbody luminosity for the underlying source at t ≃ 703 days (Table 2). Overall, additional theoretical works are needed to constrain the viability of the accreting stellar-mass BH scenario as an explanation of the underlying source.
On the other hand, a geometrically thin accretion disk around an IMBH or a SMBH with sub-Eddington accretion rate (Ṁ ≲ 10 −2 M ⊙ yr −1 ) can reasonably explain the UV-bright underlying source. These accretion rate implied for an IMBH would also be roughly consistent with the wind-reprocessed framework for the prompt emission of AT 2018cow that followṀ ∝ t −4/3 , though the accretion rate may not be enough to cause significant precession (Equation 3). The IMBH or SMBH scenario would imply that AT 2018cow was a TDE. Previous studies have disfavored the TDE hypothesis for various reasons, most notably the difficulty in explaining the dense CSM (Margutti et al. 2019;Huang et al. 2019), the existence of such BH at the outskirt of the galaxy without any signs of a coincident massive host system (Lyman et al. 2020), and the mass limit M BH < 850 M ⊙ derived from the NICER QPO (Pasham et al. 2021). However, some evidence may suggest otherwise. For example, environmental studies (Roychowdhury et al. 2019;Lyman et al. 2020) have noted that a faint tidal tail in the host galaxy (also see Figure 1) that traces star-forming activities around AT 2018cow could be evidence of recent dynamical interaction. One may speculate that there could be a chance that a straggling IMBH or SMBH was left behind from this dynamical interaction. Another example is the low-frequency Xray QPO discovered by Zhang et al. (2022) from XMM-Newton and Swift observations. Zhang et al. (2022) suggested that the low-frequency QPO is more consistent with IMBHs and SMBHs, and the NICER QPO limit could be relaxed by introducing the IMBH or SMBH in a compact binary. Therefore, the TDE scenario for AT 2018cow may be worth revisiting along with the hypothesis that the underlying source was an accreting IMBH or SMBH.

SUMMARY & CONCLUSION
In this study, we examined the UV-bright transient underlying source at the precise position of AT 2018cow revealed by the three HST observations taken ∼2-4 years post-discovery (t ≃ 703, 1119, 1453 days). The HST observation at t ≃ 1453 days, which we requested after independently discovering the underlying source, showed significant fading in the UV bands relative to the observations at t ≃ 703 days ( Figure 5). This establishes the transient nature of the source, which could be cause either by an intrinsic (i.e., emission associated with AT 2018cow) and/or extrinsic (i.e., increased absorption along the line of sight) effect.
The underlying source is bright (L UVO,min ∼ 10 39 erg s −1 ) and exceptionally blue (F336W−F555W = −1.3) with an unconstrained peak further in the UV (λ peak ≲ 2358Å). The blue spectrum at t ≃ 703 days can be described by a spectral index of α = 1.99 (similar to the Rayleigh-Jeans tail) or by a blackbody with a high temperature (T ≳ 10 5 K) and a small radius (R ≲ 20 R ⊙ ). A flatter spectrum at t ≃ 1453 days with α = 1.66 could be an indication of of cooling (to T ∼ 6×10 4 K) and expansion (to R ∼ 30 R ⊙ ) of a blackbody, or an increase in extinction (with a color excess of E B−V ≃ 0.072) assuming the Cardelli extinction law (Cardelli et al. 1989) with R V = 3.1.
We considered five origins of the properties and evolution of this peculiar UV-bright underlying source: (i) significant contribution from a star cluster, (ii) increased extinction from newly-formed dust along the line of sight (iii) ejecta-CSM interaction, (iv) magnetar spin down, and (v) a remnant accretion disk around a BH. We disfavored significant contribution from a star cluster based on comparisons with BPASS and LEGUS clusters because for this scenario to work, the underlying source has to contain both an extremely young cluster and a transient source bluer than the Rayleigh-Jeans tail. We found that although dust formation appears reasonable in the context of AT 2018cow, the fading was unlikely purely due to dust extinction because of the extremely blue color already observed.
We additionally ruled out ejecta-CSM interaction involving the known radio-producing CSM from modeling the expected radiation, and magnetar spin down with B ∼ 10 15 G based on the energy output. However, we cannot rule out ejecta-CSM interaction involving a denser CSM component (e.g., a previously ejected hydrogen-rich envelope) or magnetar spin down with B ≲ 10 14 G, and additional modeling would be required to constrain these possibilities.
Finally, we found that a precessing accretion disk at a high inclination angle can reasonably explain the color, brightness, and evolution of the HST SEDs. However, a major uncertainty is the type of BH at the center of the accretion disk. A stellar-mass BH would require super-Eddington accretion over multiple years with an accretion rate possibly declining slower than the pre-dictedṀ ∝ t −4/3 and additional mechanisms to explain the lack of X-ray emission. On the other hand, while an IMBH or a SMBH could naturally explain the lack of associate X-ray emission with an inferred accretion rate consistent with the wind-reprocessed framework for AT 2018cow (paper I; Chen et al. 2023) that followṡ M ∝ t −4/3 , this would appear to violate the limit of M BH < 850 M ⊙ from the NICER QPO (Pasham et al. 2021) and its existence at the location of AT 2018cow would still be difficult to explain.
Putting together all the pieces, including results from paper I (Chen et al. 2023), we find that central engine and ejecta-CSM interaction are still the preferred power sources that can coherently explain both the luminous FBOT AT 2018cow and the remnant UV-bright slow-evolving transient underlying source. However, we note that for ejecta-CSM interaction to fully explain the observations, multiple CSM components are necessary: (I) a dense CSM shell (Ṁ ∼ 1 M ⊙ yr −1 ) is required to power the fast-rising luminous peak of AT 2018cow at t ∼ 1 day (Xiang et al. 2021;Pellegrino et al. 2022), (II) dense aspherical CSM (unknown density and distribution) is required to sustain the optically thick rapidlyfading prompt emission over t ∼ 20 − 60 days, (III) relatively less-dense CSM (Ṁ ∼ 10 −6 − 10 −4 M ⊙ yr −1 ) is required to power the radio emission up to t ∼ 600 days (Ho et al. 2019;Nayana & Chandra 2021), and (IV) dense and extended CSM (unknown density and distribution) farther away from the transient (R ≳ 10 17 cm), likely the previously ejected hydrogen envelope, is required to power the UV-bright slow-fading transient over t ∼ 700−1500 days. While these CSM components could hypothetically exist, all but the radio-producing CSM are poorly constrained both observationally and theoretically. Additional theoretical works, and perhaps followup observations for the underlying source, can help constrain or rule out possible CSM components.
In the context of our analyses, we favor the accreting BH scenario because some of the expected phenomena from this scenario can reasonably explain the observations of AT 2018cow and the underlying source. As shown in paper I (Chen et al. 2023), the fading prompt emission and the associated peculiar thermal properties can be explained by continuous wind outflow driven by an accreting central engine, and as argued in this paper, the evolution of the remnant accretion disk can naturally give rise to the underlying source. This would support the hypothesis that AT 2018cow and the class of luminous FBOTs may form entirely new class of BH transients powered predominantly by accretion. However, there are challenges faced by this interpretation, such as the requirement of disk precession and the uncertainty in the mass of the BH involved. These challenges could either be counterarguments to this hypothesis or interesting constraints that may contain new information on AT 2018cow and BH transients. Late-time evolution of transient accretion disks could be an interesting topic for future theoretical studies not only to examine the underlying source of AT 2018cow but also to explore potential transient phenomena related to accreting BHs and potentially reveal new classes of UV-bright BH transients similar to the underlying source.
Our studies highlight the importance of late-time observations, which for AT 2018cow, provided significant constraints on the late thermal properties (paper I; Chen et al. 2023) and led to the discovery of an unprecedented underlying transient for an FBOT years post-discovery (this paper). These merits justify similar late-time observations and monitoring for nearby peculiar transients through powerful telescopes such as the HST and JWST. For new discoveries of nearby "Cow-like transients", late-time observations will be crucial in probing remnant emission possibly associated with ejecta-CSM interaction or an accretion disk and if the host galaxy experienced similar dynamical interaction that could leave behind a straggling IMBH or SMBH.
Lastly, we mention that our study demonstrates the need for next-generation UV telescopes because a UVbright transient such as the underlying source, which was not explicitly predicted by previous studies for AT 2018cow, was barely recognized through the HST. Without NUV monitoring from the HST, the underlying transient source would have been completely missed based on the optical and (the lack of) X-ray emission. Next-generation UV telescopes will be crucial for expanding the transient phase space and exploring potential new UV transients such as the underlying source of AT 2018cow.