
\section{List of Figures}

\section{Introduction} \label{sec:intro}
Massive stars ($\gtrsim$8 M$_{\odot}$) end their life in core collapse after nuclear burning is no longer sustainable.  Often the collapse of a massive star results in the ejection of the stellar envelope and a luminous core-collapse supernova (SN), but this may not always be the case. Indeed, one surprising claim is that standard type IIP SNe, for which progenitor stars have been identified in deep pre-explosion data, arise from red supergiant (RSG) stars with only modest initial masses $\lesssim$18 M$_{\odot}$ \citep[e.g.,][]{smartt09,smartt15},  even though RSGs as a population extend to $\gtrsim$25 M$_{\odot}$ \citep{humphreys79}. The apparent lack of high-mass RSG progenitors to SN IIP suggests that either they are systematically missed by transient surveys (e.g., they explode preferentially in heavily obscured regions) or that the highest mass RSGs do not end their life as SNe. However, the value of the mass cutoff at $\sim$18\,$M_{\odot}$, its statistical validity and its interpretation are controversial due to the complexities of connecting the limited observations to uncertain models \citep[e.g.,][]{davies07,smith11,walmswell12,kochanek12,davies18a,davies20a,davies20b,kochanek20}.


The putative lack of SN type II progenitors with initial masses $\gtrsim$18 M$_{\odot}$ has led to the hypothesis that more massive stars cannot yield successful SN explosions, instead ending their lives as ``failed'' SNe that collapse directly into black holes \citep[e.g.,][]{kochanek08}, possibly with an associated fast, blue \citep{kashiyama15} or faint, $\approx$yr-long transient \citep{piro13,lovegrove13}. Producing simulations that yield successful high-energy SN explosions has been a long-standing difficulty.  Modern simulations and theoretical work have found general agreement with this picture of a large fraction of failed SNe from high-mass RSGs, in that a majority of lower-mass RSG stars have compact cores that are more likely to produce successful explosions, while more massive ($\gtrsim$16--$20\,M_{\odot}$) stars may be more likely to collapse into black holes \citep[e.g.,][]{oconnor11,ugliano12,sukhbold14,pejcha15,ertl16,sukhbold16}. Though not strictly monotonic in mass, the general trend may be related to a transition in the star's core between convective and radiative carbon burning near the end of its life \citep{sukhbold19}.

Verification that massive RSG stars can directly collapse into black holes without an accompanying SN would shed light on the observed compact remnant population \citep{kochanek14} and point to the stellar-mass progenitors of the black hole mergers observed by gravitational wave experiments \citep{lvc_O1O2,GWTC2}.  The ideal observational signature of a failed SN would manifest as a bright RSG star that winks out of existence (or nearly so) from one epoch to the next, and stays in that state permanently.  With no accompanying luminous SN, the only explanation for the disappearance of such a massive star would be a direct collapse into a black hole, although some fading, optical/near-IR emission may remain from a weak envelope ejection or fallback accretion onto the newly formed black hole \citep{lovegrove13,perna14}.  

Multiple projects have searched for disappearing massive stars.  The first and most comprehensive of these has utilized the Large Binocular Telescope (LBT) to re-image nearby galaxy fields ($\lesssim$10\,Mpc) with a cadence of $\approx$months over the last decade \citep{gerke15,adams17a,neustadt21}. This search has identified the best known candidate for a failed SN, dubbed NG6946-BH1 \citep{adams17b}, which was an apparent $\sim$25 M$_{\odot}$ ($\sim$10$^{5.5}$ L$_{\odot}$) RSG which brightened on a year timescale before nearly vanishing in the optical and infrared (IR). Follow-up observations show some fading remnant emission at the position, consistent with expectations of fallback accretion \citep{adams17b,basinger20}. The ongoing LBT survey has identified other, potential failed SN candidates, including a blue supergiant star in M101 \citep{neustadt21}. A separate, multi-epoch archival \textit{Hubble Space Telescope} (\textit{HST}) search also identified a yellow supergiant (YSG) star that apparently vanished \citep{reynolds15}. 




A primary challenge in identifying true failed SNe is that massive, evolved stars are subject to several forms of large-amplitude variability on essentially all timescales, some which may mimic a disappearing star at optical wavelengths. Eruptive variability, like that of the luminous blue variables (LBVs; e.g., \citealp{humphreys99,smith01b,smith17}) or massive stellar mergers can form dust in their dense outflows, obscuring the surviving star in the optical, though producing bright mid-IR emission (e.g., \citealp{smith16,smith18,blagorodnova17,pejcha16a,pejcha16b}, see also, IR-dominated transients in \citealp{jencson19}). 

Massive, cool supergiants, while subject to pulsational instabilities \citep{heger97,yoon10} manifesting as periodic or semiregular variations on timescales of months to years \citep[e.g.,][]{kiss06,yang11,yang12,soraisam18}, can also undergo episodes of enhanced mass loss accompanied by observed LBV-like color changes \citep{smith04a}. These episodes may produce complex circumstellar material (CSM) like that observed for nearby, Galactic Y/RSGs that suggests prior, distinct mass ejections \citep[e.g.,][]{humphreys97,smith01a,smith04b,monnier04}. During such an event, increased opacity from a dense molecular wind (as in Mira variables, \citealp{reid02}; and for recent work on RSGs see, e.g., \citealp{davies21}) or the formation of dust can obscure star causing a dramatic visual dimming. Such a scenario has, for example, been discussed at length to explain the recent historic minimum of Betelgeuse that occurred in 2019--2020, known as the ``Great Dimming'' \citep[e.g.,][]{levesque20,montarges21}. Modern extragalactic transient searches are designed to find bright, explosive outbursts that can be seen at large distances. Dramatic fading events may also signify important instabilities in massive evolved stars, however, as with failed SNe, they are more challenging to find and our knowledge of their impact on stellar evolution is limited. 

In this paper, we describe the first results from a program to search for failed SNe using a set of new \textit{HST} data directed at nearby galaxies that have both a long history of archival \textit{HST} imaging, and have hosted multiple core-collapse SNe. We have identified a massive, cool supergiant star in M51 which underwent a significant dimming in 2019, preceded by a roughly year-long brightening phase, both of which pointed to a potential failed SN. As the first compelling candidate from our search, we named the object M51-DS1 ($=$ M51 ``Dimming Star'' 1)
and began a dedicated follow-up campaign. Our subsequent \textit{HST} and ground-based near-IR imaging showed that the star had rebrightened in 2021 ruling out a failed SN explanation for the dimming event. Nevertheless, we take a detailed look at the available data for the object as both a false-positive in our failed SN search and as a prime example of exceptional variability for a massive, evolved supergiant. 


\section{A New \textit{HST} Search for Failed SNe}\label{sec:FSN_search}

As described above, the clear observational signature of a failed SN is the permanent disappearance of a massive star. In order to observe more examples of the failed SN phenomenon, we chose nearby galaxies that are both proven core-collapse SN producers and have multiple epochs of \textit{HST} imaging over the last $\sim$30 years.  To do this, we first identified galaxies that have hosted two or more core-collapse SNe in the modern era using the Asiago Supernova Catalogue \citep{asiago}, supplemented by SNe found in the last several years using a query of the Transient Name Server\footnote{https://www.wis-tns.org/}. 
We cross-matched this list of galaxies with those that had two or more existing epochs of \textit{HST} data in the $F814W$ filter, and we reobserved the resulting 31 galaxies (all with $D\lesssim$45 Mpc) over 41 orbits with the Advanced Camera for Surveys (ACS) in Cycle 26 (PI: D.\ Sand, PID: 15645), matching the footprint of previous \textit{HST} imaging as best as possible.  The $F814W$ filter was chosen as our discovery filter because it is among the most sensitive to the supergiant progenitors to failed SNe, and has been commonly used throughout \textit{HST}'s lifetime with ACS, the Wide Field Planetary Camera 2 (WFPC2), and the Wide Field Camera 3 (WFC3); $F814W$ is also less sensitive to extinction in regions of star formation than bluer bands.  Each newly acquired $F814W$ image is a full orbit with a standard sub-pixel dither sequence, and has a sufficient depth ($F814W\approx27.2$\,mag with a signal-to-noise ratio [S/N] of $\approx$5) to observe massive stars with initial stellar masses $\gtrsim$10 M$_{\odot}$ throughout the distance range we cover, although some of the archival imaging we use is only sensitive to initial stellar masses of $\gtrsim$15 M$_{\odot}$.
 
We identify massive star candidates that have vanished using difference imaging, subtracting older \textit{HST} epochs with our newly acquired data. We outline the general procedure in Section~\ref{sec:obs} as we describe the identification of M51-DS1.  Key to our strategy is the wide range of data available in both the \textit{HST} and ground-based data archives for each of our target galaxies, which allows us to establish a baseline of existence for the progenitor star of any given failed SN candidate stretching back years prior to disappearance.  We will provide further details about our strategy and expected number of failed SNe in a future work, as analysis is ongoing.

\section{Observations and Data Processing} \label{sec:obs}

Here we describe the observations and data processing associated with the identification of M51-DS1 as a failed SN candidate, along with subsequent follow-up observations.

\subsection{HST Image Processing and Subtraction}
To search for large amplitude variables and ``disappearing'' luminous stars, we performed image subtraction with the available archival
WFC3/UVIS and ACS/WFC images in the $F814W$ filter that overlap with the footprint of our Cycle 26 program observations. For M51, this includes
ACS/WFC frames taken from 2005 January 12--22 (PI: S.\ Beckwith, PID: 10452) and a large campaign consisting of 34 epochs taken between 2016 October 5 and 2017 February 17 to study massive star variability (PI C.\ Conroy, PID: 14704; see \citealp{conroy18}; hereafter C18). 

We used the Cycle 26 ACS/WFC $F814W$ observations of M51 as template images for subtraction with each of the archival frames. These include two spatially overlapping visits taken on 2019 May 24 and 31, respectively, designed to maximize coverage of M51 and overlap with the archival images. Each visit consisted of four dithered 564~s exposures. We first downloaded the \texttt{CALACS} calibrated and charge-transfer-efficiency (CTE)-corrected \texttt{flc} frames for both visits from the Mikulski Archive for Space Telescopes. We processed the images using the \texttt{AstroDrizzle} software package, including automated cosmic ray rejection, sub-pixel alignments with \texttt{TweakReg}, and final combination into a single, drizzled mosaic at the native WFC pixel scale of $0\farcs05$ with an effective PSF (ePSF) of $\approx2$\,pixels. 

We then processed the available archival $F814W$ frames from ACS/WFC in a similar manner. 
The 2005 observations consisted of a 6-pointing, mosaicked mapping of M51, with each pointing composed of four dithered, individual 340\,s exposures. We aligned each of the 24 individual frames to the drizzled Cycle 26 template mosaic using \texttt{TweakReg}, typically using between $\approx$300 and $4000$ common stars depending on the size of the overlap region with the template.
For the 2016--2017 data, each of the 34 epochs consists of 4 individual exposures totalling 2200\,s. We again aligned each frame to the template using \texttt{TweakReg}, where the number of common stars used for alignment varied between $\approx$100-500 stars. We achieved ($\mathrm{rms} \lesssim 0.1$ WFC pixels in the $x$ and $y$ directions) for every frame. We then ran the aligned images through \texttt{AstroDrizzle}. For the purposes of image subtraction, we produced one large drizzled mosaic of the 2005 images and split the 2016--2017 into three batches consisting of observations taken between 2016 October -- 2017 January, 2017 February--May, and 2017 June--September. We used the 2019 template image as a reference in \texttt{Astrodrizzle} so that all the output images were drizzled onto the same pixel grid. 

To prepare the drizzled images for subtraction, we ran each through the Astromatic packages \texttt{SExtractor} and \texttt{PSFEx} \citep{bertin96,bertin11} to generate source catalogs and models of the point-spread function (PSF) for each image. The images and corresponding noise maps generated from the inverse-variance-map output of \texttt{AstroDrizzle} were then converted to units of electrons on the scale of the 2019 template using the zeropoint and exposure time information in the image headers. The background level of the images were also estimated and removed using the minimum value from a median filter on a 15$\times$15 grid of the region of overlap with the template. 

We then performed image subtraction using the ZOGY algorithm \citep{zackay16} on the background-subtracted images, with the noise maps and PSF model postage stamps as additional inputs. The subtraction produces a difference image, the difference image PSF, and a Scorr image (Equation 25 in \citealp{zackay16}), which is a match-filtered S/N image optimized for point source detection. In order to use the Scorr image for the detection of variable objects, we normalize the Scorr values down to have a standard deviation of 1 over the image. This is necessary to account for the contribution of correlated pixel noise inherent in our drizzled images, which is not included in the calculation of the Scorr image with ZOGY. Finally, SExtractor is run on both the ``positive'' (archival-minus-template) and ``negative'' (template-minus-archival) difference images to build source catalogs of both fading and brightening objects. These are cross-matched with the corresponding, normalized Scorr image to select significant detections of $\mathrm{S/N} \geq 5$ (those objects with peak normalized Scorr values $\geq$5 in at least one pixel). 

\subsection{Candidate Identification}\label{sec:candidates}
To select initial candidates, we require that a source be detected at $|\Delta \lambda L_{\lambda}| \gtrsim 10^4~L_{\odot}$ in the ``positive'' difference images from at least two archival $F814W$ frames compared to the 2019 template, corresponding to an absolute magnitude in the difference image photometry of $M_{F814W} \leq -6.38$\,mag (Vega system) and indicating a luminous object that has faded significantly. This cut-off in the change in luminosity is the same as that used for the LBT search \citep{gerke15,adams17a,neustadt21}, which is sufficient to exclude lower mass large-amplitude variables including R CrB stars and Miras \citep[$M_I>-6$\,mag,][]{tisserand09,soszynski09}. We also require that the two detections be separated in time by at least 1\,month. This establishes an archival baseline for a given source as a persistent object, i.e., to reject real astrophysical transient events such as novae as well as chance cosmic ray hits and other image processing and subtraction artifacts.

We then visually inspect the position of the object in all available $F814W$ \textit{HST} imaging to further reject any obvious image artifacts and examine the full history of the source in that band. A strong failed SN candidate will appear as a star-like source that is present in all archival frames of sufficient depth prior to its first significant fading, after which the source will not appear to rebrighten. We note that we do not require that a source actually completely ``disappear'', because a possible companion or unrelated, nearby star in a crowded region may continue to be detectable at the location even after a true failed SN. Any source that passes this visual inspection will be assigned a ``DS'' name (for ``Dimming Star'') and flagged for detailed analysis, including PSF-fitting photometry of the available \textit{HST} imaging in other filters (see Section~\ref{sec:dolphot} for details). 

A subset of the sequence of available \textit{HST} $F814W$ imaging for M51-DS1, our first compelling candidate (as defined by the criteria outlined above), is shown in the bottom row of panels in Figure~\ref{fig:hst_im}. The star was clearly detected in the 2005 and 2016 ACS/WFC images (also in 1995 and 2008 WFPC2 images not shown). In the 2019 image obtained as part of our search the star has seemingly vanished and is clearly detected as a fading candidate in our difference images. 

\begin{figure*}
\centering
\includegraphics[width=0.7\textwidth]{NGC5194-DS1_HST_imaging.pdf}
\caption{\label{fig:hst_im}
HST imaging of M51-DS1 and the surrounding region. From left to right along the bottom row of panels, we show the sequence of ACS/WFC $F814W$ images from 2005, 2016, 2019, the 2016--2019 difference image, and the most recent 2021 image. The source underwent a dramatic fading, nearly disappearing in the 2019 frame (highlighted in the red-bordered panels), before rebrightening in 2021. The large, upper-left panel shows the color-composite ACS/WFC images from 2005 ($F435W$ in blue, $F555W$ in green, and $F814W$ in red) of the location of M51-DS1 along an inner spiral arm and prominent dust lane $\approx$36$\arcsec$ east of the galaxy nucleus. The positions of the nearest two \ion{H}{2} regions from \citet{croxall15} are indicated with yellow circles of 1$\arcsec$ radius and labelled with the derived $E(B-V)$ values for each one. A zoomed-in view of the immediate $1\farcs5 \times 1\farcs5$ (62.4~pc on a side) region around M51-DS1 is shown in the upper-rightmost panel. Below this, down the right side of the figure, we then show the same region in the ACS/WFC $F606W$ images from 2016 and 2021, and the WFC3/IR $F110W$ images from 2012 and 2021, demonstrating the dimming of the source across visible and near-IR wavelengths.
}
\end{figure*}

\subsection{Archival imaging, Follow-up Observations and Photometry}\label{sec:imaging}


Following our identification of M51-DS1 as a disappearing star and promising candidate for a failed SN, we began examining the wealth of multiband, archival \textit{HST} imaging, including the ACS/WFC $F814W$ images described above and the complimentary $F435W$, $F555W$, and $F658N$ (2005) and F606W (2016--2017) from the same programs. There is also coverage with WFC3/IR ($0\farcs13$ pixels) in $F110W$ and $F128N$ (Pa$\beta$) taken on 2012 September 4 (PI: J.\ Koda, PID: 12490). As shown in Figure~\ref{fig:hst_im}, M51-DS1 is notably red in appearance in the available archival imaging, and the optical source appears to be associated with a bright, relatively isolated star in the near-IR frames. Of the available WFC3/UVIS observations, we use only the 2012 observations in the redder $F673N$ and $F689M$ filters (PI: K.\ Kuntz, PID: 12762) where the star is clearly identifiable in the drizzled images. Similarly for WFPC2, while there are numerous prior observations in many filters that cover the location on the WF cameras, we consider only the $F814W$ images from 1995 (PI: R.\ Kirshner; PID: 5777) and 2008 (PI: M.\ Meixner, PID:  11229) given the red color of the star and the increasing likelihood of contamination in bluer filters for the larger pixels of the WF chips ($0\farcs0996$).

We were also awarded an \textit{HST} Cycle 28 mid-cycle proposal for additional deep, optical and near-IR imaging of the field (PI: J.\ Jencson, PID: 16508) in order to verify the extreme fading of this source. The images for this program were obtained with ACS/WFC in $F606W$ and $F814W$ (total exposure time of 2208\,s each filter) on 2021 April 28--29 and with WFC3/IR in $F110W$ and $F160W$ on 2021 June 6 (1200\,s exposure per filter). These new images revealed a surviving star at the location of M51-DS1 that had clearly rebrightened at $F814W$ from the 2019 minimum (see the bottom right corner of Figure~\ref{fig:hst_im}). While this rules out a failed SN scenario for this object, the star is notably fainter at each of $F606W$, $F814W$, and $F110W$ compared to the corresponding pre-disappearance imaging. 

\subsubsection{\texttt{DOLPHOT} Photometry}\label{sec:dolphot}
We processed the available \textit{HST} imaging mentioned above with \texttt{DOLPHOT} \citep{dolphin00,dolphin16} to obtain PSF-fitting photometry of M51-DS1 and the surrounding stars in its vicinity. As inputs to \texttt{DOLPHOT}, we use the CTE-corrected \texttt{flc} frames for ACS/WFC and WFC3/UVIS observations, \texttt{flt} frames for WFC3/IR observations, and the standard exposure level \texttt{c0f} and data quality \texttt{c1f} frames for WFPC2. Each input image was also run first through \texttt{Astrodrizzle} to flag cosmic ray hits. We employ the parameter settings used for the \textit{HST} PHAT survey \citep{dalcanton12,williams14}. These settings were also used by C18 to build their photometry catalogs of M51, to which we can directly compare our results. 

For most of the data, we ran \texttt{DOLPHOT} separately for each instrument and filter combination over the available epochs, using the corresponding drizzled images as references for alignment. Generally, \texttt{DOLPHOT} produced excellent alignments better than $\approx$0.2\,pixels rms for ACS/WFC and WFC3/UVIS, and $\approx$0.3\,pixels rms for WFC3/IR and WFPC2. Owing to the large quantity of ACS/WFC $F814W$ and $F606W$ data (primarily from the C18 program), we had to run them in batches. In order to obtain consistent results across all runs for a given filter, we use the ``warmstart'' option in \texttt{DOLPHOT}, in which we use the output star catalog from an initial batch where M51-DS1 is well detected to fix the positions of stars in the subsequent runs for a given filter.  \texttt{DOLPHOT} computes and applies aperture corrections to a radius of $0\farcs5$ for the reported photometry. We then applied the appropriate corrections to infinite apertures for each instrument and filter combination.\footnote{See \citep{bohlin16} for ACS; for WFC3: \url{https://www.stsci.edu/hst/instrumentation/wfc3/data-analysis/photometric-calibration}; for WFPC2: \url{https://www.stsci.edu/itt/review/dhb_2010/WFPC2/wfpc2_dhb.pdf}} 
To estimate the photometric uncertainties, we compute the rms deviations of the measurments for the individual frames that comprise a given \textit{HST} observing visit of all stars within a 100\,pixel radius of M51-DS1 as a function of magnitude, added in quadrature with the nominal statistical uncertainty reported by \texttt{DOLPHOT}. 

For all the data that we ran, a ``good'' star (``object type''$=$1) was detected by \texttt{DOLPHOT} at the location of M51-DS1 as measured in the 2016 $F814W$ ACS frames within $\approx$0.5 ACS pixels, with the exception of the 2005 $F435W$ and $F555W$ ACS images. In these bluer filters, photometry of M51-DS1 is limited by crowding of nearby sources. We therefore adopt the magnitudes of the nearest detect star (1.3 and 0.8 pixels away for $F435W$ and $F555W$, respectively) as limits. Our final \texttt{DOLPHOT} photometry is presented in Table~\ref{table:HST_phot}. Throughout this work, all photometry is reported on the Vega magnitude system.  

Our measurements of M51-DS1 at $F814W$ from the ACS/WFC 2016--2017 data are typically about 0.05--0.1\,mag fainter than those reported for the star in the C18 catalog, though the offset is comparable with our measurement uncertainties. The offset is somewhat larger at $F606W$, with our measurements systematically at $\approx$0.3\,mag fainter that those of C18. In experimenting with \texttt{DOLPHOT} runs, we found that the choice of the alignment reference image as well as the specific images included in a given batch affected the final star list. In our final catalogs, there were stars immediately adjacent to M51-DS1 that were not in the C18 catalog. The simultaneous fitting of these additional stars by \texttt{DOLPHOT} may explain our lower flux measurements. 

\begin{deluxetable*}{ccccc}
\tablecaption{HST \texttt{DOLPHOT} Photometry \label{table:HST_phot}}
\tablehead{\colhead{UT Date} & \colhead{MJD} & \colhead{Inst.} & \colhead{Band} & \colhead{App.\ Magnitude\tablenotemark{a}} \\ 
\colhead{} & \colhead{} & \colhead{} & \colhead{} & \colhead{(mag)} }
\startdata
1995 Jan 15.76 & 49732.76 & WFPC2   & $F814W$ & 22.03 (0.20) \\
2005 Jan 21.36 & 53391.36 & ACS/WFC & $F435W$ & $>24.8$      \\
2005 Jan 21.37 & 53391.37 & ACS/WFC & $F555W$ & $>25.0$      \\
2005 Jan 21.38 & 53391.38 & ACS/WFC & $F814W$ & 22.25 (0.11) \\
2008 Mar 31.14 & 54556.14 & WFPC2   & $F814W$ & 22.79 (0.20) \\
2012 Sep 04.82 & 56174.82 & WFC3/IR & $F110W$ & 20.09 (0.17) \\
2012 Sep 04.83 & 56174.83 & WFC3/IR & $F128N$ & 19.54 (0.09) \\
2016 Oct 05.25 & 57666.25 & ACS/WFC & $F606W$ & 24.51 (0.10) \\
2016 Oct 05.31 & 57666.31 & ACS/WFC & $F814W$ & 22.44 (0.05) \\
2016 Oct 14.25 & 57675.25 & ACS/WFC & $F606W$ & 24.54 (0.09) \\
\enddata
\tablenotetext{a}{Vega magnitudes. 1$\sigma$ uncertainties are given in parentheses.}
\tablecomments{Table~\ref{table:HST_phot} will be published in its entirety in the machine-readable format. A portion is shown here for guidance regarding its form and content.}
\end{deluxetable*}

\subsubsection{Ground-based Imaging and Photometry}
We obtained follow-up imaging in the near-IR with the MMT and Magellan Infrared Spectrograph (MMIRS; $0\farcs2$ pixels; \citealp{mcleod12}) on the 6.5-m MMT Observatory telescope on Mt.\ Hopkins at the Smithsonian's Fred Lawrence Whipple Observatory. We obtained $JHK_s$ imaging over a two-night run on UT 2021 February 23 and 24, consisting of dithered sequences alternating between the target position in the central regions of M51 and an offset blank-sky field every few minutes to allow for accurate subtraction of the bright near-IR sky background. We reduced the images using a custom pipeline\footnote{Adapted from the MMIRS imaging pipeline developed by K.\ Paterson available here: \url{https://github.com/CIERA-Transients/Imaging_pipelines}} that performs standard dark-current subtraction, flat-fielding, sky background estimation and subtraction, astrometric alignments, and final stacking of the individual exposures. 

We also obtained $J$ and $K_s$-band imaging with the Near-Infrared Imager\footnote{\url{http://www.gemini.edu/instrumentation/niri}} (NIRI) on the 8-m Gemini N Telescope on Maunakea through a Directors Discretionary Time program (PI: J.\ Jencson; PID: GN-2021A-DD-101). The images were obtained withe the f/6 camera ($0\farcs117$ pixels) over multiple nights from 2021 April 1--6, again using sequences of dithered images alternating between the target position and an offset sky position. We also downloaded NIRI $JHK$ images from the Gemini Observatory archive
covering the location of M51-DS1 and taken with in a similar manner the same camera setup on 2005 June 27 (PI: S.\ Smartt; PID: GN-2005A-Q-49). We reduced all the images using \texttt{DRAGONS}\footnote{\url{https://dragons.readthedocs.io/en/v2.1.1/index.html}}, a Python-based platform for reducing Gemini data, and following the procedures for extended sources outlined in the NIRI imaging reduction tutorial\footnote{\url{https://dragons.readthedocs.io/projects/niriimg-drtutorial/en/stable/}}. 

With the large field-of-view of the MMIRS imager ($6.9 \times 6.9$~arcmin), we were able to calibrate the photometric zeropoints using aperture photometry of relatively isolated stars in the MMIRS images with cataloged $JHK_s$-band magnitudes in the Two Micron All Sky Survey (2MASS; \citealp{skrutskie06}). We then derived a model of the ePSF for each image by fitting bright, isolated stars using the \texttt{EPSFBuilder} tool of the \texttt{photutils} package in \texttt{Astropy}. We performed PSF-fitting photometry at the location of M51-DS1 as well as for a set of approximately 60 stars spread across the images with varying degrees of crowding and background galaxy emission. We include a low-order, two-dimensional polynomial in the fit to account for the spatially varying background for each star, taking care to avoid over-fitting the data. We adopt the rms error of the fit residuals, scaled by a factor of the square root of the reduced $\chi^2$ (typically $\gtrsim 1$) for the fit, as the nominal statistical uncertainty per pixel, and multiply by the effective footprint, or number of ``noise pixels'', of the ePSF\footnote{A derivation of this quantity is provided by F.\ Masci here: \url{http://web.ipac.caltech.edu/staff/fmasci/home/mystats/noisepix_specs.pdf}} to obtain an estimate of the statistical uncertainty for each flux measurement. We used the set of 2MASS calibration stars to derive aperture corrections ($\lesssim0.1$\,mag in all three filters) to place the PSF-fitting magnitudes on the scale of the image photometric zeropoints. We adopt the statistical flux uncertainty, summed in quadrature with the rms error of the stars used in estimations of the zeropoint and ePSF aperture correction, as the total uncertainty in our final magnitudes. Owing to the limited number of isolated 2MASS stars, even with the large field-of-view of MMIRS, the zeropoint rms (typically $\approx$0.1\,mag) dominates the error budget. 

The field-of-view of the NIRI images is smaller ($\approx2\arcmin\times2\arcmin$), and there were not enough isolated 2MASS stars in the central regions of M51 to do a direct calibration. We instead cross-calibrated our PSF-photometry of stars in the NIRI image, performed in the same manner as described above, to the set of $\approx$15--20 common stars with the corresponding MMIRS image in the same filter. For our measurements of M51-DS1, we adopt the statistical uncertainty from the PSF-fitting (as above), summed in quadrature with the zeropoint uncertainty (from the standard deviation of the individual stars used in the cross-calibration) as our measurement uncertainty. 


The image quality of all the ground-based near-IR data was between $\approx$0$\farcs5$ and $0\farcs8$ FWHM---a factor of $\approx$2--3 times that of WFC3/IR---generally with better seeing in the redder $H$ and $K$ bands. As noted in Section~\ref{sec:imaging}, M51-DS1 is relatively isolated in the near-IR, with the nearest comparably bright object in the WFC3/IR images at a separation of $\approx$1$\arcsec$. Still, our ground-based photometry may suffer some contamination from nearby unrelated sources, particularly in the $J$ band. The 2005 and 2021 $J$-band fluxes are a bit higher than those from the comparable $F110W$ images, though the discrepancies are $\lesssim1.7\sigma$. The 2021 $H$-band flux is fully consistent with the comparable $F160W$ measurement, and we also expect any contamination to be small or negligble in the $K$-band. Our final photometry of M51-DS1 from the ground-based near-IR images is provided in Table~\ref{table:GB_nearIR}.

\begin{deluxetable*}{ccccc}
\tablecaption{Ground-based Near-IR Photometry \label{table:GB_nearIR}}
\tablehead{\colhead{UT Date} & \colhead{MJD} & \colhead{Tel./Inst.} & \colhead{Band} & \colhead{App.\ Magnitude\tablenotemark{a}} \\ 
\colhead{} & \colhead{} & \colhead{} & \colhead{} & \colhead{(mag)} }
\startdata
2005 Jun 27.34 & 53548.34 & Gemini/NIRI & $J$   & 19.61 (0.14) \\
2005 Jun 27.36 & 53548.36 & Gemini/NIRI & $H$   & 19.14 (0.16) \\
2005 Jun 27.38 & 53548.38 & Gemini/NIRI & $K$   & 18.54 (0.21) \\
2021 Feb 23.36 & 59268.36 & MMT/MMIRS     & $J$   & 20.04 (0.11) \\
2021 Feb 23.49 & 59268.49 & MMT/MMIRS     & $H$   & 19.26 (0.09) \\
2021 Feb 24.40 & 59269.40 & MMT/MMIRS     & $K_s$ & 18.73 (0.13) \\
2021 Apr 01.45 & 59305.45 & Gemini/NIRI & $J$   & 20.04 (0.08) \\
2021 Apr 02.41 & 59306.41 & Gemini/NIRI & $K_s$ & 18.83 (0.12) \\
\enddata
\tablenotetext{a}{Vega magnitudes on the 2MASS system. 1$\sigma$ uncertainties are given in parentheses.}
\end{deluxetable*}

\section{Analysis of M51-DS1} \label{sec:analysis}

\subsection{Extinction, Distance, Host Galaxy Environment, and Metallicity}\label{sec:M51_env}
M51-DS1 is located at an R.A. and decl.\ of
$13^{\mathrm{h}}29^{\mathrm{m}}56\fs16, +47\degr11\arcmin47\farcs8$ (J2000.0), near the inner regions of the nearby, star-forming galaxy M51. Throughout this work, we assume a distance modulus for M51 from \citet{mcquinn16} of $m-M = 29.67 \pm 0.02$ (statistical) $\pm 0.07$ (systematic; \citealp{rizzi07})~mag ($D = 8.58$\,Mpc) based on the luminosity of the tip of the red giant branch (TRGB) method, and that the systematic uncertainties associated with calibrating this method dominate over the statistical measurement uncertainties. We adopt the value from the NASA/IPAC Infrared Science Archive (IRSA) for the Galactic extinction toward the position of M51-DS1 of $E(B-V) = 0.03$~mag, based on the \citet{schlafly11} recalibration of the \citet{schlegel98} dust maps, and assuming a standard \citep{fitzpatrick99} reddening law with $R_V = 3.1$.

The star is $\approx$36$^{\arcsec}$ from the galaxy nucleus along the outer edge of a densely populated spiral arm and prominent dust lane (see Figure~\ref{fig:hst_im}). We thus expect the region may also be subject to substantial and spatially-variable extinction from the host. \citet{croxall15} presented spectroscopic, gas-phase chemical abundances for 29 \ion{H}{2} regions in M51 including extinction estimates derived from Balmer-line ratios as part of the CHemical Abundances of Spirals (CHAOS) program. The two nearest \ion{H}{2} regions to M51-DS1 from this study are labelled in Figure~\ref{fig:hst_im}. One is centered on the dark dust lane with $E(B-V) = 1.0$\,mag, while the other is along the inside edge of the spiral arm in a region that appears visually similar to that of M51-DS1 with $E(B-V) = 0.4$\,mag. \citet{wei21} recently derived extinction maps of the M51 system from UV/optical and IR photometry, from which we again find M51-DS1 to be associated with a region of higher extinction with $A_V \approx 1.0$--$1.2$\,mag, or $E(B-V) \approx 0.3$--0.4\,mag for $R_V = 3.1$. These estimates are based on spatial averages in seeing-limited spectra ($\approx$1$\arcsec$; for CHAOS) or lower resolution space-based imaging (e.g., $\approx$3$\arcsec$ with \textit{Spitzer} in \citealp{wei21}) and while they are useful guideposts, the appearance of the dust lanes in the higher resolution \textit{HST} images indicates that there is likely substantial variation in the foreground host extinction on smaller spatial scales. In our analysis below, we therefore consider the implications of various amounts of additional reddening due to dust in the interstellar medium of M51, ranging from $E(B-V)$ = 0 mag to 1.6 mag.

From their set of \ion{H}{2} regions, \citet{croxall15} found O/H abundances in the central regions of M51 to be approximately solar or somewhat supersolar (assuming solar abundances of $12+\log(\mathrm{O/H})_{\odot} = 8.69$--8.78, \citealp{asplund09,ayres13}). Similarly, the metallicity map derived by \citet{wei21} suggests somewhat supersolar metallicities of $\log(Z/Z_{\odot}) \approx 0.0$--0.3\,dex in the vicinity of M51-DS1 along the inner spiral arm. 

\subsection{Optical/near-IR Light Curves and Colors}\label{sec:lcs}
The light curves of M51-DS1, dating back to 1995, are shown in Figure~\ref{fig:hst_lcs}. The magnitudes displayed and discussed here have been corrected only for Galactic extinction. The star was detected dozens of times at $F814W$ between 1995 and 2017, exhibiting some variability at the level of $\approx$0.7\,mag between $F814W = 22.74 \pm 0.20$ and $22.08 \pm 0.06$\,mag ($M_{F814W} = -6.9$ to $-7.6$\,mag). During 2016--2017, the evolution of the source was tracked at high cadence in both $F814W$ and $F606W$ with ACS/WFC (see inset in Figure~\ref{fig:hst_lcs}). The source displays a smooth, slow rise in both bands at a relatively constant, red color between $F606W - F814W = 1.6$ and 1.9\,mag, where the spread in the color measurements is comparable to their uncertainties. 

\begin{figure*}
\centering
\includegraphics[width=\textwidth]{M51DS1_HST_lcs.pdf}
\caption{\label{fig:hst_lcs}
The light curves of M51-DS1 in $F606W$, $F814W$, and $F110W$ from archival and new \textit{HST} images (filled symbols) going back more than 26\,yr to 1995. The region enclosed in the dash-lined box of high-cadence measurements from $F814W$ and $F606W$ imaging taken between 2016--2017 is shown in more detail in the inset zoom-in panel. Ground-based near-IR measurements from 2005 and new 2021 images are shown as open symbols. All photometry has been corrected for Galactic extinction to M51. 
}
\end{figure*}

Subsequently, the star underwent a dramatic fading in the 2019 ACS/WFC $F814W$ observations, dropping to at least $F814W > 24.2$ ($>-5.5$\,mag) in our \texttt{DOLPHOT} measurement, $>2.2$\,mag below the brightest prior detection. Though \texttt{DOLPHOT} nominally obtained a detection (Section~\ref{sec:dolphot}), the source looks to have nearly completely disappeared (Fig~\ref{fig:hst_lcs}; bottom, center panel) at $F814W$, and given the high degree of crowding at this location, we treat this measurement as an upper limit. Then, in our 2021 follow-up observations approximately two years later, the source has rebrightened partially to $F814W = 22.97 \pm 0.09$\,mag ($-6.6$\,mag; similar to the prior $F814W$ observed in 2008), and displays a slightly redder color at $F606W - F814W = 1.97 \pm 0.16$\,mag compared to the 2016--2019 measurements.

As expected given its red optical color, the source was brighter in the near-IR at $F110W = 20.05 \pm 0.17$\,mag ($M_{F110W} = -9.6$\,mag) in the 2012 WFC3/IR image. In our recent 2021 \textit{HST} follow-up imaging, around the time of the partial recovery seen in the optical, the star is also fainter at $F110W = 20.44 \pm 0.18$ ($M_{F110W} = -9.2$\,mag) with a red near-IR color of $F110W - F160W =  1.0 \pm 0.2$\,mag. The ground-based imaging with NIRI and MMIRS tells a similar story. In the 2005 NIRI imaging, the star is detected at $J = 19.58 \pm 0.14$ ($M_J = -10.1$\,mag), with $J-K = 1.05 \pm 0.25$\,mag. Subsequently, in our 2021 ground-based follow-up imaging, the star is again fainter in the near-IR filters with NIRI (MMIRS) compared to the pre-2019 levels at $J = 20.01 \pm 0.08$\,mag ($20.01 \pm 0.11$\,mag) and $J - K_s = 1.20 \pm 0.14$\,mag ($1.29 \pm 0.17$\,mag). 

\subsubsection{Near-IR Photometric Classification and Bolometric Luminosity}\label{sec:nearIR_Lbol}
At $M_K = -11.1$\,mag (Galactic extinction correction only), the star is well above the TRGB, and moreover, is brighter than any asymptotic giant branch (AGB) stars identified in nearby galaxies including the Large and Small Magellanic Clouds (L/SMC), M31, and M33 \citep[see, e.g.,][and references therein]{cioni06,boyer11,massey21}. Its near-IR color of $J-K_s \approx 1.0$--1.2\,mag is also consistent with the range often used to identify RSGs, though it may be somewhat bluer given the likelihood of significant foreground host extinction (see Section~\ref{sec:M51_env}). The $K$ band is particularly useful as a luminosity indicator for RSGs, both because the effects of extinction are reduced compared to the optical and bluer near-IR bands, and because the bolometric correction, ${\rm BC}_K$, is found empirically by \citet{davies18a} to be constant across early- to late-M spectral types for cool supergiants in Milky Way and LMC star clusters. Thus, assuming an M-type spectrum ($T_{\rm eff} \lesssim 3700$\,K) for M51-DS1 and adopting their value of ${\rm BC}_K = 3.0$\,mag, we obtain bolometric luminosities in the range of $\log (L/L_{\odot}) \approx 5.15$--5.37 depending on the amount of extinction assumed in excess of the Milky Way foreground from $E(B-V) = 0.0$ up to 1.6\,mag (the highest value found to produce good fits in our modelling of the SED in Section~\ref{sec:SED_models}). The value for zero extinction, $\log (L/L_{\odot}) = 5.15$, can likely be viewed as a robust lower limit, as any foreground extinction or an earlier intrinsic spectral type ($<$M0 with ${\rm BC}_K < 3.0$), will both increase the inferred bolometric luminosity of the star. We discuss the possible location of the star in the Hertzsprung--Russell Diagram (HRD) and inferred evolutionary state in comparison with the results of our SED fitting in Section~\ref{sec:HRD}. 

\subsection{Color-magnitude diagrams}\label{sec:CMDs}
In Figure~\ref{fig:HST_CMDs}, we present the F606W--F814W color-magnitude diagram (CMDs) derived from our \texttt{DOLPHOT} photometry on the 2021 ACS/WFC images to place M51-DS1 in the context of the nearby stellar population. Correcting only for Milky Way extinction, M51-DS1 appears as the most luminous red star ($F606W - F814W \gtrsim 2$\,mag) at $F814W$ within a projected distance of 50 ($1\farcs2$) or 100\,pc ($2\farcs4$) at its 2016--2017 level. Notably as well, it is one of only two stars within 50\,pc that appears so red---the other is indicated by the magenta circle in both of the figure panels. Compared to the single-star, non-rotating, solar-metallcity stellar tracks from the Mesa Isochrones and Stellar Tracks models (MIST; \citealp{choi16,choi17}), its location in the CMD would correspond to an RSG with an initial mass of $\approx 12$--15~$M_{\odot}$ (2016--2017), or $\approx 8$--10~$M_{\odot}$ (2021), and an age below $\approx$15 or 30~Myr, respectively. As noted in Section~\ref{sec:M51_env}, however, there is likely to be significant, and spatially variable foreground extinction in the region. For values of $E(B-V)$ between $0.4$ and $1.0$\,mag---in the range inferred from previous analysis of nearby \ion{H}{2} regions by \citep{croxall15}---the position of the star would correspond to higher masses up to $\approx$20--25~$M_\odot$, and younger ages $\lesssim$7--10~Myr. Our analysis of the full SED of M51-DS1 in Section~\ref{sec:SED_models} suggests even higher extinction values of $E(B-V) = 0.8$--1.6\,mag, possibly pointing to significant circumstellar extinction, corresponding to intial masses as high as $\approx$25--40~$M_\odot$ and implying an age as young as $\approx$5~Myr in the CMD. 

\begin{figure*}
\centering
\includegraphics[width=0.48\textwidth]{CMD_image.pdf}
\includegraphics[width=0.51\textwidth]{NGC5194-DS1_ACS_WFC_F606W_F814W_F814W_20210428_evenwarmstart_CMD.pdf}
\caption{\label{fig:HST_CMDs}
\textit{Left:} Immediate environment of M51-DS1 (red ``x''-symbol) in the 2005 ACS/WFC images (color-composite of $F435W$, $F555W$, and $F814W$, as in Figure~\ref{fig:hst_im}) showing a mix of red and blue stars and dark regions that are likely obscured by dust. Radii corresponding to 20 ($0\farcs48$) and 50 pc ($1\farcs2$) projected distances from M51-DS1 are indicated.
Sources indicated by multicolor symbols correspond to those objects marked in the CMD in the right panel. \textit{Right:} The CMD of the region around M51-DS1 from the 2021 ACS/WFC $F606W$ and $F814W$ images. Stars within 20, 50, and 100 pc (projected; $0\farcs48$, $1\farcs2$, and $2\farcs4$, respectively) from the location of M51-DS1 are shown as orange, blue, and gray points, respectively. Additional multicolor markers correspond to those objects indicated in the left panel. We also show single-star, solar metallicity, non-rotating stellar evolutionary tracks from MIST as the black and gray dashed curves for a range of masses between 8 and 40~$M_\odot$. Measurements for M51-DS1 are shown as the large stars with varying degrees of foreground host extinction, $E(B-V)$, indicated by different colors as labeled.
We include several M51-DS1 epochs of $F606W$ and $F814W$ imaging including 2016--2017 and 2021. The nominal 2019 $F814W$ \texttt{DOLPHOT} measurement (uncertainty) is indicated as the red horizontal line (shaded bar), though we consider this as an upper limit (see main text) as indicated by the downward arrow. 
}
\end{figure*}

As is evident from the tricolor image shown in Figure~\ref{fig:HST_CMDs}, the immediate region surrounding M51-DS1 is populated by a range of blue-, yellow- and red-appearing stars. This is reflected as a notable spread in the $\mathrm{F606W} - \mathrm{F814W}$ color between $\approx -0.3$ and $1.6$\,mag for the bulk of the stars in the CMD within the 20- and 50-pc radii. If one assumes a coeval population of stars undergoing isolated evolution, this could be attributed to effects of foreground dust, possibly implying variable extinction in the range $E(B-V) \approx 0.0$--$2.0$\,mag. Interpretation of the CMD as a whole, however, is complicated by additional factors on top of the uncertain extinction, including the possibility of multiple, overlapping stellar populations of different ages, and the effects of binary interaction---known to be increasingly common among more massive stars \citep[e.g.,][]{sana12,moe17}---on the location of a given source in the diagram. 

We look briefly at a few individual sources of particular note. Within 50~pc, the apparently brightest object at $F814W$ in the CMD, along with two somewhat fainter objects, are associated with what appears to be a partially resolved star cluster (indicated with cyan squares) and may not represent secure detections of individual stars. After this, there are two bright, yellowish stars (indicated with yellow diamonds), at $M_{\mathrm{F814W}} \approx -7.4$\,mag and $\mathrm{F606W} - \mathrm{F814W} = 0.7$\,mag, which, in the absence of additional foreground extinction, would correspond to stars of at least $\approx$12--15~$M_{\odot}$ evolving across the Hertzsprung gap. Within a smaller radius of 20~pc ($0\farcs48$), the most luminous objects after M51-DS1 top out at $M_{\mathrm{F814W}} \approx -6.5$\,mag, including a relatively blue (presumably unreddened) source at $\mathrm{F606W} - \mathrm{F814W} = -0.3$\,mag (indicated in both the image and CMD by a blue circle). This would correspond to a $\approx$30~$M_{\odot}$ star at the end of its main-sequence lifetime and implying an age $\lesssim$6~Myr. This is consistent with the higher mass and younger age inferred for M51-DS1 above when one allows for significant foreground and/or circumstellar extinction ($E[B-V] \gtrsim 0.8$\,mag) and assuming it is part of the same, coeval population. 

\subsubsection{Comparison to Luminous, Red, Large-amplitude Variable Stars}
Regardless of the magnitude of the foreground extinction to M51-DS1, the track of possible locations it occupies in the CMD in Figure~\ref{fig:HST_CMDs} essentially aligns with the red-end points of the MIST evolutionary tracks for stars of initial masses $\gtrsim 15~M_{\odot}$. RSGs in this region are expected to be susceptible to pulsational instabilities driven by partial ionization of hydrogen in their envelopes (e.g., \citealp{li94,yoon10}, and see region of ``supergiant instabilities'' in Figures~13 and 16 of C18). 

To compare M51-DS1 to other such luminous, red variable stars, we selected stars from the catalog of M51 luminous variable star light curves presented in C18. Of the more than $\approx$72000 objects in their catalog (including only the Milky Way extinction corrections already applied), we pulled out the 90 objects having $M_{\rm F814W} \leq -7.0$\,mag at any epoch (1995, 2005, and 2016--2017 data), average colors (including epochs within 0.5\,mag of the $F814W$ peak) of $F606W - F814W \geq 1.5$\,mag, and maximum $F814W$ amplitudes $\Delta F814W \geq 1.0$\,mag. We examined both the images and light curves of these objects by eye to remove false positives such as foreground proper motion stars and stars with large amplitude variability in the light curves that is not apparent in visual inspection of the images, likely indicating erroneous photometry. We also cull the sample by rejecting photometric points with large uncertainties ($\sigma_m > 0.2$\,mag) and stars where the color evolution is dominated by light curve scatter, usually in the fainter $F606W$ measurements. We thus select a clean sample of 8 luminous, red, large-amplitude variables, and show their evolution in the CMD in Figure~\ref{fig:LRVs_CMD}. 

\begin{figure*}
\centering
\includegraphics[width=0.49\textwidth]{M51-DS1_Betelgeuse_LRVs_CMD.pdf}
\includegraphics[width=0.49\textwidth]{M51-DS1_Betelgeuse_LRVs_DCMD.pdf}
\caption{\label{fig:LRVs_CMD}
CMD of luminous ($M_{F814W} \leq -7.0$\,mag), red ($F606W - F814W \geq 1.5$), high-amplitude ($\Delta F814W \geq 1.0$\,mag) variable stars in M51 from the catalog of C18 (small circles connected by dashed lines) compared to the 2016--2017, 2019, and 2021 data of M51-DS1 (red stars). The plotting colors of the catalog variable points correspond to their photometric color evolution during fading fading events: to the blue (blue), approximately constant (yellow), and to the red (dark orange). We also show the $M_I, V-I$ evolution of Betelgeuse during its 2019 dimmming event from AAVSO, and the location of the progenitor of NGC6946-BH1 from \citep{adams17b}. As in Figure~\ref{fig:HST_CMDs}, we show single star tracks of various masses from MIST.
}
\end{figure*}

These 8 stars can then be separated according to their color evolution: Those that evolve to the red as they fade, those that remain at a relatively constant color, and those that evolve to the blue (shown as orange, yellow, and blue circles in Figure~\ref{fig:LRVs_CMD}, respectively). The distinctions between these groups are seen more clearly in the differential CMD shown in the right panel of the figure, in which the average $F814W$ magnitudes and $F606W - F814W$ colors (using points within 0.5\,mag of the $F814W$ peak) have been subtracted out for each source. Interestingly, the ``yellow''-class objects appear near the end of the RSG evolutionary tracks at the Hayashi limit (though we note that only foreground Milky Way extinction corrections have been applied), possibly indicating that their near-vertical evolution in the CMD may indeed be attributed to large-amplitude, radial pulsations that are theoretically expected for such stars. Between the 2016--2017 and 2021 data, the evolution of M51-DS1 appears most similar to this group of stars. Unfortunately, during the $F814W$ minimum in 2019, we lack any information on its color evolution. 

The red-evolving objects notably become much redder than the end of the MIST tracks at $F606W - F814W > 2.5$\,mag. While some of these stars may suffer additional host reddening not accounted for here, their color evolution is nonetheless striking. A natural explanation of this would be the formation of dust, possibly during an episode of enhanced mass loss from the star, that temporarily obscures and reddens the light from the photosphere. Alternatively, an increase in the opacity of the strong molecular absorption bands, in particular of TiO, that dominate RSG spectra can produce substantial swings in the optical colors without a correspondingly large decrease in temperature. The $F814W$ and $F606W$ filters of interest here notably both contain strong TiO bands. \citet{davies21} recently explored the effect of winds on RSG spectra, finding that the addition of winds increases the molecular absorption, effectively shifting the stars to later spectral types at a fixed photospheric temperature and producing excess mid-IR emission without needing to invoke dust. 

These effects have been discussed extensively regarding the recent, historic photometric minimum of the RSG Betelgeuse, dubbed the ``Great Dimming'', in which the star faded by more than a magnitude in the $V$ band and then recovered over the course of $\approx$5~months between the end of 2019 and early 2020 \citep{guinan19,guinan20,dupree20}. Multi-epoch, high-resolution imaging revealed the appearance of an optically-dark patch over part of the stellar surface during the dimming and exclude a transiting dust cloud \citep{montarges21}. UV observations of chromospheric variability also indicate the formation of a dense, outflowing structure---possibly enhanced by a timing coincidence with outward motion of the pulsational cycle---just before the event \citep{dupree20}. The most likely explanation for the dimming, then, was an increase in the opacity of this ouflowing material related to the formation of molecules (see, e.g., evidence from TiO-band monitoring and tomography in \citealp{harper20,kravchenko21}, and discussion in \citealp{davies21}) and possibly also large-grain dust \citep[e.g.,][]{levesque20}.
As shown in Figure~\ref{fig:LRVs_CMD}, the evolution of Betelgeuse in $V$ and $I$ from the American Association of Variable Star Observers\footnote{\url{https://www.aavso.org/}} (AAVSO) as it recovers from near the photometric minimum of the 2019--2020 dimming is remarkably similar to the ``red''-class variables, though smaller in amplitude. We suggest that these objects may represent similar or even more extreme instances of episodic, enhanced mass-loss in RSGs, though we reserve a full analysis for future investigations. This may signal that many RSGs experience Betelgeuse-like dimming events during this phase of their evolution.

We briefly comment on the two blue-evolving sources shown in Figure~\ref{fig:LRVs_CMD} as well. Given their red colors, we caution that the apparent evolution of these stars may be affected by contamination from nearby objects during their minima at $F606W$. Still, their evolution is reminiscent of that typical of bluer massive stars, namely the characteristic S Doradus variations between a hot (visibly fainter), quiescent state and a cool (visibly brighter) outburst state \citep[e.g.,][]{humphreys94}.  Theoretically, massive stars in the bluer portion of the CMD are expected to be susceptible to instabilities arising in radiation-dominated envelopes that approach the local Eddington limit (\citealp[e.g.,][]{paxton13,jiang15,owocki15} and see the ``radiation-dominated instabilities'' region in Figures~13 and 16 of C18). We suggest that the ``blue''-evolving variables could be intrinsically bluer stars, possibly LBVs or yellow supergiants experiencing LBV-like episodes \citep[as in][]{smith04a}, that are now reddened by the dust in the foreground ISM or formed in their own circumstellar environments during prior eruptive mass-loss events. A detailed examination of these objects is outside the scope of the present work. 

Finally, having laid out the modes of variability observed (and expected) for luminous, red-appearing stars, we return to the variability of M51-DS1. As noted above, where color information is available for the 2016--2017 and 2021 epochs, M51-DS1 tracks with the mostly vertical evolution of the yellow-colored variables in Figure~\ref{fig:LRVs_CMD}, consistent with expectations for pulsational instabilities occuring in massive ($\gtrsim 15 M_{\odot}$) RSGs, though large-grain circumstellar dust may also produce extinction that is relatively gray \citep[see, e.g.,][]{scicluna15,massey06,haubois19}. The long period of brightening seen between 2016--2017 may correspond at most to one half of a full pulsation cycle, requiring that any periodicity be $\gtrsim$2\,yr in duration. Given the star's high luminosity ($M_K = -11.1$\,mag), this is consistent with expectations from the RSG period--luminosity relation observed in the Milky Way, L/SMC, M33 and M31 \citep{kiss06,yang12,soraisam18}. The maximum amplitude of this variability (excluding the 2019 minimum), is $\approx$0.9\,mag, only slightly larger than pre-2019 variability observed at the level of $\Delta F814W \approx 0.7$\,mag. At an exceptional $\Delta F814W \gtrsim 2.1$\,mag, it is unlikely, though not strongly excluded, that the 2019 minimum is part of the inferred (semi)regular pulsation cycle ($\lesssim10$\% chance of finding the source in the bottom portion of a pulsation only once in the 6 light curve samples over 26 years). We therefore suggest that the 2019 dimming event of M51-DS1 represents a relatively rare occurrence, possibly an exceptional mass-loss event and more extreme cousin of ``Great Dimming'' of Betelgeuse. Though we lack the necessary color information during the dimming event to confirm this for M51-DS1, a statistical analysis of the ensemble of massive variables in M51 with the wealth of archival and newly obtained \textit{HST} imaging in a future study will establish the frequency and range of properties (i.e., mass-loss rates, durations, dust masses) of episodic mass loss in RSGs. 

\subsection{SED Analysis} \label{sec:SEDs}
Figure~\ref{fig:SEDs} shows the multi-epoch SEDs of M51-DS1, constructed from the available \textit{HST} and ground-based photometry. The photometric magnitudes were converted to band-luminosities, $\lambda L_{\lambda}$, using the zero point fluxes and effective wavelengths available with the \texttt{pysynphot} package \citep{pysynphot} for \textit{HST} instrument and filter setups, and those compiled by the Spanish Virtual Observatory (SVO) Filter Profile Service\footnote{Documentation for the SVO Filter Profile Service is available at \url{http://ivoa.net/documents/Notes/SVOFPSDAL/index.html} and \url{http://ivoa.net/documents/Notes/SVOFPSDAL/index.html}} for the ground-based observations \citep{svo1,svo2}. 

\begin{figure*}
\centering
\includegraphics[width=0.32\textwidth]{NGC5194_DS1_SEDs.pdf}
\includegraphics[width=0.32\textwidth]{NGC5194_DS1_SED_GRAMS_2017_2012_2005.pdf}
\includegraphics[width=0.32\textwidth]{NGC5194_DS1_SED_GRAMS_2021.pdf}
\caption{\label{fig:SEDs}
\textit{Left}: Multi-epoch SEDs of M51-DS1 constructed from visible and near-IR photometry. Space-based measurements from \textit{HST} are shown as filled symbols, while ground-based near-IR measurements are open symbols. Upper limits are indicated with downward arrows. \textit{Center:} The pre-2019 SED constructed from the 2017 ACS/WFC, 2012 WFC3/UVIS and IR, and 2005 NIRI photometry (black diamonds) is shown along with GRAMS RSG models that provide good fits ($\chi^2 \leq \chi^2_{\mathrm{min}} + \Delta \chi^2_{90}$; see text) to the data. The best fitting model, described by the parameters listed in the upper-left corner of the panel, is shown as the orange, thick-dashed curve, with the corresponding synthetic photometry points in the available bands shown as red squares. The input (unreddened) stellar photosphere to the model is shown as the blue dashed curve. \textit{Right:} Same as the center panel, but for the 2021 SED constructed from our ACS/WFC, WFC3/IR and MMIRS photometry. 
}
\end{figure*}

Prior to its dramatic dimming in 2019, the source was characterized by a very red SED peaking in the near-IR between 1--2\,$\mu$m at $\lambda L_{\lambda} \approx 10^{5}~L_{\odot}$. The variability seen in the light curves (Section~\ref{sec:lcs}) is also reflected in the SEDs constructed from the 2005, 2012, and 2016--2017 data, though the overall color and shape of the SED remains largely constant. In 2019, the source faded by $\Delta \lambda L_{\lambda} \gtrsim 2.6\times10^4~L_{\odot}$ at $F814W$, a factor of $\gtrsim 7$ or $\gtrsim 85$\% of the flux in that band. Without multi-band observations in 2019, we are unable to constrain the shape of the SED or total drop in bolometric luminosity at that epoch. During our follow-up 2021 observations, the source has partially recovered in flux at $F814W$. The 2021 SED also peaks in the near-IR at $\lambda L_{\lambda} \approx 8.2\times10^4~L_{\odot}$ in $F160W$. While overall similar in color and shape, the source is dimmer across all the available optical and near-IR bands compared to the pre-2019 SEDs.

\subsubsection{SED Modelling}\label{sec:SED_models}
To estimate physical parameters of the star, we attempted to fit the SEDs with the Grid of Red supergiant and Asymptotic Giant Branch ModelS (GRAMS; \citealp{sargent11,srinivasan11}). This suite of radiative transfer models consists of a base grid of 1225 spectra from spherically symmetric shells of varying amounts of silicate dust \citep[][appropriate for RSGs]{ossenkopf92} around stars of constant mass-loss rates computed using the dust radiative transfer code \texttt{2-Dust} \citep{ueta03}. The published grid uses input PHOENIX model photospheres \citep{kucinskas05,kucinskas06} for 1\,$M_{\odot}$ stars (model spectra can be scaled for more luminous and massive, i.e., supergiant, stars) with effective temperatures, $T_{\mathrm{eff}}$, between 2100 and 4700\,K, and at a fixed sub-solar metallicity of $\log(Z/Z_{\odot}) = -0.5$  and a fixed surface gravity of $\log g = -0.5$. The amount of circumstellar dust is characterized in terms of the optical depth at 1~$\mu$m, $\tau_1$, from which a dust mass-loss rate, $\dot M_{\mathrm{d}}$ can be inferred assuming a wind speed of $v_w = 10$\,km\,s$^{-1}$. An additional parameter in the GRAMS grid is the inner radius of the dust shell, $R_{\mathrm{in}}$. We found that our results were largely insensitive to this parameter and we chose to fix it at $R_{\mathrm{in}} = 11.0 R_*$, where $R_*$ is the stellar radius, to reduce the number of free parameters.

To model the source prior to its near-disappearance in 2019, we combined \textit{HST} measurements from the 2017 ACS/WFC observations in $F606W$ and $F814W$ (MJD 58011) and the 2012 WFC3/UVIS observations in $F673N$ and $F689M$ (MJD 56020--56027) and IR observations in $F110W$ and $F128N$ (MJD 56174), together with the 2005 ground-based NIRI measurements in $J$, $H$ and $K$ (MJD 53548). Given the observed variability of the source during this time period, we experimented with fitting different combinations of the available pre-2019 data, but found that these epochs together gave the best results. We also included additional reddening (\citealp{fitzpatrick99} law with $R_V =3.1$) for a range of $E(B-V)$ values between 0.0 and 3.0\,mag to account for the likely significant and uncertain host extinction (see Section~\ref{sec:M51_env}). We then fit each (reddened) model spectrum in the expanded grid (99225 models in total) by computing $\chi^2$ values between the data points and synthetic photometry on the spectra, weighted by data point uncertainties and allowing for an overall scaling factor of the flux as a free parameter. In our analysis, we will consider ``good'' models as those with $\chi^2 \leq \chi^2_{\mathrm{min}} + \Delta \chi^2_{90}$, where $\chi^2_{\mathrm{min}}$ is minimum value of $\chi^2$ found for all the models, and $\Delta \chi^2_{90}$ is the 90\%-tile of the $\chi^2$ distribution with degrees of freedom $\nu$. 

The results of this fitting procedure for the pre-2019 SED are shown in the center panel of Figure~\ref{fig:SEDs}. The best-fitting model ($\chi^2/\nu = 3.12$) has $T_{\mathrm{eff}} = 4300$, $E(B-V) = 1.28$ and $\log(L/L_{\odot}) = 5.60$. With $\nu = 5$ degrees of freedom, $\Delta \chi^2 = 9.24$, and we find the set of 257 good models with $T_{\mathrm{eff}} = 3700$--4700~K, $E(B-V) = 0.56$--1.56\,mag and $\log(L/L_{\odot}) = 5.35$--5.75. Unsurprisingly, $T_{\mathrm{eff}}$ and $E(B-V)$ are strongly positively correlated for the set of good models, and we note that the range of $T_{\mathrm{eff}}$ they span extends to the maximum value in the grid. The $\chi^2$ distribution is well behaved around the minimum at $T_{\mathrm{eff}} = 4300$~K in the range we tested, although highly reddened, warmer models could, in principle, also provide acceptable fits to the SED.  These results point to a very luminous, cool ($\approx$K or early-M spectral type) supergiant star, near the empirical Humphreys-Davidson limit ($\log[L/L_{\odot}] = 5.5$--5.8; \citealp{humphreys79,davies18b}). In comparison to single-star evolutionary tracks from the Mesa Isochrones and Stellar Tracks models (MIST; \citealp{choi16,choi17}), this would correspond to a star with an initial mass of $\approx$24--40~$M_{\odot}$  with terminal age $\lesssim 8$~Myr (see Figure~\ref{fig:HRD} and further discussion in Section~\ref{sec:HRD}). 

We fit the post-fading SED, constructed from the 2021 ACS/WFC ($F606W$, $F814W$) WFC3/IR ($F110W$, $F160W$), and MMIRS ($J$, $H$, $K_s$) imaging, using the same procedure and show the results in the rightmost panel of Figure~\ref{fig:SEDs}. In this case, the best-fitting model ($\chi^2/\nu = 4.78$) has $T_{\mathrm{eff}} = 4300$, $E(B-V) = 1.64$ and $\log(L/L_{\odot}) = 5.52$. The set of 125 good models ($\nu = 3$; $\Delta \chi^2_{90} = 6.25$) now has $T_{\mathrm{eff}} = 3900$ -- 4700~K, $E(B-V) = 1.0$ -- 1.9\,mag and $\log(L/L_{\odot}) = 5.31$--5.65. Overall, results for the 2021 SED suggest a slightly dimmer star than the pre-fading data, though there is substantial overlap in the allowed parameter ranges. 

While the GRAMS grid includes models with high values of circumstellar extinction up to $\tau_1 \approx 60$, our fitting results for both the pre- and post-fading data favor $\tau_1 \lesssim 0.9$ and corresponding to dust mass-loss rates below $\dot M_{\mathrm{d}} \lesssim 10^{-7}~M_{\odot}$~yr$^{-1}$. This is in accord with modern estimates of mass-loss rates for normal RSG winds in this luminosity range \citep{beasor20}. At the same time, the models favor relatively high values of foreground host (i.e., ISM-like) extinction. Importantly, longer wavelength data ($\gtrsim 3$--5\,$\mu$m and especially at the 10\,$\mu$m silicate feature) that are sensitive to emission from warm dust would be required to directly constrain the amount of circumstellar material (CSM), so our ability to distinguish circum- and interstellar extinction should be viewed cautiously. Notably too, the modelling results for the 2021 data favor higher values of the foreground extinction to account for the redder colors of the source compared to the pre-fading source. An increase in the amount of intervening ISM material is presumably unphysical, indicating that other effects not captured by the models may be at play in the observed SED evolution. These may include the effects of increased molecular opacity during episodes of enhanced mass loss \citep[see, e.g.,][]{davies21}, changes in composition or the grain-size distribution of circumstellar dust, and/or the effects of non-spherical geometry in stellar surface (e.g., large convective bubbles and cold spots) or the CSM.  

Lastly, the input PHOENIX model photospheres of the published GRAMS grid are at a metallicity of $\log(Z/Z_{\odot}) = -0.5$, originally chosen to be similar to that of the LMC, while the local environment of M51-DS1 indicates that models at higher, somewhat super-solar, metallicities ($\log[Z/Z_{\odot}] = 0.0$--0.3; see Section~\ref{sec:M51_env}) would be more appropriate. Though metallicity will affect the photospheric spectrum in the optitcal/near-IR, we do not expect this to be significant compared to the effects of of an external dusty wind or the possible effects of a complex geometry mentioned above. We examined this by also fitting the pre-fading SED with bare PHOENIX photospheres at low, solar, and high metallcities ($\log[Z/Z_{\odot}] = -0.5, 0.0, +0.5$), again including the amount foreground extinction $E(B-V)$ as a free parameter. In general we find similar results regardless of metallicity. The best fitting models have $T_{\mathrm{eff}} = 4300$ at low and solar metallicities and $T_{\mathrm{eff}} = 4500$ at high metallicity, all well within the preferred parameter ranges found for the GRAMS models. 

\subsection{Location in the H--R Diagram and Evolutionary State}\label{sec:HRD}
We examine possible locations M51-DS1 in the Hertzsprung--Russell diagram (HRD; Figure~\ref{fig:HRD}) and discuss implications for the inferred initial mass and evolutionary state of the star. Specifically, we directly compare the results of applying a simple $K$-band bolometric correction \citep[][assuming an M-type RSG; see Section~\ref{sec:nearIR_Lbol}]{davies18a} to those of our SED model fitting (Section~\ref{sec:SED_models}). Overall, the luminosity of the star inferred from applying $BC_K$ is lower than that from the SED fitting, falling in the range $\log (L/L_{\odot}) \approx 5.15$--5.37---notably above the range inferred for the growing collection of SN progenitors and comparable to that inferred for the failed SN candidate NGC\,6946-BH1---depending on the value assumed for the foreground extinction between $E(B-V) = 0.0$ and 1.6\,mag. The SED modelling prefers somewhat warmer ($T_{\mathrm{eff}} = 3700$--4700~K) photospheric models at relatively high extinction ($E(B-V) = 0.56$--1.56\,mag) and luminosities in the range $\log(L/L_{\odot}) = 5.35$--5.75. Given the constraints on ISM environment discussed in Section~\ref{sec:M51_env}, the high $E(B-V)$ values $\gtrsim 1$\,mag would likely point to significant circumstellar extinction.

\begin{figure}
\centering
\includegraphics[width=0.5\textwidth]{M51-DS1_HRD.pdf}
\caption{\label{fig:HRD}
HRD showing possible locations of M51-DS1. The result for best fitting model for the pre-2019 SED is shown as the star symbol, while the range of good models are represented by the small circles. The color of each point corresponds to the value of the extinction (in excess of the Galactic foreground), $E(B-V)$, for that model as indicated by the color bar. The ``x''-symbols show the locations of M51-DS1 assuming early, mid, and late M-type spectra and applying the $K$-band bolometric corrections of \citet{davies18a} to the 2005 NIRI measurement, and with the same mapping of color to $E(B-V)$. We show stellar evolutionary tracks from MIST (non-rotating, solar metallicity) for a set of massive stars in the range $M = 8$--40\,$M_{\odot}$ as black curves for comparison. We also indicate the locations of Betelgeuse (purple diamond; \citealp{dolan16}), the RSG progenitor of NGC6949-BH1 (thin black diamond; model P5 from \citealp{adams17a}), and the collection of directly detected SN progenitors of Types II (light gray squares) and IIb (dark gray circles) from \citet{smartt15}.
}
\end{figure}

Even for an early M-type star ($T_{\mathrm{eff}} \approx 3700$), where there is overlap between the two methods, the inferred luminosity from SED fitting is about a factor of 1.3 higher than from applying $BC_K$ for the same extinction of $E(B-V) = 0.8$\,mag. This modest discrepancy could be resolved if the effective extinction law were grayer---either from a shallower ISM law (i.e., $R_V > 3.1$) for the central regions of M51 or perhaps related to properties of the stellar wind or circumstellar dust not captured by the models (e.g., nonspherical geometry or nonisotropic scattering by dust grains)---resulting in a smaller correction at bluer wavelengths and a lower bolometric luminosity for a given model. Alternatively, using a smaller $BC_K$ at earlier spectral types (as found for cool supergiants in the Magellanic Clouds; 
\citealp[e.g.,][]{elias85,levesque06,davies18b}) would also bring the luminosity estimates closer to agreement. 

Altogether, we infer that M51-DS1 is likely a luminous RSG at $\log(L/L_{\odot}) \approx 5.2$--$5.3$ (for $T_{\mathrm{eff}} \lesssim 3700$~K and reasonable values of $E(B-V)=0.4$--1.3) corresponding to a highly evolved star of initial mass $M \approx 19$--22\,$M_{\odot}$, or 
somewhat warmer ($T_{\mathrm{eff}} = 3700$--4700~K)---possibly a YSG on a post-RSG track back to higher temperatures
---at higher luminosities up to $\log(L/L_{\odot}) \approx 5.4$--$5.7$. This more extreme case would suggest an initial mass as high as $M \approx 26$--40\,$M_{\odot}$ for a star near the empirical Humphreys-Davidson limit in an exceptionally short-lived and rare state of stellar evolution.

\section{Summary and Conclusions} \label{sec:summary}
We have presented a detailed analysis and characterization of a remarkable dimming event of a luminous, cool supergiant star in the nearby, star-forming galaxy M51, found as part of a new search for failed SNe as ``disappearing'' massive stars with \textit{HST}. The object, which we call M51-DS1, was detected dozens of times in the wealth of archival \textit{HST} imaging for more than two decades since 1995 before undergoing a dramatic, near-total disappearance in the 2019 ACS/WFC $F814W$ images obtained as part of our search. Thus meeting our criteria for a strong failed SN candidate (outlined in Section~\ref{sec:candidates}), we conducted follow-up \textit{HST} imaging observations with a Cycle 28 mid-cycle program using ACS/WFC and WFC3/IR, supplemented with ground-based near-IR imaging from Gemini-N/NIRI and MMT/MMIRS, to attempt to confirm the disappearance of the star. The star was found to have partially rebrightened to near its pre-disappearance levels across optical and near-IR, ruling out the terminal collapse of the star in a failed SN. As summarized below, our analysis indicates instead that an isolated episode of enhance mass loss on a very massive, cool Y/RSG star---potentially a rare, more extreme cousin of the recent ``Great Dimming'' of Betelgeuse---can explain its observed properties. 

M51-DS1 is located in the inner regions of M51 along the edge of a spiral arm. Previous studies of the host galaxy suggest an environment of approximately solar or somewhat supersolar metallicity ($\log[Z/Z_{\odot}] \approx 0.0$--0.3) and significant foreground extinction ($E[B-V] \approx 0.4$--1.0\,mag). Photometry of stars in the immediate vicinity indicate the presence of a young population ($\lesssim 6$\,Myr old) within a 20~pc radius. M51-DS1 itself is the brightest red source in this region, implying an initial mass of at least $\approx$30\,$M_{\odot}$ if it is part of the same population. It is difficult, however, to disentangle possible confounding effects, including a mix of stellar populations of different ages and high spatial variation in the foreground extinction toward individual sources up to $E(B-V) \approx 2.0$\,mag. At $M_K = -11.1$\,mag in 2005 and with $J-K\approx1.0$--1.2\,mag (Galactic extinction correction only), the star is consistent with a luminous RSG, and applying the $K$-band bolometric correction of \citet{davies18a} yields $\log(L/L_{\odot}) = 5.2$--5.3, implying an initial mass $M = 19$--22\,$M_\odot$ (assuming an M-type supergiant with $T_{\mathrm{eff}}\lesssim3700$~K and foreground extinction $E[B-V]=0.4$--1.3\,mag). Modelling of the SED suggests higher temperatures ($T_{\mathrm{eff}}\approx3700$--$4700$) and luminosities ($\log[L/L_{\odot}] = 5.4$--5.7), pointing to a more massive (initial mass $M = 26$--40\,$M_\odot$) YSG or post-RSG star that would be consistent with the inferred age of the very young stellar population in the immediate vicinity and possibly requiring significant circumstellar extinction. 

The extended, 26-yr light curve of M51-DS1 indicates a history of variability with $\Delta F814W \approx 0.7$--0.9\,mag, consistent with semiregular variations typical of Y/RSGs and expected theoretically for massive cool stars above $\gtrsim$15$\,M_{\odot}$. At $\Delta F814W \gtrsim 2.1$\,mag, we find that the 2019 minimum is unlikely to be part of the typical variation cycle for this star, but rather an uncommon and exceptionally deep dimming event---notably, more than $\gtrsim$1.5\,mag deeper than the 2019-2020 historic minimum of Betelgeuse in the comparable $I$-band. In the context of other large-amplitude modes of variation seen in luminous red stars in M51, we suggest that the 2019 event of M51-DS1 could be associated with an enhanced episode of mass loss from the star, in which increased opacity from the dense, molecular wind or the formation of new dust grains temporarily obscured the star.

This discovery highlights a central challenge for the definitive identification of failed SNe, namely, that the massive stars in question exhibit optical variability that can mimic a disappearing star. In particular, episodic mass loss in massive, cool supergiants remains poorly understood \citep[e.g.,][]{smith14}, both in terms of the physical mechanisms that drive it and critical properties such as the frequency, possible duration, and mass-loss rates of the most extreme events---an extended mass loss event could plausibly obscure a massive star in the optical for years. While Galactic examples of episodic mass loss from cool supergiants, like Betelgeuse's Great Dimming, can be studied in real time and in exquisite detail, such events will be few and far between. Still, significant progress can be made in mining the ever-growing trove of archival data from both space-based (e.g., \textit{HST}, \textit{Spitzer}) and numerous ground-based time-domain surveys to uncover and characterize the population of such sources in nearby galaxies, both as possible contaminants in ongoing failed SN searches and as direct probes of episodic mass-loss in massive, evolved stars. At the same time, a continued, concerted effort to monitor existing and newly uncovered failed SN candidates is vital, especially observations in the near-IR and mid-IR, e.g., soon with the \textit{James Webb Space Telescope} and later with the \textit{Nancy Grace Roman Space Telescope}, to test for surviving obscured stars and direct signatures of a potential obscuring wind or dust. 


\begin{acknowledgements}
Time domain research by D.J.S.\ is also supported by NSF grants AST-1821987, 1813466, 1908972, \& 2108032, and by the Heising-Simons Foundation under grant \#2020-1864. 
J.S.\ acknowledges support from NASA grant HST-GO-15645.003-A and the Packard Foundation.
Research by S.V.\ is supported by NSF grants AST-1813176 and AST-2008108. E.R.B.\ is supported by NASA through Hubble Fellowship grant HST-HF2-51428 awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA, under contract NAS5-26555. 

Based on observations made with the NASA/ESA Hubble Space Telescope, obtained at the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS5-26555. These observations are associated with programs \# HST-GO-15645, 16508, 10452, 14704, 12490, 12762.Support for program \# HST-GO-15645 and \# HST-GO-16508 was provided by NASA through a grant from the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS5-26555.

Observations reported here were obtained at the MMT Observatory, a joint facility of the University of Arizona and the Smithsonian Institution.

Based on observations obtained at the Gemini Observatory (Programs GN-2005A-Q-49,GN-2021A-DD-101), which is operated by the Association of Universities for Research in Astronomy, Inc., under a cooperative agreement with the NSF on behalf of the Gemini partnership: the National Science Foundation (United States), National Research Council (Canada), CONICYT (Chile), Ministerio de Ciencia, Tecnolog\'{i}a e Innovaci\'{o}n Productiva (Argentina), Minist\'{e}rio da Ci\^{e}ncia, Tecnologia e Inova\c{c}\~{a}o (Brazil), and Korea Astronomy and Space Science Institute (Republic of Korea). 

We acknowledge with thanks the variable star observations from the AAVSO International Database contributed by observers worldwide and used in this research.

\end{acknowledgements}

\facilities{HST (ACS, WFC3, WFPC2), Gemini:Gillett (NIRI), MMT (MMIRS), AAVSO}


\software{\texttt{AstroDrizzle}, \texttt{TweakReg} (\url{http://drizzlepac.stsci.edu/}; \citealp{hack12}), \texttt{SExtractor}, \texttt{PSFEx} (\url{https://www.astromatic.net/software/};
\citealp{bertin96,bertin11}), \texttt{DOLPHOT} (\url{http://americano.dolphinsim.com/dolphot/}; \citealp{dolphin00,dolphin16}, \texttt{DRAGONS} (\url{https://dragons.readthedocs.io/en/v2.1.1/index.html}), \texttt{Astropy} (\url{https://www.astropy.org/}; \citealp{2013A&A...558A..33A,astropy}), \texttt{photutils}, \texttt{EPSFBuilder} (\url{https://photutils.readthedocs.io/en/stable/}; \citealp{bradley20_photutils})}


