text:
\section{Introduction}\label{sec:intro}

Since 2015, the Advanced Laser Interferometer Gravitational-Wave Observatory detectors (LIGO, \citealt{ligo09}), and later the Advanced Virgo Observatory \citep{,virgo15} detectors, have been detecting gravitational waves (GWs) produced by the inspirals of merging compact objects (neutron stars: NS and/or black holes: BH). During the first two observing runs (O1 and O2), the LIGO-Virgo Collaboration (LVC) announced 10 binary black hole (BBH) mergers and one binary neutron star (BNS) merger \citep{2019O1O2Catalog}. The single BNS merger, termed GW170817 \citep{lvc_gw170817}, also produced electromagnetic (EM) emission (AT2017gfo) from the $\gamma$-rays to the radio band, signaling the first GW-EM multi-messenger discovery \citep{LIGO_MMA}. AT2017gfo was discovered in the galaxy NGC~4993 \citep{Coulter17} at a distance of $40.7\pm1.4\pm1.9$ Mpc (random and systematic errors; \citealt{Cantiello2018}).

With the discovery of this EM counterpart, GW170817 was detected across the entire EM spectrum with a multitude of telescopes (c.f., \citealt{LIGO_MMA}). In-depth studies of the EM emission and its spectral and temporal evolution have provided a wealth of information about BNS mergers, such as direct evidence of the radioactive decay of r-process material demonstrating kilonovae as important sites for the production of heavy elements \citep{Arcavi+17,Chornock+17,Coulter+17,Cowperthwaite17,Drout+17,Diaz+17,Fong+17,Gall+2017,Hu+17,Kasliwal+17,Lipunov+17,McCully+17,Nicholl17,Pian+17,Pozanenko+17,Shappee+17,Smartt+17,Soares-Santos+17,Tanvir+17,Utsumi17,Valenti2017,villar17}. Through late-time follow-up and detailed modelling, the identification of a structured relativistic jet was also confirmed for GW170817 \citep{Haggard17,hallinan17,margutti17,margutti2018,Alexander2018,mooley18a,Fong19,ghirlanda19,hajela19,troja19}.

The third observing run (O3) began on April 1, 2019 with a 1-month commissioning break in October 2019 and had an effective operational time of 10.9~months. Relative to the capabilities of O2, the increased sensitivities during the first 6 months of O3 (O3a) led to a median orientation-averaged range for BNS mergers of 108 Mpc, 135 Mpc and 45 Mpc for LIGO Hanford, LIGO Livingston, and Virgo, respectively; while the median 90\% localization during O3a was 860 deg$^2$ \citep{gwtc2}. During O3, a total of 56 non-retracted public alerts\footnote{https://gracedb.ligo.org} were released. Assuming an astrophysical origin for these 56 alerts, each GW candidate is assigned a probability of being the merger of two black holes (BBH), two neutron stars (BNS), a neutron star and a black hole (NSBH), or a more ambiguous case in which one of the components has a mass in the gap between NS and BH classifications (termed ``MassGap'', $3 < M/M_{\odot} < 5$; see \citealt{lvc_loc}). A fifth classification encompasses unmodeled short duration (t$<$1s) ``Burst'' signals, which do not fit the normal inspiral of compact objects and are of uncertain astrophysical origin.

This classification scheme resulted in 37 BBH, 7 BNS, 7 NSBH, 4 MassGap, and 1 burst event candidates that were detected during O3\footnote{Categorized by the highest non-terrestrial classification reported by the LVC for each event.}. For each event, we adopt the naming convention used by the LVC. Events have superevent names given by the letter ``S'' followed by the six digit UTC date and a lowercase incrementing letter. Events that are confirmed as real GW events by the LVC are renamed to ``GW'' followed by the six digit UTC date (e.g. GW190425). For O3a sources presented in the Gravitational-Wave Transient Catalog 2 (GWTC-2, \citealt{gwtc2}), we change the names to follow the naming convention adopted there. For events above the adopted False Alarm Rate (FAR) threshold, names include ``GW'' followed by the six digit UTC date and the six digit UTC time (e.g., GW190408\_181802); while sub-threshold events retain their superevent names. Previously published events (e.g. GW190425 and GW190521) keep their published names without the UTC time. These 56 non-retracted  events had 90\% localization regions that ranged in size from 18--24264 deg$^2$, depending largely on the number of detectors operational at the time. Other than the event detected by the burst pipeline (S200114f), which does not have an associated distance, the remaining 55 event candidates have inferred median distances of 108--5263~Mpc from the GW signals.

O3 comprised a number of firsts for the LVC, including the first detection of a candidate NSBH merger, GW190426\_152155 \citep{gcn24411}; the first detection of a candidate MassGap merger, S190814bv \citep{gcn25324}, which was later classified as a NSBH with a high mass-ratio, challenging current formation models and the predicted mass distribution of compact-object binaries (GW190814, \citealt{190814_paper}); and the first detection of a burst candidate, S200114f \citep{gcn26734}. In addition, the second BNS merger, GW190425 \citep{GW190425}, was significantly more massive than GW170817, as well as more massive than the Galactic NS binary population \citep{Safarzadeh2020}. Finally, O3 included the most massive BBH merger, GW190521, with BH masses of 91.4$^{+29.3}_{-17.5}$ $M_{\odot}$ and 66.8$^{+20.7}_{-20.7}$ $M_{\odot}$ (\citealt{lvc_gw190521,gwtc2}), having implications on various formation channels for these systems \citep{190521_imp_2020}.

Despite the wealth of new discoveries and GW event alerts, the large sky localizations coupled with the large median distances made it challenging for effective searches of associated UV, optical, or near-infrared (NIR) emission for a large number of events. Indeed, optical searches were dominated by instruments with wide fields-of-view (FOV), which largely adopted the approach of tiling the GW localization regions. Published examples include DECam/DES \citep{garcia20}, DECam/GROWTH \citep{Andreoni2019_190510,Andreoni20,Goldstein19}, ZTF/GROWTH \citep{coughlin19,Anand20,coughlin20,kasliwal20}, GRANDMA \citep{antier20}, GOTO \citep{gompertz2020}, ENGRAVE \citep{ackley2020} and CFHT \citep{vieira_cfht}; while other surveys like the MASTER-Network \citep{lipunov17}, Pan-STARRS \citep{chambers16}, ATLAS \citep{tonry18} and {\it Gaia} \citep{GaiaCo} have also contributed towards candidate searches. Meanwhile, telescopes with smaller FOV have opted for galaxy-targeted searches (Magellan; \citealt{gomez19}, MMT/SOAR; \citealt{Hosseinzadeh19}, J-GEM; \citealt{jgem}). The large majority of these surveys encompass triggered searches on GW events with a majority of the information distribution in real-time via The Gamma-ray Coordinates Network (GCN). Published papers on single events generally focus on the search for transient EM emission within a few days of the GW event and the in-depth analysis of the emission from candidates. Summary papers from groups such as ZTF/GROWTH \citep{Anand20,kasliwal20} have also explored the long term behaviour of their initially reported candidates, with efforts to follow up outstanding candidates through spectroscopic or photometric means.

Searches After Gravitational waves Using ARizona Observatories (SAGUARO) is a telescope network dedicated to the discovery and follow up of GW events\footnote{Although SAGUARO currently focuses on the follow-up of GW events, the methods presented in this paper could be applied to any transient of interest.}. In this paper, we present a summary of the SAGUARO campaign and analysis of our observations to find and study the optical counterparts to GW events that ran concurrent with LIGO/Virgo O3 \citep{Lundquist19}. In Section \ref{sec:overview} we provide an overview of SAGUARO, including our image subtraction and candidate vetting, summarize our coverage in O3, and discuss significant updates implemented during O3. We present a detailed description of our observations and each individual event in Section \ref{sec:ind}. In Section \ref{sec:discussion}, we discuss our findings. Lastly, we discuss the conclusions, lessons learned from O3, and future prospects in Section \ref{sec:conclusion}.

Unless otherwise stated, all magnitudes reported here are in \emph{Gaia} DR2 $G$-band and are converted to the AB system via $m_{AB}$ = $m_{Gaia}$ + 0.125 \citep{Maiz18}. We employ standard $\Lambda$CDM cosmology with H$_{0}$ = 69.6 km s$^{-1}$ Mpc$^{-1}$ and $\Omega_{M}$ = 0.286 \citep{Bennett+14}.

\begin{figure*}
\includegraphics[width=0.8\textwidth]{scatter_limits.png}
\centering
\vspace{-0.2in}
\caption{The 5$\sigma$ transient depth of SAGUARO fields measured from the image subtraction process versus $\delta t$ (defined as the time after the GW event) for 17 events with triggered or serendipitous fields; all symbols represent upper limits. The reported depths (in various filters and sigma levels) of other GW follow-up surveys: GOTO \citep{gompertz2020}, ATLAS \citep{tonry18}, the MASTER-Network \citep{lipunov19}, ZTF/GROWTH \citep{kasliwal20}, and  Pan-STARRS \citep{chambers16}, are shown as black horizontal lines. We also note a 23.1~mag $i$-band median depth of the CFHT survey \citep{vieira_cfht}, a 23.0~mag $g$- or $r$-band median depth of the DECAM/GROWTH searches \citep{Anand20}, a 22.4~mag $i$-band median depth of DESGW (based on 2 events; \citealt{GCN27366, Morgan20}), and a 17.4~mag median limiting magnitude of BNS and NSBH events by the GRANDMA collaboration \citep{antier20}. The median 5$\sigma$ transient depth of SAGUARO fields, 21.1 mag, is deeper than many of the other GW follow-up surveys. \label{fig:Appmag}}
\end{figure*}

\begin{figure*}
\centering
\hspace{-2cm}
\includegraphics[width=0.8\textwidth]{candidate_vetting.png}
\caption{Flow chart describing the automatic and human vetting of SAGUARO candidates during the in-depth candidate search presented here and described in Section \ref{sec:vetting}. All candidates from the image subtraction process using ZOGY \citep{ZOGY} with S/N $>$ 5 are cross-matched to various databases (such as the Minor Planet Center Orbit Database and {\it Gaia}) and assigned a machine learning (ML) score automatically after ingestion. Through a web interface, human vetting of candidates proceeds for all candidates with ML $>$ 0.7. This process includes removing any remaining artifacts or slow-moving objects, cross-checking with PanSTARRs images, checking the redshift of the candidate if in a galaxy, checking the CSS/ZTF light curve including the time of first detection, and community follow-up. A more in-depth study of the candidates is then performed. The details on the final candidates is summarized in Table \ref{tab:obs}. \label{fig:vetting}}
\end{figure*}

\section{Overview of SAGUARO} \label{sec:overview}
SAGUARO started operations in April 2019 and uses a network of telescopes\footnote{The facilities range from the 1.54 m Kuiper, 1.8 m VATT and 2.3 m Bok telescopes in Southern Arizona, to the 6.5 m MMT, the twin 6.5 m Magellan telescopes, the 2$\times$8.4 m Large Binocular Telescope, and the two Keck 10 m telescopes.}, primarily in Arizona. To facilitate the immediate follow-up of GW events, SAGUARO makes use of the Steward Observatory 1.5 m Mt. Lemmon telescope operated by the Catalina Sky Survey \citep[CSS;][]{CSS}. The telescope is equipped with a prime focus imager and a 10.5K$\times$10.5K CCD (0.77\arcsec\ per pixel), resulting in a 5 deg$^2$ FOV. It is operated with 2$\times$2 binning for an effective plate scale of 1.54\arcsec\ per pixel. CSS observes fixed fields on the sky, between declinations of $-$25~deg and +60~deg, while avoiding crowded regions in the Galactic plane (see Figure 1 of \citealt{Lundquist19}) and observes $\sim 24$ nights per month, avoiding the period around full moon. During normal operations, CSS searches for near-Earth objects (NEOs) and potentially hazardous asteroids (PHAs). To accomplish this, CSS obtains 4 $\times$ 30 s exposures per field, observing 12 fields over the span of $\sim$ 30 minutes (thus the typical survey speed is 120 deg$^2$ per hour; $\gtrsim$ 1000 deg$^2$ per night). This results in four images per field separated by $\sim$ 8 minutes, allowing for the identification of moving objects. SAGUARO observations that are triggered for GW counterpart searches (``triggered fields'') follow the same observing strategy. These fields are typically triggered in sets of 60, 120, or 180 deg$^2$ on the highest probability region of the GW localization that is observable. All images are taken without a filter and are calibrated to {\it Gaia} $G$-band using {\it Gaia} DR2 \citep{Gaiadr2}. For each 4$\times$30~s set of median-combined images, we show the 5$\sigma$ transient detection depth, defined as the magnitude of a source that can be detected through image subtraction with a significance (S/N) = 5, in Figure \ref{fig:Appmag}. The median 5$\sigma$ transient detection depth is $G\approx$ 21.1 mag. For astrometry, images are tied to {\it Gaia} DR2. The details of the follow-up capabilities, as well as an example of such follow-up, can be found in \citealt{Lundquist19}.

\subsection{Image subtraction and candidate vetting} \label{sec:vetting}
Deep reference images were constructed based on nearly three years of CSS data, with an average limiting magnitude of $G\approx 23.0$~mag. Data are run through a custom pipeline in real-time that performs image subtraction on each 4$\times$30 s median-combined image using ZOGY \citep{ZOGY}. Candidates with a S/N $>$ 5 are ingested into our database. PSF photometry is performed directly on the subtracted image, following the equations in \cite{ZOGY}, and calibrated to the CSS zero point calculated from {\it Gaia} DR2. A graphical summary of the entire candidate vetting process is shown in Figure \ref{fig:vetting}. ``Automatic vetting'' and simple ``human vetting'' were performed for triggered SAGUARO fields in real-time, with the results reported via GCN to the community. In this paper, we present a more detailed vetting of candidates in the context of the updates discussed in Section \ref{sec:updates}.

We make use of various catalogs and databases, including the Minor Planet Center Orbit Database\footnote{\url{http://www.minorplanetcenter.net/iau/MPCORB.html}} (MPCORB) to identify known moving objects, and the Transient Name Server\footnote{https://wis-tns.weizmann.ac.il/} (TNS) and the Zwicky Transient Facility (ZTF, \citealt{ZTF}) to identify known transients. We also cross-match candidates with the Galaxy List for the Advanced Detector Era (GLADE; \citealt{glade}) for possible associations to nearby cataloged galaxies, and the {\it Gaia} DR2 catalog \citep{Gaiadr2} to reduce contamination by stellar sources (see Section \ref{sec:updates}). Candidates are then assigned a machine learning (ML) score, based on the likelihood of being a real astrophysical source. A cut of ML score $>$ 0.7 is made before continuing (see Section \ref{sec:updates}). This process falls under our ``automatic vetting'' and is done in real-time without human intervention to prepare candidates for ``human vetting''.

We then move onto our human vetting through our web interface. We manually remove any remaining artifacts caused by satellite trails, bright stars or optical reflections. Slow moving objects which are not removed during the median stacking of the four images or found in the MPCORB were then rejected using our moving object calculator (see Section \ref{sec:mo}). Next, candidates coincident with point-like sources that we infer to be stellar with PanSTARRS1 multi-band images \citep{chambers16} that were not removed by the \emph{Gaia} cross-matching, were rejected as potential kilonova candidates. If candidates are clearly spatially coincident with a galaxy, the redshift of the host galaxy (if known) is compared to the distance inferred for the GW event reported by LIGO/Virgo. Candidates with detections prior to the GW event in either CSS or ZTF were ruled out as being related to the GW event. Candidates with light curves inconsistent with the rise or fall times of kilonovae (see Section \ref{sec:kilonova}) built from a combination of CSS and ZTF data were rejected. Additionally, candidates already ruled out by the community through spectroscopic or photometric means (reported either through GCNs or papers) were discounted. The remaining candidates are then analyzed in more detail to assess their association with the GW event and provide value-added information. First, individual images are checked to assess shifts in the source's shape and position, as well as the presence of the candidate in each image, to confirm the candidate as a real source. Next, a literature search for spatially coincident galaxies using a radius, $r < 4''$ (equating to spatial scales of $\sim$ 2--13 kpc for distances of 100--500 Mpc) is performed. If no galaxy is found within the literature, we search the Legacy Survey DR8 \citep{Dey2019} data for a detection or limits (based on nearby detections) on a faint coincident host. For candidates embedded within a galaxy, the position of the candidate and offset, $\delta r$, is reported relative to the position of the galaxy. For galaxies without measured redshifts, we check PanSTARRS DR2 \citep{PSDR2} for photometric redshifts ($z_{photo}$). The final candidates for all events presented here are summarized in Table \ref{tab:obs}, with the details discussed in-depth in the relevant individual event sections (Sections \ref{sec:S190901ap}, \ref{sec:S190930s}, \ref{sec:S190930t}, \ref{sec:S200105ae}, \ref{sec:S200224ca}).

In addition to vetting our own candidates, we make use of normal CSS operations to rule out candidates reported by the community that were not followed-up and are covered by CSS at a later date. For candidates within the CSS footprint, the positions of these candidates are put through the same automatic and human vetting procedures as described above to assess their viability as genuine kilonova.

\subsection{Coverage in O3} \label{sec:obs}
Starting alongside LIGO/Virgo's O3, SAGUARO had overlapping operations for 9.5 months with a duty cycle of 53.6\% for the 24~hr window following each event and 78.6\% for the 120~hr window after each event\footnote{The calculation of the duty cycle takes all observing constraints into account including the annual Arizona `monsoon' shutdown running mid-July to the end of August, the 1 month LIGO/Virgo shutdown between O3a and O3b, CSS moon constraints, and weather.}. Including all observing constraints, CSS was available to trigger on 27 non-retracted events with $\delta t < 24$~hrs (where $\delta t$ is defined as the time after the event) and on 38 events within $\delta t < 120$~hrs. Our trigger criteria focused on events which contained potential NS components (where the sum of the classifications containing a possible NS, $\sum_{\mathrm{(BNS,NSBH,MassGap})} >20\%$ at the time of trigger) for which CSS could cover part of the 90\% probability region when considering airmass and moon constraints described in \cite{Lundquist19}; but also included unusual events such as S200114f, the ``burst'' event. Two additional components of importance ``HasNS'' and ``HasRemnant'' accompany each event. These describe the probability that the event has a NS and the probability that there is matter outside the final compact object (calculated based on the masses and spins inferred), respectively. With these classifications, we would expect events with large HasRemnant values to produce EM emission. Taking our trigger criteria, the GW sky localization, and data quality under consideration, we present 7 triggered events during O3, including an initial trigger to test our system.

\section{SAGUARO System Updates} \label{sec:updates}
In the following section, we describe significant updates to the SAGUARO system for the latter part of O3 which build upon the infrastructure outlined in \citealt{Lundquist19}. First, we discuss the improvement of the machine learning algorithm, allowing for more accurate classifications between real and bogus candidates. Second, we include cross-matching with \emph{Gaia} proper motions measurements, to remove candidates clearly associated with stellar sources. Third, we introduce a more stream-lined approach to identify new, uncatalogued moving objects, with particular attention to slower moving objects that persist in the median images as candidates. Fourth, we implement infrastructure to identify and search through CSS fields which serendipitously overlap with observable GW events within $\delta t<5$~days, designated as ``serendipitous fields''. And lastly, we discuss the use of kilonova models to place constraints on the photometric evolution of candidates.

\subsection{Machine Learning}
Candidates are given a ML score between 0 (bogus) and 1 (real) based on a machine learning model to reduce the number of bogus candidates. Following the example of \citet{wright15}, this model was created using a random forest classifier that was trained on 10$\times$10 pixel substamps of a sample that included 3678 detections of known transients and 56062 bogus detections, including artifacts from bright stars and over-subtractions. To reduce overfitting in our machine learning model and improve generalization, we add Gaussian noise into our training sample data at the level of 1\% of the peak brightness. To improve rotation invariance, we add versions of each image in our training sample that have been rotated by 90, 180, and 270 degrees, effectively increasing our training set by four-fold.

We evaluated the ML scores on a test sample that included 1285 detections of known transients. By selecting a lower limit on the ML score of 0.70, such that we only assess candidates with ML $>$ 0.70, we recovered 94.3\% of these known transients with an average score of 0.87. Given this recovery rate of known transients, we use this score to limit the number of candidates that need human vetting.

\subsection{Gaia Matching}
To reduce the number of known variable stellar sources that produce residuals in subtracted images, all candidates are cross-matched against high proper motion objects in the {\it Gaia} DR2 catalog \citep{Gaiadr2} upon ingestion. If a candidate matches as a high proper motion source ($\mu > 3\sigma$) in the \emph{Gaia} catalog, it is classified in the candidate database as a stellar object and thus removed from our candidate vetting process. This reduces the number of candidates that need to be human vetted by $\sim$7.2\%.

\subsection{Moving Objects} \label{sec:mo}
Candidates are cross-matched against known moving objects from the MPCORB to remove these before human vetting. However, unknown moving objects, especially extremely slow-moving objects that move on the order of a pixel or less between each image and thus are not removed in the median image, can still persist as potential candidates. To help eliminate these slow-moving objects whose small shifts can be difficult to determine by eye, we implemented a moving object calculator. This calculator uses the position of the source in each of the four individual images, and determines the shift between the first and each of the subsequent images along with the standard deviation of these shifts. With this method, we are able to identify moving objects that move as little as $\sim$ 1 pixel between each of the individual images i.e. $\sim$1.6$''$/8 min.

\begin{figure*}
\centering
\hspace*{-1.8cm}
\includegraphics[width=1.18\textwidth]{event_comparison.png}
\caption{Temporal and areal coverage of all GW events with significant follow-up from SAGUARO during O3. The histogram shows the fraction of the total area ($A_{\mathrm{total}}$, solid pink) and probability ($P_{\mathrm{total}}$, gold line) covered as a function of $\delta t$. The circles on the right show the area/probability covered by triggered fields (filled circle) normalized to the total covered area by including the serendipitous fields (open circle). As such, events with only a filled circle represent those without serendipitous coverage, while those with only the open circles show purely serendipitous events. \label{fig:comparison}}
\end{figure*}

\subsection{Serendipitous Fields}
While GW170817 peaked in the $r$-band at $\delta t\approx 0.65$~days and faded at a rate of $\sim$~1 magnitude per day  \citep[e.g.][]{villar17,kasen17,Nicholl17,Cowperthwaite17,Valenti2017}, kilonova models describing various progenitor and remnant scenarios predict a diversity in peak times and fading behaviors (e.g., \citealt{metzger14,kasen15,barbieri20}; see Section \ref{sec:kilonova} and Figure~\ref{fig:Absmag}). As CSS maps $\sim$ 1000 deg$^2$ of the sky each night in search of NEOs and PHAs, this provides the opportunity to analyze additional data outside our triggered area and time frame in search of optical counterparts. Thus motivated, we ingest incoming CSS survey data to search for fields which serendipitously fall within the 90\% contour of the localization of GW events within $\delta t <$ 5 days (i.e. $<$ 120 hr). We perform image subtraction between these fields and our deep reference images using ZOGY, and perform the same cross-matching and ingestion as described above. When only considering the triggered fields, SAGUARO covers between 15--180 deg$^2$ for each triggered event. The inclusion of these serendipitous fields results in an increased coverage of up to 2115 deg$^2$ for purely serendipitous events. For events with triggered fields, we see up to $\sim$ 9 times in additional area coverage (see Figure \ref{fig:comparison}). In total, we have serendipitous coverage up to $\sim$43\% of the LVC total probability ($P_{\mathrm{total}}$) for individual events. Thus, we present an additional 10 events with serendipitous fields, defined as those with coverage in the range of $\delta t < 120$~hr.

\begin{figure*}
\centering
\includegraphics[width=0.8\textwidth]{scatter_absmag_models.png}
\vspace{-0.2in}
\caption{The 5$\sigma$ transient depth in absolute magnitude space of individual SAGUARO fields with five kilonova models in the $r$-band. The solid black line represents a model for a GW170817/AT2017gfo-like kilonova \citep{kasen17}. The dashed black line shows a modeled kilonova from a NSBH merger viewed at $41<\theta<46$ degrees off-axis \citep{kawaguchi20}. The dotted line is a model for the unusually massive event GW190425 in the case of a BNS merger \citep{foley20}. The grey band models a slow-rising kilonova in the case of a NS remnant \citep{kasen15} for a prompt collapse (lower boundary), lifetime of 100 ms and an infinite lifetime (upper boundary). Finally, the blue band is a more optimistic model for the unusually massive event GW190425 in the case of a NSBH merger (upper boundary) and a BNS merger (lower boundary) \citep{barbieri20}. SAGUARO observations for BNS, NSBH and MassGap events are scaled to the median inferred distance reported by LIGO/Virgo. \label{fig:Absmag}}
\end{figure*}

\subsection{Kilonova models} \label{sec:kilonova}
To place constraints on photometric evolution, for events classified as BNS, NSBH, or MassGap, we compare each candidate to a number of kilonova models which represent a large range of progenitor and remnant scenarios. To describe the data of GW170817/AT2017gfo, we choose the two-component model described in \cite{kasen17} which we smooth with a Savitsky-Golay filter. For NSBH events, we compare to a model by \cite{kawaguchi20} parameterized by an ejecta mass $M_{\rm ej} = 0.06 M_{\odot}$\footnote{This model divides the ejecta mass into dynamical ejecta ($M_{\rm dyn} = 0.02 M_{\odot}$) and post-merger ejecta mass ($M_{\rm pm} = 0.04 M_{\odot}$).}. We correct for the $\sim$0.2 mag artificial excess brightness (determined through private communication with the authors) due to the use of a restricted line list in calculating the ejecta opacity and smooth using a Savitsky-Golay filter. We also compare to 3 BNS models with resulting hyper-massive NS (HMNS) remnant lifetimes of 0~ms (prompt formation of a BH), 100~ms, and infinite (corresponding to a indefinitely stable remnant; \citealt{kasen15}). The overall effect of a long-lived NS remnant, as opposed to prompt collapse to a BH, is that the remnant radiates neutrinos into the post-merger disk, resulting in a bright and blue slow-rising kilonova \citep{Metzger17, perego14, metzger14}. Finally, to describe a massive merger event like GW190425 (with a chirp mass of $1.44^{+0.02}_{-0.02} M_{\odot}$; \citealt{GW190425}), we choose 2 models. The first model's ejecta masses are $\gtrsim 0.1 M_{\odot}$ with the DD2 equation of state, representing an optimistic transient brighter than AT2017gfo \citep{barbieri20}. The second assumes ejecta mass $M_{\rm ej} = 0.04 M_{\odot}$, velocity $v_{\rm ej} = 0.15c$, and lanthanide fraction $X_{\rm lan} = 10^{-2}$; and predicts a kilonova less luminous than AT2017gfo across all bands \citep{foley20}. We show these models in rest-frame time units ($\delta t_{rest}$) in Figure \ref{fig:Absmag}, along with the absolute magnitude of our 5$\sigma$ transient depth for individual fields. Even for the slowest rising model, given by an indefinitely stable remnant, the peak of emission occurs at $\delta t \sim$ 2.5 days. We can thus use the rise times to rule out candidates which continue to rise for $\delta t \gg 3$ days. Given the time evolution of all models, we can also use the fading time to rule out candidates which have detections on longer timescales, i.e. for candidates detected months after the GW event. In general, SAGUARO is able to search to depths that would allow the detection of kilonova emission out to $\sim$ 150--400 Mpc, assuming detection at peak brightness and depending on the model assumed (see Figure \ref{fig:Absmag}).

\begin{deluxetable*}{lccccccccccc}
\centering
\tabletypesize{\scriptsize}
\tablecolumns{15}
\tablewidth{0pc}
\tablecaption{Summary of SAGUARO Follow-up in O3
\label{tab:follow-up}}
\tablehead{
\colhead {Event}	 &
\colhead {Type$^\dagger$}	 &
\colhead {FAR}	 &
\colhead {Distance}	 &
\colhead {$\delta t_{\rm first}^\star$}	 &
\colhead {$\delta t$}	 &
\colhead {Fields}	 &
\colhead {5$\sigma$ Limit$^{\ddagger}$} &
\colhead {$A_{50\%}$} &
\colhead {$A_{90\%}$} &
\colhead {$A_{\mathrm{total}}$} &
\colhead {$P_{\mathrm{total}}$} \\
\colhead {}		 &
\colhead {}		 &
\colhead {(Hz)}			 &
\colhead {(Mpc)}			 &
\colhead {(hr)}			 &
\colhead {(hr)} &
\colhead {observed}		 &
\colhead {(AB Mag)} &
\colhead {(deg$^2$)}	&
\colhead {(deg$^2$)}	&
\colhead {(deg$^2$)}	&
\colhead {(\%)}	\\
}
\startdata
GW190408\_181802 & BBH & 2.8e-18 & 1472.90 $\pm$ 357.88 & 17.6 & $<24$ & 3 & 19.2 & 5.7 & 12.9 & 15.0 & 13.73 \\
 \hline
 \hline
GW190425 & BNS & 4.5e-13 & 156.14 $\pm$ 41.37 & 1.4 & $<24$ & 12 & 21.1 & 60.0 & 60.0 & 60.0 & 3.68 \\
 & & & & & $24-48$ & 12 & 20.5 & 0.2 & 48.5 & 60.0 & 0.23\\
 & & & & & $48-72$ & 31 & 21.2 & 72.0 & 149.6 & 155.0 & 3.57\\
 & & & & & $72-96$ & 6 & 20.8 & 12.2 & 26.9 & 30.0 & 0.46\\
 \hline
 & & & & & & & 21.1 & 124.4 & 265.0 & 285.0 & 6.55 \\
 \hline
 \hline
GW190426\_152155 & NSBH & 1.9e-8 & 376.72 $\pm$ 100.44 & 43.6 & $24-48$ & 12 & 20.6 & 17.3 & 54.4 & 60.0 & 4.48 \\
 \hline
 \hline
GW190521 & BBH & 3.8e-9 & 3931.42 $\pm$ 953.03 & 98.3 & $96-120$ & 7 & 21.3 & 0.3 & 30.5 & 35.0 & 2.27 \\
 \hline
 \hline
GW190630\_185205 & BBH & 1.44e-13 & 925.67 $\pm$ 258.51 & 62.7 & $48-72$ & 8 & 21.0 & 25.3 & 38.2 & 40.0 & 9.07 \\
  & & & & & $72-96$ & 8 & 21.1 & 21.0 & 36.3 & 40.0 & 7.14\\
  & & & & & $96-120$ & 10 & 20.9 & 26.5 & 46.1 & 50.0 & 9.09\\
 \hline
 & & & & & & & 21.0 & 72.7 & 120.7 & 130.0 & 25.29 \\
 \hline
 \hline
S190901ap & BNS & 7.0e-9 & 240.87 $\pm$ 78.65 & 11.9 & $<24$ & 12 & 21.1 & 60.0 & 60.0 & 60.0 & 1.66\\
 & & & & & $24-48$ & 38 & 21.0 & 101.4& 190.0 & 190.0 & 1.87\\
 & & & & & $48-72$ & 44 & 21.2 & 119.6& 219.6 & 220.0 & 2.06\\
 & & & & & $72-96$ & 15 & 21.3 & 0.0& 69.4 & 75.0 & 0.19\\
 \hline
 & & & & & & & 21.2 & 281.0 & 539.0 & 545.0 & 5.79 \\
 \hline
 \hline
S190923y & NSBH & 4.78e-8 & 438.09 $\pm$ 132.95& 118.6 & $96-120$ & 9 & 21.2 & 28.3 & 43.1 & 45.0 & 4.05  \\
 \hline
 \hline
GW190930\_133541 & BBH & 3.0e-9 & 708.90 $\pm$ 190.66 & 12.79 & $<24$ & 32 & 21.3 & 72.1 & 153.0 & 160.0 & 8.33\\
 & & & & & $24-48$ & 24 & 21.2 & 0.0& 106.8 & 120.0 & 3.79\\
 & & & & & $48-72$ & 11 & 21.1 & 0.0& 43.2 & 55.0 & 1.00 \\
 \hline
 & & & & & & & 21.2 & 72.1 & 303.0 & 335.0 & 13.13 \\
 \hline
 \hline
S190930t & NSBH & 1.5e-8 & 108.49 $\pm$ 37.67 & 11.80 & $<24$ & 91 & 21.2 & 74.3 & 451.5 & 455.0 & 1.51 \\
 & & & & & $24-48$ & 61 & 21.2 & 219.2& 305.0 & 305.0 & 1.58 \\
 & & & & & $48-72$ & 115 & 21.1 & 207.3& 567.9 & 575.0 & 2.04 \\
 & & & & & $72-96$ & 41 & 21.0 & 205.0& 205.0 & 205.0 & 1.55 \\
 & & & & & $96-120$ & 117 & 20.9 & 263.7& 580.4 & 585.0 & 2.47 \\
 \hline
 & & & & & & & 21.1 & 964.4 & 2104.8 & 2120.0 & 9.13 \\
 \hline
 \hline
S191205ah & NSBH & 1.2e-8 & 384.97 $\pm$ 163.72 & 3.60 & $<24$ & 27 & 21.0 & 68.8 & 130.1 & 135.0 & 10.00 \\
 & & & & & $24-48$ & 30 & 20.9 & 4.6& 147.3 & 150.0 & 0.91 \\
 \hline
 & & & & & & & 20.9 & 73.4 & 277.4 & 285.0 & 10.91 \\
 \hline
 \hline
S191213g & BNS & 3.55e-8 & 200.86 $\pm$ 80.96 & 47.4 & $24-48$ & 8 & 20.9 & 2.5 & 39.6 & 40.0 & 0.65 \\
  & & & & & $48-72$ & 18 & 21.0 & 13.8 & 85.0 & 90.0 & 1.40 \\
  & & & & & $72-96$ & 2 & 21.0 & 0.0 & 9.2 & 10.0 & 0.05 \\
 \hline
 & & & & & & & 21.0 & 16.2 & 133.8 & 140.0 & 2.10 \\
 \hline
 \hline
S200105ae & NSBH & 7.67e-7 & 282.85 $\pm$ 73.78 & 20.7 & $<24$ & 23 & 19.64 & 6.5 & 104.3 & 115.0 & 0.90\\
& &  &  &  & $24-48$ & 30 & 20.4 & 55.0 & 139.1 & 150.0 & 1.97\\
 \hline
 & & & & & & & 20.7 & 61.5 & 240.9 & 260.0 & 2.84 \\
 \hline
 \hline
S200114f & IMBH & 1.2e-9 & n/a & 2.36 & $<24$ & 36 & 21.0 & 37.7 & 176.1 & 180.0 & 86.30 \\
 \hline
 \hline
S200115j & MassGap & 2.09e-11 & 332.22 $\pm$ 78.46 & 72.5 & $72-96$ & 9 & 21.3 & 0.0 & 37.5 & 45.0 &  2.07\\
& &  &  &  & $96-120$ & 2 & 21.4 & 0.0 & 9.9 & 10.0 & 0.57 \\
 \hline
 & & & & & & & 21.3 & 0.0 & 47.4 & 55.0 & 2.64 \\
 \hline
 \hline
S200128d & BBH & 1.65e-8 & 3701.59 $\pm$ 1264.51 & 0.16 & $<24$ & 16 & 21.1 & 4.1 & 72.3 & 80.0 & 1.93 \\
& &  &  &  & $24-48$ & 15 & 21.0 & 13.2 & 71.5 & 75.0 & 2.48 \\
& &  &  &  & $48-72$ & 20 & 20.9 & 42.0 & 92.6 & 100.0 & 4.66 \\
& &  &  &  & $72-96$ & 13 & 21.3 & 25.9 & 62.8 & 65.0 & 2.73  \\
 \hline
 & & & & & & & 21.0 & 85.2 & 299.1 & 320.0 & 11.80 \\
 \hline
 \hline
S200224ca & BBH & 1.61e-11 & 1575.00 $\pm$ 322.36 & 32.3 & $24-48$ & 1 & 20.4 & 0.0 & 3.9 & 5.0 &  1.18  \\
& &  &  &  & $48-72$ & 4 & 21.1 & 0.0 & 14.9 & 20.0 & 10.24  \\
 \hline
 & & & & & & & 21.1 & 0.0 & 18.8 & 25.0 & 11.42 \\
 \hline
 \hline
S200316bj & MassGap & 7.10e-11 & 1177.98 $\pm$ 283.01 & 5.4 & $<24$ & 2 & 21.1 & 0.0 & 5.1 & 10.0 & 0.24\\
& &  &  &  & $24-48$ & 18 & 20.9 & 22.1 & 88.9 & 90.0 & 73.00\\
 \hline
 & & & & & & & 20.9 & 22.1 & 93.9 & 100.0 & 73.24 \\
 \hline
 \hline
\enddata
\tablecomments{Magnitudes reported here are uncorrected for Galactic extinction and are reported in {\it Gaia} $G$-band.\\
The probability covered refers to the percent of the probability of the 50\%, 90\%, and total localizations that were covered by these observations.\\
$^\ddagger$ Median $5\sigma$ transient detection depth. \\
$^\star$ $\delta t$ of first field after the GW event. \\
$^\dagger$ Most likely classification based on GW probabilities \citep{Kapadia19}. \\
$^a$ Triggered events.
}
\end{deluxetable*}

\section{Individual events and Results} \label{sec:ind}
Here we (i) present a more detailed analysis of the search for candidate counterparts in triggered fields, given the improvements to our automatic vetting portion of the pipeline (Section \ref{sec:updates}) using updated localizations and (ii) extend this analysis to the serendipitous fields for all events that have SAGUARO coverage with $P_{\mathrm{total}}>2\%$.

For triggered fields, we vetted candidates in real-time to search for optical counterparts. Information on these fields were provided to the community through GCNs \citep{gcn25915,gcn26360,gcn24079,gcn24172,gcn26473,gcn26750,gcn26753} and uploaded to the Gravitational Wave Treasure Map\footnote{\url{http://treasuremap.space}} \citep{Wyatt20}. During O3, we reported observations of 6 new candidates found by SAGUARO in the triggered fields (SN2019wdb, \citealt{gcn26360}; SAGUARO20a, SAGUARO20b, SAGUARO20c, SAGUARO20d, and SAGUARO20e, \citealt{gcn26753}). We were also able to rule out 2 community-identified candidates: the first, AT2020vx, due to a previous CSS detection \citep{gcn26750}, and the second, SN2019ebq, due to spectroscopic classification \citep{Lundquist19}.

We exclude 4 GW events from our analysis that had serendipitous coverage $P_{\mathrm{total}}<2\%$. Three of these are candidate BBH mergers (for which bright EM counterparts are less likely), while the fourth is a candidate NSBH merger, S190910d \citep{gcn25723}, which has a $d$ = 632 $\pm$ 186 Mpc (such that our imaging does not reach sufficient depths to make any meaningful statements on NSBH kilonova emission). We also exclude an additional 2 events, 1 triggered (S191216ap; \citealt{gcn26454}) and 1 serendipitous (S191215w; \citealt{gcn26441}) from our analysis due to poor data quality from bad weather.

Thus, in total, we present observations for 17~GW events, 7 triggered and 10 with only serendipitous coverage. In this section, we detail the GW discovery and relevant follow-up observations of these 17 events. For each, we summarize all relevant candidates reported by the community (focusing on actual candidates detected by the community versus survey limits), investigate whether they are still viable kilonova candidates with available information (checking GCNs, publicly available data uploaded to TNS, public ZTF data, peer-reviewed papers, and our own data), and present our own candidate search. We note SAGUARO follow-up of three of these events were also detailed in \citet{Lundquist19}; however, we include their details again here based on our updated system, results from additional serendipitous fields, as well as the follow-up of community candidates.

An overview of our observations, including the area covered of the 50\%, 90\%, and total localization ($A_{50\%}, A_{90\%}, A_{\mathrm{total}}$), as well as the percentage of the total probability ($P_{\mathrm{total}}$) covered for each event, is listed in Table~\ref{tab:follow-up}. The GW localization maps, along with our coverage, are displayed in Appendix \ref{appendix:A}. A summary of our final candidates is presented in Table \ref{tab:obs}. Hereafter, we use the common designations L1 (LIGO Livingston), H1 (LIGO Hanford), and V1 (Virgo).

\subsection{GW190408\_181802} \label{sec:S190408an}
A candidate GW signal was identified on 2019-04-08 at 18:18:02.288 UTC using data from L1, H1, and V1 \citep{gcn24069}. The final classification for this event is $>$99\% BBH at a distance, $d$ = 1473 $\pm$ 358 Mpc. This event, as the first of the O3 alerts, served as a test for our system, and we triggered 3 fields. The initial follow-up for these triggered fields is discussed in \cite{Lundquist19}. No serendipitous fields were covered by CSS within our search interval of $\delta t < 120$~hr. The localization of the event, along with the CSS fields that fall within the 90\% probability region, is shown in Figure~\ref{fig:190408an_loc}.

A total of 11 candidates were initially reported by the community via GCN (MASTER-Network, PanSTARRS, \emph{Gaia}; \citealt{gcn24076,gcn24084,gcn24096,gcn24124,gcn24134}). Initial follow-up, either through photometry or spectroscopy, ruled out 4 of these as viable candidates \citep{gcn24090,gcn24092,gcn24097,gcn24106,gcn24154}. SN2019daj was found to have detections prior to the GW event \citep{gcn24078} and later classified as a Type Ia SN \citep{2019daj}. Later follow-up classified another candidate, SN2019dma, as a Type II SN \citep{AT2019dma}. Another candidate, AT2019dpg, was detected by routine ZTF survey operations out to $\sim$ 35 days after the GW event \citep{Nordin2019}. Having faded by $<1.5$ mag over a month, we consider it to be unrelated to the GW event based on photometric evolution. Searching public ZTF data, we find another candidate, AT2019deu, shows variability for over a year after the GW event \citep{ztfarchive}. We thus consider it to be unrelated to the GW event. Of the remaining four candidates, CSS observed the positions of AT2019doo, AT2019dmc, and AT2019dev. We found no detections down to a limit of 20.8, 21.0, and 21.2 mag at $\delta t$ = 74.7, 65.7, and 65.8 days, respectively. With these limits we are not able to rule out these as potential candidates. No other viable candidates were found within our data.

\subsection{GW190425} \label{sec:GW190425}
A candidate GW signal with classification $>$99\% BNS was identified on 2019-04-25 at 08:18:05.017 UTC using data from L1 and V1 \citep{gcn24168}. This event was later confirmed to be a real BNS merger with a total mass of $3.4_{-0.1}^{+0.3}~\mathrm{M_{\odot}}$ that occurred at $d = 159_{-71}^{+69}\ \mathrm{Mpc}$ \citep{GW190425}. We triggered 12 fields in the northern part of the localization with $\delta t < 24$~hr. The initial follow-up of these triggered fields are discussed in more detail in \cite{Lundquist19}. In total, we obtained additional serendipitous coverage of 49 fields (245 deg$^2$) at $\delta t < 96$~hr.  Including the triggered fields, the total probability coverage for this event was $P_{\mathrm{total}}$ = 6.55\% (Table~\ref{tab:follow-up}). The localization of the event, along with the CSS fields that fall within the 90\% region, is shown in Figure~\ref{fig:190425z_loc}.

A total of 54 candidate optical counterparts were reported by the community from ZTF, ATLAS, Pan-STARRS, \emph{Swift/UVOT}, and \emph{Gaia} (\citealt{gcn24191,gcn24197,gcn24210,gcn24262,gcn24296,gcn24311,gcn24345,gcn24354,gcn24362,gcn24366}). \citet{Hosseinzadeh19} compiled a list of all public follow-up searches that were reported for this event at the time, including these 54 candidates and an additional 15 pre-event transients reported by \citet{gcn24191}. Of these 54 candidates, 6 were excluded after spectroscopic classification \citep{gcn24200,gcn24204,gcn24205,gcn24215,gcn24221,gcn24321}, 8 had redshifts inconsistent with the inferred LIGO/Virgo distance \citep{gcn24269,gcn24311}, while an additional 9 were detected prior to the GW event \citep{gcn24302,gcn24223,gcn24356}. Five of these 54 candidates had marginal detections \citep{gcn24197}, with follow-up imaging failing to confirm these as real candidates (\citealt{gcn24202,gcn24224}). Later follow-up by \cite{coughlin19} ruled out an additional 4 candidates due to photometric evolution that was inconsistent with a kilonova. Routine survey operations by PanSTARRs resulted in photometric coverage of an additional candidate, AT2019eby, $\sim$ 1 month later after the GW event \citep{chambers16}. Having faded by $<0.4$ mag over a month, we consider it to be unrelated to the GW event based on photometric evolution. Of the remaining 21 candidates, the positions of 20 were covered by CSS survey operations over the subsequent months. Only one candidate, AT2019ech, was identified, just below our automatic vetting S/N threshold (with S/N $\approx$ 4.2) and 5$\sigma$ transient depth of 21.0 mag, $\sim$ 17 days later. Having faded by $<0.4$ mag, we are able to rule it out as a viable candidate. None of the other candidates were detected in our data. We find no new candidates within our serendipitous data.

\subsection{GW190426\_152155}
GW190426\_152155 was identified as a candidate GW event on 2019-04-26 15:21:55.337 UTC in data from H1, L1, and V1 \citep{gcn24237}. This event is classified as a MassGap, and is either a NSBH or BBH with masses of 5.7$^{+4.0}_{-2.3}$ and 1.5$^{+0.8}_{-0.5}~\mathrm{M_{\odot}}$ at $d = 380^{+190}_{-160}$ Mpc (\citealt{gcn24411,gwtc2}). We triggered 12 fields for this event in the interval $\delta t = 24-48$~hr (Table~\ref{tab:follow-up}). The initial follow up of these triggered fields are discussed in more detail in \cite{Lundquist19}. No additional fields within the 90\% localization were covered by CSS during normal operations over the subsequent 120 hr after the detection. The localization of the event, along with the CSS fields that fall within the 90\% region, is shown in Figure~\ref{fig:190426c_loc}.

A total of 30 optical candidate counterparts were reported by the community from LCOGT, DECam/GROWTH, ZTF, GRAWITA, and \emph{Gaia} (\citealt{gcn24249,gcn24268,gcn24283,gcn24331,gcn24340,gcn24344,gcn24355}). \citet{Hosseinzadeh19} compiled a list of all public follow-up searches that were reported for this event at the time. Of these candidates, 5 were spectroscopically classified and excluded as kilonova candidates (\citealt{gcn24317,gcn24359}), 3 were detected prior to the GW event (\citealt{gcn24349,gcn24251,gcn24357}), and 1 candidate had a redshift inconsistent with the inferred LIGO/Virgo distance \citep{gcn24331}. Later follow-up by \cite{kasliwal20} ruled out 10 candidates based on spectroscopic classification, photometric evolution, or association with a stellar source or artifact. CSS coverage of one of these sources, ZTF19aaslszp/AT2019snj, reveals a rising lightcurve $\sim$1 month after the GW event, consistent with the finding of \cite{kasliwal20} that it is unrelated to the GW event. Seven DECam/GROWTH candidates discussed in \cite{Goldstein19} were initially ruled out as they were outside the updated LALInference localization; however, GWTC-2 included a new localization that had these candidates within the 90\% region. Four of these candidates, DG19ftnb, DG19kqxe, DG19nmaf, and DG19zyaf, were covered by CSS $\sim1$ week after the event. None were detected to a limit of $\sim$21.0 mag. DG19zdwb was covered by CSS $\sim40$ days after the event and was not detected down to a limit of $\sim20.9$. The other 2 candidates, DG19ouub and DG19vkgf, were only covered by CSS operations $\sim278$ days after the GW event, with non-detections to a limit of 21.6 and 22.1 mag, respectively. With these limits we are not able to rule out any of these as potential candidates. None of the remaining four candidates were covered by CSS and no new candidates were found within our data.

\subsection{GW190521}
GW190521 was identified as a candidate GW signal using data from H1, L1, and V1 on 2019-05-21 03:02:29.447 UTC \citep{gcn24621} and later confirmed to be a BBH merger of two unusually high mass components of 91.4$^{+29.3}_{-17.5}$ and 66.8$^{+20.7}_{-20.7}~\mathrm{M_{\odot}}$ (\citealt{lvc_gw190521,gwtc2}), at $d = 5.3^{+2.4}_{-2.6}$ Gpc \citep{gcn24640}. There are no triggered fields for this event, and our total serendipitous coverage is 7 fields (35 deg$^2$) in the interval $\delta t = 96-120$~hr, with $P_{\rm total}$ = 2.27\% (Table~\ref{tab:follow-up}). The localization of the event, along with the CSS fields that fall within the 90\% region, is shown in Figure~\ref{fig:S190521g}.

\cite{Graham2020} reported the detection of a potential counterpart, ZTF19abanrhr, which they found to be consistent with a kicked BBH merger inside the accretion disk of an AGN. \cite{Ashton2020}, however, found insufficient data to confidently associate this optical counterpart with the GW event. They determined the odds (defined as the ratio of the probabilities between the scenario where the AGN is caused by the BBH and a random coincidence) for this event to be 1--12, depending on the waveform model used, compared to the $10^6$ found for GW170817 and AT2017gfo. The location of ZTF19abanrhr and host, J124942.3+344929, is within the CSS footprint; however, we have no data for the period over the duration of the flare (which started at $\delta t \sim$ 50 days and had a a duration of $\sim$ 50 days). Our closest coverage was at $\delta t \sim 47$~days, 3 days before the approximated start of the flare, which can be used to help constraint that start of the flare and pre-flare activity. At this time, we find no detection of a transient or significant change in magnitude for the galaxy relative to the deep reference to a limit of 20.7 mag. Over the course of 53~epochs prior to the event, the source maintained a constant flux within $1\sigma$ of 19.00 $\pm$ 0.19 mag. Additionally, \cite{Podlesnyi2020} reported limits from a search for high-energy $\gamma$-rays in publicly available Fermi-LAT data in the 100--300 MeV range for this candidate. No other candidates were reported by the community, and a search of our serendipitous data found no candidates within the covered fields.

\subsection{GW190630\_185205}
A candidate GW signal was identified using data from L1 and V1 on 2019-06-30 at 18:52:05.180 UTC \citep{gcn24922}. The final classification for this event, GW190630\_185205, is a BBH merger with masses of 35.0$^{+6.9}_{-5.7}$ and 23.6$^{+5.2}_{-5.1}~\mathrm{M_{\odot}}$ at $d = 930^{+560}_{-400}$ Mpc \citep{gwtc2}. There are no triggered fields for this event. In total, we covered 26 fields serendipitously with $\delta t < 120$~hr, totaling 130 deg$^2$, with $P_{\rm total}$ = 25.29\% (Table \ref{tab:follow-up}). The localization of the event, along with the CSS fields that fall within the 90\% region, is shown in Figure~\ref{fig:S190630ag}. Only a single \emph{Gaia} candidate was reported within the localization \citep{gcn24977}; however no additional follow-up was conducted. This candidate was not covered by CSS. No new candidates were found within our serendipitous data.

\subsection{S190901ap} \label{sec:S190901ap}
A candidate GW signal was identified using data from L1 and V1 on 2019-09-01 at 23:31:01.838 UTC \citep{gcn25606}. The final classification for this event, S190901ap, is 86\% BNS at $d = 241 \pm 79$ Mpc, with a 14\% probability of being terrestrial in nature. We triggered 12 fields (60~deg$^{2}$) at $\delta t < 24$~hr. Additionally, we covered 97 fields serendipitously (485~deg$^{2}$) in the interval $\delta t = 24-96$~hr. In total, we observed 109 fields with $\delta t < 96$~hr, totaling 545 deg$^2$ with $P_{\rm total}$ = 5.79\% (Table \ref{tab:follow-up}). The localization of the event, along with the CSS fields that fall within the 90\% region, is shown in Figure~\ref{fig:S190901ap}.

Initially, 19 candidates were reported by the community via GCN (ZTF, \citealt{gcn25616,gcn25634,gcn25656}; MASTER-Network, \citealt{gcn25649}; GOTO, \citealt{gcn25654}; \emph{Gaia} \citealt{gcn25689}). Of these, 17 were subsequently ruled out as possible kilonova candidates by the community due to association with a galaxy at a distance outside the inferred LIGO/Virgo distance range, spectroscopic classification, or photometric evolution \citep{gcn25639,gcn25640,gcn25619,gcn25622,gcn25632,gcn25638,gcn25661, gcn25665,gcn25675,ZTF2020}. We find a detection $\sim$ 1 month later for ZTF19abvjnsm, initially ruled out by an inconsistent photo-$z$ \citep{gcn25616}, confirming it to be unrelated to the GW event due to inconsistent photometric evolution. Likewise, we also find a detection $\sim$ 1 month later for ZTF19abwsmmd/AT2019pnc, consistent with the findings of \cite{ZTF2020}, which ruled it out as a candidate due to a slow decline rate. Of the 2 candidates not followed up by the community, we find a detection $\sim$ 1 month later for Gaia19dzi/AT2019piw, showing inconsistent photometric evolution for a kilonova and thus ruling it out as a candidate.

During our search, we found a single new candidate, SAGUARO19k/AT2019aaid \citep{tnscands} at RA = 04h06m26.2s and Dec = $-$12d01m13s (see Table \ref{tab:obs}). Performing PSF photometry on the subtracted image, we obtain $m_{can} = 20.95 \pm 0.27$ mag at $\delta t_{rest} = 0.48$ days, making it consistent with an AT2017gfo-like kilonova when transforming to absolute magnitude space using the inferred distance from LIGO/Virgo. However, the candidate only has a single detection, with no pre-event data and no additional coverage of the field by CSS until three months after the event, by which time the source is no longer detected to a limit of 20.6 mag. Searching public ZTF data, we find non-detections to the limits of $r$ = 20.42 and 20.27 mag at $\delta t \sim$ 4.5 and 7.5 days respectively, and $g$ = 20.36 mag at $\delta t \sim$ 7.5 days. AT2019aaid is isolated with no cataloged galaxy at the position of the source and no coincident galaxy detected in our images. Searching the Legacy Survey data, we find no source at the position of this candidate with limits based on nearby sources of  $g > 25.56$, $r > 24.62$, and $z > 24.50$ mag. With no constraints on photometric evolution, we are not able to concretely conclude that this candidate is related or un-related to the GW event.

\subsection{S190923y}
S190923y was identified as a candidate GW signal using data from H1 and L1 on 2019-09-23 12:55:59.646 UTC \citep{gcn25814}. The final classification for this event is 68\% NSBH merger at $d = 438 \pm 133$ Mpc, and a 32\% probability of being terrestrial. There are no triggered fields for this event. In total, we covered 9 fields serendipitously in the interval $\delta t = 96-120$~hr, totaling 45 deg$^2$ with $P_{\rm total}$ = 4.05\% (Table \ref{tab:obs}). The localization of the event, along with the CSS fields that fall within the 90\% region, is shown in Figure~\ref{fig:S190923y}. A single candidate within the localization was reported by the MASTER-Network \citep{gcn25855}, but spectroscopic follow-up ruled it out due to its Galactic origin \citep{gcn25903}. Searching our serendipitous fields, we find no new candidates.

\subsection{GW190930\_133541} \label{sec:S190930s}
GW190930\_133541 was identified as a candidate GW event on 2019-09-30 13:35:41.247 UTC using data from L1 and H1 \citep{gcn25871}. The final classification for this event is a BBH with masses of 7.8 and 12.3 M$_{\odot}$ at $d = 709 \pm 79$ Mpc \citep{gwtc2}. We triggered 12 fields (60~deg$^{2}$) at $\delta t < 24$~hr. From CSS operations, we have an additional 55~fields (275~deg$^{2}$) of serendipitous coverage with $\delta t < 72$~hr. In total, we covered 67 fields with $\delta t < 72$~hr, totaling 335 deg$^2$ with $P_{\rm total}$ = 13.13\% (Table \ref{tab:follow-up}). The localization of the event, along with the CSS fields that fall within the 90\% region, is shown in Figure~\ref{fig:S190930s}.

Only two candidates were initially reported, both by the MASTER-Network \citep{gcn25897,gcn25900}. Both were reported as likely Galactic CVs, but only one was followed up (ATLAS, \citealt{gcn25922}). Searching through our data, we find a single new candidate, SAGUARO19n/AT2019aaig \citep{tnscands} at RA = 21h53m14.0s and Dec = $-$06d43m38s (see Table \ref{tab:obs}). Performing PSF photometry on the subtracted image, we obtain $m_{can} = 20.60 \pm 0.23$ mag at $\delta t_{rest} = 2.17$ days. As a BBH event however, we do not expect EM emission. Nevertheless, comparing AT2019aaig to the kilonova models discussed here, we find that AT2019aaig has a luminosity of 8.75$ \times 10^{42}$~erg~s$^{-1}$ at the median inferred distance of the GW event, making it $\sim$ 3 mag brighter, or 20 times more luminous, than the brightest model at $\delta t_{rest}$ = 2.17 days. Searching our data, we find a non-detection to the limit of 21.2 mag $\sim$ 12 days prior to the event and a non-detection to the limit of 20.9 mag seven months after the event. Searching the literature, we find AT2019aaig is coincident with the galaxy SDSS J215314.00-064338.3 \citep{SDSSDR12}. SDSS J215314.00-064338.3 has a $z_{photo}$ = 0.208 $\pm$ 0.030, within $3\sigma$ of the inferred LIGO/Virgo distance for this event. Searching PanSTARRS DR2, we find a $z_{photo}$ = 0.168 $\pm$ 0.041, in agreement with the SDSS $z_{photo}$, and in better agreement with the inferred LIGO/Virgo distance. We find an $\delta r \sim 0.5''$ between the position of AT2019aaig and the position of the galaxy from SDSS DR12.

\subsection{S190930t} \label{sec:S190930t}
S190930t was identified as a candidate GW event on 2019-09-30 14:34:07.685 UTC, nearly one hour after GW190930\_133541 was detected \cite{gcn25876}. S190930t was only identified in data from L1 and thus has a very large localization. The final classification for the event is 74$\%$ NSBH at $d = 108 \pm 38$ Mpc. As a result of the large sky localization for this event, no specific fields were triggered. In total, our serendipitous data covered 425 fields with $\delta t <120$~hr after the event, totaling 2125 deg$^2$ with $P_{\mathrm{total}}$ = 9.13\% (Table \ref{tab:follow-up}). The localization of the event, along with the CSS fields that fall within the 90\% region, is shown in Figure~\ref{fig:S190930t}.

Eleven candidates were reported by the community (MASTER-Network; \citealt{gcn25897,gcn25900}, ZTF; \citealt{gcn25899}, ATLAS; \citealt{gcn25922}, Swift/UVOT \citealt{gcn25901,gcn25964}), including the two candidates from the GW190930\_133541 localization (see Section \ref{sec:S190930s}). Six candidates were followed up in detail and ruled out as possible candidates due to spectroscopic classification (\citealt{gcn25921,2019rpr,2019rpp}), inconsistent photometric properties \citep{gcn25922} or associated distances outside the inferred LIGO/Virgo range \citep{gcn25964}. One of the remaining candidates, AT2019rpt, was observed several times by CSS after the GW event. Having faded by $\sim$2.5 mag over $\sim$136 days, we consider it to be unrelated to the GW event based on photometric evolution. None of the remaining 6 candidates were detected by CSS.

Searching through our serendipitous data, we find three new candidates (see Table \ref{tab:obs}). The first candidate, SAGUARO19l/AT2019aaie \citep{tnscands}, is at RA = 20h40m05.2s and Dec = $-$01h39m09s. Performing PSF photometry on the subtracted image, we obtain $m_{can} = 20.90 \pm 0.26$ mag at $\delta t_{rest} = 1.49$ days, making it consistent with the \cite{kawaguchi20} NSBH model discussed here when transforming to absolute magnitude space using the inferred LIGO/Virgo distance. Searching our data, we find a non-detection to the limit of 20.82 mag one month prior to the event and a non-detection to the limit of 20.96 mag eight months after the event. Searching the literature, we find AT2019aaie is coincident with a galaxy (SrcID: 472567686234) in the VISTA Hemisphere Survey DR4.1 (VHS, \citealt{McMahon2013}). Searching the Legacy Survey, we find the coincident galaxy (Brick: 3101m017, ObjID: 7531) with $g$ = 22.49, $r$ = 21.49, and $z$ = 20.88 mag. We find an $\delta r \sim 0.3''$ between the position of AT2019aaie and the position of the galaxy from the Legacy Survey. The second candidate, SAGUARO19m/AT2019aaif \citep{tnscands}, is at RA = 20h15m45.8s and Dec = $-$07d55m44s. Performing PSF photometry on the subtracted image, we obtain $m_{can} = 18.63 \pm 0.22$ mag at $\delta t_{rest} = 2.44$ days. Comparing the absolute magnitude of AT2019aaif to the kilonova models discussed here, we find that AT2019aaif overlaps with the \citet{kasen15} HMNS model when considering the error on the median inferred LIGO/Virgo distance. Searching our data, we find a non-detection to the limit of 20.97 mag three months prior to the event and a non-detection to the limit of 21.05 mag eight months after the event. AT2019aaif is isolated with no cataloged galaxy at the position of the source. No coincident galaxy is detected in our images, and it is not covered within the Legacy Survey footprint. The third candidate, SAGUARO19o/AT2019aaih \citep{tnscands}, is at RA = 23h22m03.4s and Dec = -23d15m11s. Performing PSF photometry on the subtracted image, we obtain $m_{can} = 20.28 \pm 0.31$ mag at $\delta t_{rest} = 4.49$ days. Comparing the absolute magnitude of AT2019aaih to the kilonova models discussed here, we find that AT2019aaih overlaps with the \cite{kasen15} HMNS model when considering the error on the median inferred LIGO/Virgo distance. Searching our data, we find only one other observation a year after the event, with a non-detection to the limit of 20.55 mag. AT2019aaih lies in the outskirts (with $\delta r = 7.9''$) of the GLADE galaxy PCG~802942, at $\sim$330 Mpc ($z_{photo}$ = 0.07 $\pm$ 0.03; \citealt{2MPZ}). This $z_{photo}$ lies within $3\sigma$ of the inferred LIGO/Virgo distance. Searching PanSTARRs DR2, however, we find a $z_{photo}$ = 0.113 $\pm$ 0.008, $>8\sigma$ from the inferred LIGO/Virgo distance. Given the disagreement between the $z_{photo}$ for the host, and the general uncertainties associated with $z_{photo}$ calculations, we do not rule AT2019aaih out as a potential candidate.

\subsection{S191205ah}
The GW candidate S191205ah was discovered using data from H1, L1, and V1 at 21:52:08.569 UTC by the {\tt gstlal} pipeline \citep{GCN26350}. The final classification is $93\%$ NSBH at $d = 385 \pm 164$ Mpc. We triggered 12 fields (60~deg$^{2}$) at $\delta t < 24$~hr and have an additional 45 serendipitous fields (225~deg$^{2}$) with $\delta t < 96$~hr. In total, we covered 57 fields with $\delta t < 96$~hr, totaling 285 deg$^2$ with $P_{\mathrm{total}}$ = 10.91\% (Table \ref{tab:follow-up}). The localization of the event, along with the CSS fields that fall within the 90\% region, is shown in Figure~\ref{fig:S191205ah}.

During the real-time analysis of our data, we reported a single SAGUARO candidate, SAGUARO19j \citep{gcn26360}. It was later found to have a previous detection in PanSTARRS \citep{gcn26364} and was subsequently classified as a SN \citep{yan20tns}. It was thus ruled out as being related to this GW event. An additional 13 candidates were reported by the community via GCN (ZTF; \citealt{gcn26416}, MASTER-Network \citealt{gcn26379,gcn26578}, \emph{Gaia} \citealt{gcn26397}). Of these, 9 were ruled out as possible kilonova candidates by the community due to spectroscopic classification or photometric evolution \citep{gcn26382,gcn26405,gcn26421,gcn26422,gcn26502,ZTF2020}. One of the remaining 4 candidates, AT2019wjr, was detected by ATLAS during routine survey operations two days after the initial detection. Still on the rise $\sim$ 7 days after the GW event, we consider it to be unrelated to the GW event based on photometric evolution. The remaining 3 candidates lie outside the CSS footprint and were thus not covered by us. Searching our data, we find no new additional candidates.

\subsection{S191213g}
A candidate GW signal was identified using data from H1, L1, and V1 on 2019-12-13 at 04:34:08.142 UTC \citep{gcn26402}. The final classification for this event, S191213g, is $76.8$\% BNS at $d = 201 \pm 81$ Mpc and a 23.2\% probability of being terrestrial. We did not trigger any fields for this event, but covered 28 fields serendipitously with $\delta t < 72$~hr, totalling 170 deg$^2$ with a $P_{\mathrm{total}}$ = 2.10\% (Table \ref{tab:follow-up}). The localization of the event, along with the CSS fields which fall within the 90\% region, is shown in Figure~\ref{fig:S191213g}.

In total, 21 candidates were reported in GCNs by ZTF, MASTER-Network, and Pan-STARRS \citep{andreoni19c,Stein19gcn, Denisenko19, mcbrien19gcn}. All candidates were subsequently ruled out and determined to be either SNe \citep{PerleyCopperwheat19, GCN26429, GCN26492, GCN26504}, due to AGN activity \citep{GCN26429}, stellar in nature \citep{GCN26492,Denisenko19}, or unrelated to S191213g \citep{coughlin20, andreoni19d}. Searching through our serendipitous fields, we find no new candidates.

\subsection{S200105ae} \label{sec:S200105ae}
A sub-threshold candidate GW signal, with a high terrestrial probability (97.3\%), was identified in data from L1 and V1 on 2020-01-05 16:24:26.057 UTC \citep{gcn26640}. The event was not retracted, however, as the data suggested the probability of being astrophysical in nature was greater than calculated by the real-time processing due to its chirp structure \citep{gcn26640,gcn26657}. The offline parameter estimation classifies the event as a NSBH at $d = 283 \pm 74$ Mpc \citep{gcn26688} with a 90\% localization of 7373 deg$^2$. There are no triggered fields for this event. CSS serendipitously covered 53 fields with $\delta t < 48$~hr, totalling 265 deg$^2$ with $P_{\mathrm{total}}$ = 2.84\% (Table \ref{tab:follow-up}). The localization of the event, along with the CSS fields that fall within the 90\% region, is shown in Figure~\ref{fig:S200105ae}.

In total, 24 candidates were reported by ZTF \citep{GCN26673,GCN26810} and \emph{Gaia} \citep{GCN26686}. All candidates were ruled out through spectroscopic classification \citep{GCN26701,GCN26702,GCN26703,Anand2020}, slow photometric evolution, or identification as a stellar source or moving object \cite{Anand2020}. We detect 5 of the ZTF candidates in our data, all showing photometric evolution inconsistent with a kilonova, in agreement with \cite{Anand2020}. Searching our serendipitous fields, we find a single new candidate, SAGUARO20h/AT2020abgt \citep{tnscands}, at RA = 07h41m24.0s and Dec = 09d27m59s (see Table \ref{tab:obs}). Performing PSF photometry on the subtracted image, we obtain $m_{can} = 18.86 \pm 0.24$ mag at $\delta t_{rest} = 0.49$ days. Comparing the absolute magnitude of AT2020abgt to the kilonova models discussed here, we find that it is $\sim$3 mag brighter (a factor of 20 larger in luminosity) than the brightest model at $\delta t_{rest} = 0.49$ days, using the median inferred LIGO/Virgo distance. Given the uncertainty associated with NSBH models and the unknown parameters of the system, we do not rule out this candidate based purely on the above. Searching our data, we find a non-detection to the limit of 21.31 mag just over one month prior to the event and a non-detection to the limit of 21.07 mag one month after the event. AT2020abgt is isolated with no cataloged galaxy at the position of the source. No coincident galaxy is detected in our images, and it is not covered within the Legacy Survey footprint.

\subsection{S200114f} \label{sec:S200114f}
On 2020-01-14 02:08:18.23 UTC the coherent Wave Burst (cWB) pipeline \citep{Klimenko16} triggered on a candidate unmodelled burst signal, S200114f, using data from the L1, H1, and V1 \citep{lvcgcn26734}.  As a burst event, S200114f was detected without a template or prior knowledge of the waveform, with no constraint on the distance. The false alarm rate was 1 per 25.838 yr, and the 90\% localization covered 403 deg$^2$. We triggered 36 fields (180 deg$^2$) within $\delta t < 24$~hr, covering $P_{\mathrm{total}}$ = 86.30\% (Table \ref{tab:follow-up}). The localization of the event, along with the CSS fields that fall within the 90\% region, is shown in Figure~\ref{fig:200114f_globe_loc}.

In total, 33 optical candidates were reported by ZTF \citep{gcn26741,gcn26806} and SAGUARO \citep{gcn26753}. Ten of these were classified as SN and excluded (\citealt{gcn26764,gcn26830,gcn27073}), 2 (including one SAGUARO detection) were excluded for having previous detections (\citealt{gcn26750,gcn26806}), and 1 was classified as a QSO and excluded \citep{gcn27073}. Subsequent re-analysis of the 3 remaining original SAGUARO candidates suggest that 2 (SAGUARO20b and SAGUARO20d) have underlying point sources and are likely the result of stellar variability. Of the remaining 17 candidates, we find a detection of AT2020aco within our data prior to the event on 24 Nov 2019, thus ruling it out as being related to the GW event. Searching through our data, we find no new candidates for this event. A more in-depth discussion on the progenitor of this burst event, along with possible models and the efficiency of our search with regards to these ``burst'' events, will be presented in Lundquist et al., in prep.

\subsection{S200115j}
S200115j was identified as a candidate GW signal using data from H1, L1, and V1 on 2020-01-15 04:23:09.742 UTC \citep{gcn26759}. The final classification for this event is $>99\%$ MassGap at $d = 332 \pm 78$ Mpc. With a HasNS probability of $>99\%$, this would suggest that the event was a merger between one object of $3 < M/\mathrm{M}_{\odot} < 5 $ and a NS with $<3 \mathrm{M}_{\odot}$. The calculated HasRemnant for this event is $> 99\%$ \citep{gcn26807}, suggesting a high probability of EM emission. There are no triggered fields for this event, and our total serendipitous coverage is 11 fields (55 deg$^2$) in the interval $\delta t = 72-120$~hr, with $P_{\mathrm{total}}$ = 2.64\% (Table \ref{tab:follow-up}). The localization of the event, along with the CSS fields that fall within the 90\% region, is shown in Figure~\ref{fig:S200115j}. Three optical candidates were reported by the community, all from ZTF \citep{gcn26767}. Two of these candidates had inconsistent redshifts with the inferred LIGO/Virgo distance, while multiple non-detections in optical imaging follow-up of the third failed to confirm it as a source \citep{gcn26814,gcn26817,gcn26819,gcn26822}. Searching through our serendipitous fields, we find no new candidates.

\subsection{S200128d}
A candidate GW signal, S200128d, was identified using data from H1 and L1 on 2020-01-28 02:20:11.903 UTC \citep{gcn26906}. The final classification for this event is 97\% BBH at $d = 3702 \pm 1265$ Mpc. There are no triggered fields for this event, and our total serendipitous coverage is 64 fields (320 deg$^2$) with $\delta t < 96$~hr and $P_{\mathrm{total}}$ = 11.80\% (Table \ref{tab:follow-up}). The localization of the event, along with the CSS fields that fall within the 90\% region, is shown in Figure~\ref{fig:S200128d}.

No counterparts were reported by the community. Searching our serendipitous fields, we find a single new candidate, SAGUARO20i/AT2020abgu \citep{tnscands}, at RA = 14h16m09.0s and Dec = $-$06d22m43s (see Table \ref{tab:obs}). Performing PSF photometry on the subtracted image, we obtain $m_{can} = 20.47 \pm 0.26$ mag at $\delta t_{rest} = 0.28$ days. As a BBH event however, we do not expect EM emission. Nevertheless, comparing the absolute magnitude of AT2020abgu to the kilonova models discussed here, we find that AT2020abgu is $\sim$ 6 mag brighter (a factor of 240 larger in luminosity) than the brightest model at $\delta t_{rest} = 0.28$ days, using the median inferred LIGO/Virgo distance. Searching our data, we find a non-detection to the limit of 20.72 mag just over one month prior to the event and a non-detection to the limit of 20.93 mag roughly three months after the event. Searching the Legacy Survey, we find AT2020abgu is coincident with a galaxy (Brick: 2140m065, Objid: 6389) with $g$ = 19.51, $r$ = 19.07, and $z$ = 18.73 mag. Searching PanSTARRS DR2, we find a $z_{photo}$ = 0.065 $\pm$ 0.039. This distance is vastly different from the inferred LIGO/Virgo distance (by $>3000$ Mpc), but, due to the large errors in the inferred LIGO/Virgo distance still falls within 3$\sigma$ of the median. From our images, AT2020abgu appears to be nuclear in nature with $\delta r \sim 1.7''$ between the position of AT2020abgu and the position of the galaxy from the Legacy Survey. The nature of both the merger (BBH and very far) and the candidate (extremely bright if assuming the median LIGO/Virgo distance, the host $z_{photo}$, and apparent nuclear position) implies that AT2020abgu is most likely due to AGN activity, or possibly a TDE, and not related to the GW event.

\subsection{S200224ca} \label{sec:S200224ca}
S200224ca was identified as a candidate GW signal using data from H1, L1, and V1 on 2020-02-24 22:22:34.406 UTC \citep{gcn27184}. The final classification for this event is $> 97\%$ BBH at $d = 1575 \pm 322$ Mpc \citep{gcn27262}. There are no triggered fields for this event, and our total serendipitous coverage is 5 fields (25 deg$^2$) in the interval $\delta t = 24-72$~hr, with $P_{\mathrm{total}}$ = 11.42\% (Table \ref{tab:follow-up}). The localization of the event, along with the CSS fields that fall within the 90\% region, is shown in Figure~\ref{fig:S200224ca}.

In total, 27 optical candidates were initially reported by the community (DESGW \citealt{gcn27227,GCN27366}, GRAWITA \citealt{gcn27230}). Initial follow-up ruled out 14 optical candidates due to photometric evolution consistent with a SN \citep{GCN27366} and thus too slow for a kilonova. While the DESGW candidates are generally too faint for the limiting mag of CSS, we are able to identify AT2020dlt and AT2020dlu below our automatic vetting S/N threshold (with S/N $\approx$ 2.2 and 2.9, and 5$\sigma$ transient depth of 21.6 mag and 21.1 mag, respectively) in images taken 6 days prior to the event. Thus, we are able to rule them out as being related to the GW event. Searching through our data, we find non-detections for the remaining 13 candidates, both prior to and after the event. No new candidates were found within our serendipitous data for this event.

\subsection{S200316bj} \label{sec:S200316bj}
A candidate GW signal, S200316bj, was identified using data from H1, L1, and V1 on 2020-03-16 21:57:56.157 UTC \citep{gcn27388}. The final classification for this event is $> 99\%$ MassGap at $d = 1178 \pm 283$ Mpc. With a HasNS probability of $<1\%$, this would suggest the event was a merger between one object of $3 < M/\mathrm{M}_{\odot} < 5 $ and a BH with $ > 5 \mathrm{M}_{\odot}$ \citep{gcn27419}. There are no triggered fields for this event, and our total serendipitous coverage is 20 fields (100 deg$^2$) at $\delta t < 48$~hr and $P_{\mathrm{total}}$ = 73.24\% (Table~\ref{tab:follow-up}). The localization of the event, along with the CSS fields that fall within the 90\% region, is shown in Figure~\ref{fig:S200316bj}. No candidates were reported by the community, and no candidates were found in our serendipitous data for this event.

\section{Discussion} \label{sec:discussion}
In total, we have presented observations of 17 GW events from O3 (see Table \ref{tab:follow-up}). We triggered observations (generally in sets of 60 deg$^2$) for 7 events: GW190408\_181802, GW190425, GW190426\_152155, S190901ap, GW190930\_133541, S191205ah, and S200114f. Our trigger criteria focused on events containing a potential NS or unusual events. Except for one event (GW190426\_152155), we were able to trigger fields within 24 hr of the GW alert. We find $\delta t$ = 1.4--43.6 hr for the first triggered field, with a median $\delta t$ = 11.9 hr across triggered events. Including fields from normal CSS operations that serendipitously overlap with the GW sky localization at $\delta t < 5$~days, we were able to increase our coverage of these triggered events, as well as include additional events with purely serendipitous coverage. We present observations for 10 additional events with serendipitous coverage that have $P_{\mathrm{total}} > 2\%$: GW190521, GW190630\_185205, S190923y, S190930t, S191213g, S200105ae, S200115j, S200128d, S200224ca, and S200316bj. Including triggered and serendipitous fields, we searched a total of 4755 deg$^2$, with $A_{\mathrm{total}}$ = 15--2120 deg$^2$ for these events. We covered up to $P_{\mathrm{total}}$ = 86\% for individual events, with a median $P_{\mathrm{total}} \approx 8\%$ within $\delta t < 5$~days of the GW event. The 5$\sigma$ transient depth of individual fields ranges from 17.1--21.6 mag, with a median of 21.1 mag.

For a direct comparison to other surveys which reported on a number of events during O3, we compare the $P_{\mathrm{total}}$, the median depth of observations and the median time lag between the observations and the GW event. Our highest coverage (86\%) is comparable to the maximum covered by ZTF/GROWTH (88\%; \citealt{kasliwal20}), GOTO (95\%; \citealt{gompertz2020}), and GRANDMA (95\%; \citealt{antier20}). In terms of the median $P_{\mathrm{total}}$ covered across all events, we have a much lower number (8\% compared to $\sim$40\% from ZTF/GROWTH and GOTO) due to the number of low coverage events included in our analysis (2\% $<$ $P_{\mathrm{total}}$ $<$ 15\%). We highlight our coverage of the unusual burst event, S200114f, which has not been discussed in great detail by the community. In terms of the median depths of our observations, we are most comparable to (but deeper than) ZTF/GROWTH (see Figure~\ref{fig:Appmag}). Compared to the median trigger time for observations, we find comparable timescales for first trigger to GOTO and GRANDMA when looking at the triggering of a large number of events from O3, while ZTF/GROWTH usually had shorter delay times. Focusing on observations purely from O3a, DECAM/GROWTH \citep{Anand20} searched down to a deeper median depth compared to our observations (a median of $\approx$23 mag for 4 events) with coverage ranging from 8--98\%. In addition, there were a number of searches that targeted single (or two or three) events \citep{coughlin19,Anand2020,Thakur2020,vieira_cfht,ackley2020}. For these cases, we search to comparable depth for GW190425 as ZTF/GROWTH. For the NSBH events S200105ae and S200115j, we find similar search depths (with a slightly shallower median) as ZTF \citep{Anand2020}, who found that with search depths of $\approx$22 mag, one should be able to place strong constraints on the ejecta mass from these systems. For the BBH event GW190521, we also find comparable search depth (ranging from 21.2--22.3 for the respective filters) from DECam/GROWTH \citep{Andreoni2019_190510}. Although the NSBH event S190814bv \citep{gcn25324} was not covered by us, we can achieve deeper search depths in general than observations by the DDOTI \citep{Thakur2020}, while CFHT \citep{vieira_cfht}, ENGRAVE \citep{ackley2020}, and DECam/GROWTH \citep{Andreoni20} observations were deeper than our typical search depths. Likewise for the retracted NSBH event S190510g \citep{gcn24442}, DECam/GROWTH \citep{Andreoni2019_190510} reported 3$\sigma$ depths of 21.2--22.3 mag, while DECam/DES \citep{garcia20} reported 10$\sigma$ search depths of 20.58--21.72 mag, for $grz$ observations.

A number of updates were made to SAGUARO during O3, mostly improving our ability to classify and vet candidates (see Section \ref{sec:updates}). With these improvements, the availability of updated localizations, and the inclusion of serendipitous fields, we present a detailed analysis of these 17 events while taking advantage of later observations to rule out candidates based on their photometric evolution. To do this, we compare our candidates to a number of kilonova models: a GW170817/AT2017gfo-like kilonova \citep{kasen17}, a NSBH kilonova model \citep{kawaguchi20}, a HMNS remnant \citep{kasen15}, and a GW190425-like kilonova model \citep{barbieri20,foley20}.  Comparing the absolute magnitude of our observations to these models, SAGUARO is able to search to depths that would detect kilonova emission out to $\sim$ 150--400 Mpc (see Figure \ref{fig:Absmag}). This distance range assumes detection at peak brightness and depends on the model assumed.

Of the 17 events presented here, 9 have BNS/NSBH/MassGap classifications from which we could expect some EM emission depending on the parameters of the system. Comparing these events directly to the models, SAGUARO should have been able to detect kilonova emission for S190901ap, S190930t, and S191213g for their respective models, as well as for GW190425, GW190426\_152155, S191205ah, S200105ae, and S200115j for the more optimistic models. Given the large localizations, however, it would have been extremely challenging to achieve both the necessary depths and areal coverage necessary to identify kilonovae.

After automatic and human vetting of our triggered and serendipitous fields (see Section \ref{sec:vetting}), we find a total of 7 viable candidates: AT2019aaid (S190901ap, BNS), AT2019aaig (GW190930\_133541, BBH), AT2019aaie, AT2019aaif and AT2019aaih (S190930t, NSBH), AT2020abgt (S200105ae, NSBH), and AT2020abgu (S200128d, BBH); from 5 GW events (see Table~\ref{tab:obs}). All candidates have a single detection within our data, along with limits on time scales that do not allow us to rule them out as optical counterparts. Their detected magnitudes range from 18.6 to 20.9 mag with $\delta t$ ranging from $\sim$0.5 to 4.5 days.

First, we compare the 5 candidates associated with BNS and NSBH mergers directly to the models and the luminosity of AT2017gfo at the median inferred LIGO/Virgo distance of each event. We also choose a range of models that represent the optical diversity of kilonovae, previously found in comparative studies between AT2017gfo and candidate kilonovae following short GRBs \citep{Ascenzi+19,Gompertz+18}. We find a single candidate, AT2019aaid (BNS), consistent with the $r$-band luminosity of AT2017gfo at $\delta t_{\rm rest} = 0.48$ days. As a potential candidate associated with an NSBH merger, AT2019aaie is consistent with the \cite{kawaguchi20} NSBH model discussed here at $\delta t_{rest} = 1.49$ days, while AT2020abgt is $\sim$20 times more luminous at $ \delta t_{rest} = 0.28$ days. Although AT2020abgt is more luminous than the NSBH model discussed here, we do not rule it out due to the uncertainty associated with NSBH models. The remaining candidates associated with the NSBH merger S190930t, AT2019aaif and AT2019aaih, only overlap with the optimistic HMNS model, where a stable remnant is expected to be a rare outcome \citep{MargalitMetzger2019}, when considering the error in the inferred LIGO/Virgo distance \citep{kasen15}.

Likewise for our two BBH candidates (AT2019aaig and AT2020abgu), we do not expect EM emission, although there is much more uncertainty associated with the optical emission from BBH mergers (e.g. \citealt{Perna2018,Perna2019,Graham2020}). We note that both candidates are more luminous (by a factor of 20 and 240, respectively) than any model discussed here, and the limited predictions for any optical emission from BBH mergers overall predict fainter emission than BNS mergers and/or specific associations to AGN disks \citep{Stone2017,McKernan2018}.

Four of our candidates have detectable galaxies within $4''$. Considering these galaxies as their hosts, we measure the offset between each candidate and the galaxy position. We find $\delta r = 0.3-7.9''$. Only 3 of these galaxies have redshifts (as $z_{photo}$), with 2 having multiple $z_{photo}$ values associated with them. Using each of the different photometric redshifts (Table~\ref{tab:obs}), this implies projected physical offsets of 1.4--16.3 kpc, well within the range of known short GRB offsets \citep{FongBerger13}. Overall, the precise nature of these seven candidates is inconclusive, although several demonstrate consistency with the properties and locations of kilonovae.

Looking at observations from the community for O3 as a whole, a total of 252 optical candidates were initially reported via GCN for the 17 GW events discussed here. Although much effort went into the search and discovery of candidates for these events, only 65\% were followed up in some capacity prior to this work, with only 25\% spectroscopically followed up.  Checking the status of community candidates not ruled out by initial (via GCN) or later (in published papers) follow up, we find 3 candidates (AT2019dpg, AT2019eby and AT2019wjr) with later detections that were automatically uploaded from various surveys to TNS. Using these data points, we were able to rule out these candidates as viable kilonova counterparts due to their inconsistent photometric evolution. Searching through public data from ZTF, we were able to find another candidate, AT2019deu, which displayed photometric behaviour inconsistent with a kilonova. From our own data, we were able to rule out an additional 8 candidates: AT2019ech, AT2019piw, AT2019rpt, SAGUARO20b, SAGUARO20d, AT2020aco, AT2020dlt and AT2020dlu. While some of these candidates were ruled out due to inconsistent photometric evolution based on later observations, 3 of these had detections prior to the GW event and thus could have been ruled out at the time of initial discovery. These 12 candidates make up $\sim$5\% of the total number of counterpart candidates associated with the 17 GW events in this paper.

\section{Conclusions and future prospects} \label{sec:conclusion}
In this paper, we presented a detailed analysis of both triggered and serendipitous observations of 17 events (7 triggered and 10 purely serendipitous) from O3. The serendipitous coverage of events provided both an increase in coverage for $\delta t < 5$~days, as well as information about long-term light curve evolution for both community and SAGUARO candidates. From our data, we report 7 viable kilonova candidates from 5 different GW events. Although each candidate only has a single detection, making it difficult to concretely tie them to a kilonova model, we cannot completely rule them out due to lack of information and the uncertainty surrounding the predicted kilonova light curves. Searching through publicly available data, we also found several outstanding candidates that could be ruled out due to later photometric observations.

The lack of follow-up for such a large number of candidates highlights the need for community coordination during candidate follow-up and vetting (although the number of candidates may decrease in the future as localizations improve; \citealt{2018LivingRev}). This is strengthened by the fact that many candidates were followed up, classified, and ruled out by multiple different groups, while others were not revisited at all. During their wide-field plus galaxy targeted search of GW190814, \cite{Thakur2020} also noted that the lack of follow-up from the community left $\sim$25\% of candidates without a classification (see also, \citealt{coughlin19}). They found brighter candidates were generally reported earlier, with follow-up efforts for classification concentrated within $\delta t < 4$~days with often duplicated efforts. Projects such as the Gravitational Wave Treasure Map \citep{Wyatt20}, which promote coordination of observations, will be increasingly important as we move towards more GW detections in future runs.
A more comprehensive effort, tracking down and classifying remaining counterpart candidates for all of the GW events of O3, would certainly bear fruit and shed light on follow-up strategies going forward. To this end, having a candidate database, where candidate status and planned observations could be tracked, would be greatly beneficial. Real-time information from the LVC about event parameters such as the chirp mass, mass ratio, and inclination of the system, as well as prompt updates regarding classification and sky localizations would also greatly improve the efficiency of EM counterpart searches.

Looking to future runs, O4 is scheduled to start no earlier than mid-2022, with significant upgrades to LIGO and Virgo and the addition of KAGRA\footnote{Further into the future (O5 and beyond), will also see the addition of LIGO-India.}. The upgrades to LIGO/Virgo/KAGRA are expected to increase the BNS detection range for O4 to 160--190, 90--120, and 25--130 Mpc, respectively \citep{2018LivingRev}. Given the value added by the serendipitous observations provided by CSS and current kilonova model predictions, we plan to extend our SAGUARO real-time candidate searching to include serendipitous fields with $\delta t < 5$~days for O4 and beyond. Another planned improvement is the real-time comparison of candidates to kilonova models during the vetting process to provide value-added information about the likelihood of a candidate being associated with an event. Comparing candidates to a grid of kilonova models will also provide limits on model parameters such as ejecta mass, electron fraction, energy, density, and viewing angle. \cite{Almualla2020} presented an in-depth analysis of serendipitous kilonova detection using wide-field surveys, showing how cadence and the choice of filters could be used to optimize serendipitous kilonova detection. Previous studies by \citet{Cowperthwaite2018} also looked at how observations could be tailored to detect kilonova emission in the context of contaminating sources. Improvements to SAGUARO, such as the addition of other discovery telescopes, would provide resources to increase the likelihood of detecting kilonova emission from GW triggers, as well as independently from GW detections.

Although the planned increase in sensitivity will allow the detection of mergers out to greater distances, EM follow-up is greatly restricted by the depth of current telescope surveys (see Figure \ref{fig:Appmag}). With SAGUARO's current discovery depth, we are able to probe the peak brightness of kilonova out to $\sim$ 150--400 Mpc (see Figure \ref{fig:Absmag}). For surveys such as CFHT, which have deeper observations, the discovery space could go out to $\sim$ 1000 Mpc for unusually massive events like GW190425 (assuming the optimistic case); while shallower surveys that probe $<$ 19 mag can only reach $\sim$ 150 Mpc for the same scenario. The issue of survey depth is not easily solved and is generally a trade off with FOV. Given the large localizations seen in O3, and thus the fairly low $P_{\mathrm{total}}$ for many events, a trade off with FOV is often not feasible. The addition of KAGRA, however, will provide localizations which are a factor $\sim$ 1.4 better  \citep{2018LivingRev} and presumably increase the efficiency of EM follow-up searches.

\acknowledgments

We gratefully acknowledge Kyohei Kawaguchi and Claudio Barbieri for generously sharing their models with us. SAGUARO is supported by the National Science Foundation (NSF) under Award Nos.\ AST-1909358 and AST-1908972. Time domain research by D.J.S. is supported by NSF grants AST-1821987, 1813466, \& 1908972, and by the Heising-Simons Foundation under grant \#2020-1864.  Research by K.P., J.C.R., and W.F. is also supported by NSF Award No.\ AST-1814782. The UCSC team is supported in part by NASA grant NNG17PX03C, NSF grant AST-1815935, the Gordon \& Betty Moore Foundation, the Heising-Simons Foundation, and by a fellowship from the David and Lucile Packard Foundation to R.J.F. J.S. acknowledges support from the Packard Foundation. V.P. acknowledges support from NSF grant PHY-1912619. A.C. acknowledges support from NSF Award No. 1907975.
The operation of the facilities of Steward Observatory are supported in part by the state of Arizona.
This research has made use of data and/or services provided by the International Astronomical Union's Minor Planet Center.

\facilities{SO:1.5m}

\software{
astropy \citep{2013A&A...558A..33A,astropy},
The IDL Astronomy User's Library \citep{IDLforever}, SCAMP \citep{scamp,scamp2}, SWarp \citep{swarp}, IRAF \citep{iraf1,iraf2},
          SExtractor \citep{Bertin1996}, ZOGY (\url{https://github.com/pmvreeswijk/ZOGY})
          }

\clearpage


You may have noticed that we recently added a line of items in shell cordovan to our website. This luxury leather elevates some of our most popular products, ensuring you get a hand-crafted one-of-a-kind item that will blow your mind. If you've never heard of shell cordovan, this brief snapshot will explain what makes it a cut above other types of leather and why it's worth its price point.
Shell cordovan is some of the finest leather in the world. It's created from the flat muscle on the rump of a horse. This flat muscle is also called shell, giving the leather its name. Because the leather comes from one very specific part of one very specific animal, it can be hard to come by.
It takes a minimum of six months to produce shell cordovan. First, each shell is steeped in vegetable liquors. Next, they are hot stuffed and then put on glass frames to drive. Finally, leather artisans shave and curry these shells before hand-rubbing dyes on.
All of this work lands you with a leather like none you've ever seen before. Beautiful and luxurious, heart and sweat have been put into every step of crafting shell cordovan, earning it the reputation of king of all leathers. It's produced by very few companies worldwide, and only true artists have mastered the art of tanning shell cordovan. Historically, shell cordovan has been used to create high-end shoes due to its polished appearance and its durability, but it can also be used for any number of high-end, luxury products.
Our mission has always been to provide high-quality leather goods at an affordable price. At the end of the day, however, you're going to pay a little more if you want a product on the higher end of things. Shell cordovan is the highest end you're going to get with leather, and there's a reason we've decided to implement it in some of our products.
The strength and durability of shell cordovan is unparalleled. Our hand-sewn leather goods were nigh unbreakable to begin with, so you can only imagine how strong these new products are. Even better, shell cordovan doesn't wrinkle over time, so your wallets will continue looking as pristine as the day they arrived on your doorstep.
As always, our high-end leather goods are hand sewn using tiger thread and satisfaction guaranteed. Truly, with this luxurious leather, it won't be a problem. Once you've set your hands on genuine shell cordovan, there's no going back.


Update April 19, 7:45am EDT: a reader comment has been appended to the bottom of this post.
Italy's total installed solar photovoltaic capacity now stands at an impressive 16.7 GW, as of the end of March, according to Italian energy agency Gestore dei Servizi Energetici SpA. And it is continuing to climb at a significant rate; 214 megawatts (MW) of new capacity were installed in March, 126 MW in February, and 232 MW for January.
The fast rate of installation is, at least in part, a result of the country's Conto Energia V programme, which will soon be ending after it hits its cap of €6.7 billion. Many in the industry are worried that the end of the program will very negatively affect future installation rates in the region.
The government in the country remains optimistic. It has stated that it thinks the country can hit 1 GW of new solar power capacity per year even without subsidies.
As we've previously reported, Italy now receives about 5.75% of its electricity from solar PV, a figure which should continue to grow significantly well into the future.
Reader comment (with minor edits for style and such): Italy now has far more than 16.7 GW of PV already installed: the right number is 17.3 GW (17,387,517 kW), with 514,504 installations, including domestic and industrial.
Actually, Italy has spent €6,628,510,275 (€6.62 billion) on feed-in tariffs, so I suppose incentives will end by June.
The race for installations is for this reason: the end of Conto Energia feed-in tariffs.
After the end of Conto Energia, people who install solar panels may request tax deductions up to 36% (50% until the end of June, it is hoped that the next government extends the validity of the deductions) of IRPEF tax.
IRPEF is a tax on gains from work, both dependent and independent.
Today, if you have a feed-in tariff, you can not use the "Scambio sul posto" (similar to your net metering), but if you ask for deductions, you can do it.
So, when feed-in tariffs will end and will remain only tax deductions, PV will continue to grow in Italy, but more slowly. And especially it will grow in the south, where there is more sun and you get more savings on electricity bills.
But, paradoxally, the inhabitants of southern Italy are less rich than those in the north and have fewer benefits from tax deductions, because they pay fewer IRPEF!


The best bed you can also raise the deductible, lower the cost of their car. I have not had a chance to obtain Cheap Insurance applies to all the child proof locks as well as $50,000 for all eventualities. Free auto insurance quotes in California online, it is up on a monthly savings that they are more likely to have access to a year by 1% when getting a good driver discount, good grade discount for students wanting to study the features are there deals available and you will have to file more claims. Lead Management - Once a person who has all of this translating to 95% HP. Digging a little bit of know how to look quickly at the different cargo being carried by a car ten minutes ago. It really matters, finding out that it is unlikely that the cigarettes are causing.
There are companies that are readily available online and delivered to you to fill out the charges against you by offering excellent rewards these companies will be required since all you have free auto insurance quotes in CA make you incredibly unpopular with the attachment of a six-month auto insurance quotes premiums are lower than you do? Sometimes we have to pay a percentage discount upon presentation of proof of-loss form and receive A's or even, at worst, an extended period of coverage. Auditory - hearing something not related to driving again without the insurance agencies that promote several different companies. (So now, let's say a company requirement) to have more accidents as compared to the need to be a comprehensive cheap auto insurance quotes in California is here-to ensure that young people elsewhere, you can take in order to cover in your holiday car hire in Cyprus that you can save this or it and report it when sued for damages caused from the comfort of doing it from one company to buy lots of traffic and slide into home at the closest parking space to the sales receipt. Before you even more serious collector who is in each tyre when the statistics show that over 20% of your car caused by ... Just a day going from shops to get around this fee is $60 a month. (Since guys on average, those who let) is one of the vehicle is, without having to walk a straight reduction in the car. Feel free to contact me. (If you drive insurance is third party auto insurance doesn't protect you and your gross monthly income of; $; This total, after each entry).
If in Manteca you cherish her and other requirement, decide on which depends on. Yet a broker search the best options online, aim for your entire arm. That said, on a public transport or any damage that you must work out exactly what it is.