Variability-selected intermediate mass black hole candidates in dwarf galaxies from ZTF and WISE

While it is difficult to observe the first black hole seeds in the early Universe, we can study intermediate mass black holes (IMBHs) in local dwarf galaxies for clues about their origins. In this paper we present a sample of variability--selected AGN in dwarf galaxies using optical photometry from the Zwicky Transient Facility (ZTF) and forward--modeled mid--IR photometry of time--resolved Wide--field Infrared Survey Explorer ({\it WISE}) coadded images. We found that 44 out of 25,714 dwarf galaxies had optically variable AGN candidates, and 148 out of 79,879 dwarf galaxies had mid--IR variable AGN candidates, corresponding to active fractions of $0.17\pm0.03$\% and $0.19\pm0.02$\% respectively. We found that spectroscopic approaches to AGN identification would have missed 81\% of our ZTF IMBH candidates and 69\% of our {\it WISE} IMBH candidates. Only $9$ candidates have been detected previously in radio, X-ray, and variability searches for dwarf galaxy AGN. The ZTF and {\it WISE} dwarf galaxy AGN with broad Balmer lines have virial masses down to $10^{5.5}M_\odot$ and for the rest of the sample, BH masses predicted from host galaxy mass range between $10^{5.2}M_\odot<M_{\text{BH}}<10^{7.3}M_\odot$. We found that only 5 of 152 previously reported variability--selected AGN candidates from the Palomar Transient Factory in common with our parent sample were variable in ZTF. We also determined a nuclear supernova fraction of $0.05\pm0.01$\% year$^{-1}$ for dwarf galaxies in ZTF. Our ZTF and {\it WISE} IMBH candidates show the promise of variability searches for the discovery of otherwise hidden low mass AGN.


INTRODUCTION
It is challenging to determine how the very first massive black holes formed because they have grown over time as galaxies merge and their black holes accrete (Volonteri & Begelman 2010). High redshift black hole (BH) seeds are very difficult to detect due to their low luminosities (Volonteri & Reines 2016) so we can instead constrain models of BH seed formation by studying the least massive black holes in local galaxies (Reines & Comastri 2016). These low mass analogs to black hole seeds are called intermediate mass black holes (IMBHs) and have masses in the range 100 < M BH < 10 6 M . Dwarf galaxies of stellar mass M * < 10 9.75 M are the best place to find low mass black holes because galaxy mass and black hole mass are correlated (e.g. Kormendy & Ho 2013;Woo et al. 2013). Supernova-driven stunting of black hole growth in dwarf galaxies also makes these black holes comparable to early BH seeds (Habouzit et al. 2017;Anglés-Alcázar et al. 2017).
The low mass end of galaxy population relations such as the black hole occupation fraction and the slope and scatter of the relationship between black hole mass and the stellar velocity dispersion of the host bulge (M BH − σ * , Ferrarese & Merritt (2000)) will vary depending on the formation mechanism of early black hole seeds (see Greene et al. 2020, for a review). For example, Pop III stars will produce a population of undermassive BHs in low redshift galaxies with stellar dispersion σ * < 100 km/s while direct collapse mechanisms will produce heavier black holes, resulting in a flattening of the M BH − σ relation around masses of 10 5 M (Volonteri & Natarajan 2009). Although the potential for low mass black hole populations to constrain BH seed formation histories may be limited by uncertainties in accretion prescriptions (Ricarte & Natarajan 2018) it nonetheless motivates the discovery of a large population of black holes in low mass galaxies. A recent effort by Baldassare et al. (2020a) to fill in the low mass end of the M BH −σ * relation doubled the number of black holes with measured virial masses and stellar velocity dispersions in dwarf galaxies of M * < 3 × 10 9 M to 15 and the results were in agreement with an extrapolation of the linear relationship observed for high mass galaxies. However, larger numbers of low mass AGN, particularly in galaxies of mass M * < 3 × 10 8 , are needed to more fully constrain seed models.
The fraction of IMBHs in dwarf galaxies which are wandering in their galaxy haloes, rather than occupying the nucleus, may also constrain BH seed formation mechanisms. The wandering fraction is dependent on galaxy merger history but will be substantially higher if massive black holes were produced by gravitational runaway of massive star remnants (Miller & Hamilton 2002;Portegies Zwart & McMillan 2002) due to the high frequency of IMBH ejection during these interactions (Volonteri & Perna 2005;Holley-Bockelmann et al. 2008). Ricarte et al. (2021) found that a substantial population of wandering black holes exist in the ROMULUS cosmological simulations and make up 10% of black hole mass in the local universe. Studies of IMBHs of mass 3.8 < log M BH (M ) < 7.0 in cosmological zoom-in simulations suggest that 50% of IMBHs in dwarf galaxies are wandering within 7 kpc of the galaxy center due to historical galaxy mergers (Bellovary et al. 2018(Bellovary et al. , 2021. This is supported by radio observations of dwarf galaxy AGN (Reines et al. 2020). Observational constraints on the wandering fraction of IMBHs in dwarf galaxies will help to test the accuracy of cosmological merger simulations, the effects of gravitational wave recoil on black holes in low mass galaxies and the feasibility of the gravitational runaway black hole seed formation mechanism.
A major challenge for both the search for IMBHs and for estimating the occupation and wandering fractions is that the predicted accretion rates are very low, particularly for non-nuclear AGN. Black hole accretion may be bimodal in its efficiency, causing low mass black holes < 10 5 M to grow more slowly (Pacucci et al. 2018). Bellovary et al. (2018) found that most simulated IMBHs reach a maximum bolometric luminosity of log L bol (erg/s) < 41 and are therefore very difficult to detect. Pacucci et al. (2021) used a theoretical model based on galaxy mass and the angular momentum available in nuclear regions to estimate an active fraction of 5-22% of AGN in dwarf galaxies which increases with host mass. The level of activity likely depends on the merger history of the dwarf galaxy. Kristensen et al. (2021) find that inactive dwarf galaxies in the Illustris simulations tend to have been residing in dense environments for long times, while active galaxies had commonly been in a recent (≤ 4 Gyr) minor merger.
Discovery of AGN activity in dwarf galaxies will also help us to test predictions on the importance of AGN feedback in the low mass regime. While Geha et al. (2012) found that dwarf galaxies in SDSS with no active star formation are extremely rare and that more massive neighboring galaxies were the cause of star formation quenching, Penny et al. (2018) and Dickey et al. (2019) have found evidence for internal AGN-driven quenching in a small number of dwarf galaxies. Simulations suggest that dwarf galaxies hosting overmassive black holes have lower central stellar mass density, lower H I gas content and lower star formation rates than dwarf galaxies with undermassive counterparts, suggesting that internal feedback in dwarf galaxies can quench star formation only for higher mass black holes (Sharma et al. 2019). Simulations by Koudmani et al. (2019) found that AGN outflows in dwarf galaxies only have a small effect on regulating global star formation rates compared to supernovae and sustained high-luminosity AGN with isotropic winds.
Detailed observational studies of small dwarf galaxy samples have painted a different picture. GMOS IFU observations of 8 dwarf galaxies by Liu et al. (2020) discovered the presence of high velocity outflows in 7 out of 8, with some outflows capable of expelling a portion of outflowing material from the galaxy and enriching the surrounding circumgalactic medium. Coronal line emission inconsistent with shocks was also detected from 5 of these objects Bohn et al. (2021). Larger samples of AGN in dwarf galaxies will allow for further searches for signatures of AGN feedback such as high velocity and large scale outflows.
A range of approaches have been used to obtain samples of AGN candidates in dwarf galaxies. Some studies have taken a spectroscopic approach by looking for the emission line signatures of AGN in low mass galaxies. Reines et al. (2013) found that 151 dwarf galaxies with masses 10 8.5 M < M * < 10 9.5 M amongst a sample of 25,000 had [O iii] λ5007/Hβ and [N ii] λ6550/Hα narrow emission line ratios indicative of AGN activity. Mezcua & Domínguez Sánchez (2020) found a sample of 37 dwarf galaxies with 'hidden' AGN ionization lines by spatially resolving emission from star forming regions and nuclear AGN within each dwarf galaxy using MaNGA integral field unit spectroscopy. Molina et al. (2021) used [Fe X] λ6374 coronal line emission as a signature of AGN accretion to find 81 dwarf galaxies with possible IMBHs.
Multi-wavelength approaches have also been successful in the identification of AGN in dwarf galaxies. Mezcua et al. (2018) identified 40 AGN in dwarf galaxies with stellar masses 10 7 M < M * < 3 × 10 9 M via their X-ray emission in the Chandra COSMOS-Legacy Survey, finding an AGN fraction of ∼0.4% for redshifts less than 0.3. Reines et al. (2020) found that 39 of 111 dwarf galaxies had compact radio sources at 8-12 GHz with the Karl G. Jansky Very Large Array (VLA) and determined that 13 of these could confidently be classified as AGN. They also found that 10 of the 13 radio AGN detected were spatially offset from their optical galaxy nuclei.
A new approach to the detection of IMBH candidates has been to search for AGN-like stochastic variability from low mass galaxies as a signature of the presence of a central BH. Baldassare et al. (2018) found 135 galaxies with AGN-like variability on yearly timescales when constructing light curves of ∼ 28000 galaxies with SDSS spectra, including 12 from dwarf galaxies with stellar masses M * < 3 × 10 9 M . They therefore estimated a variability fraction of 0.1% for AGN in dwarf galaxies. A similar study in 2020 used light curves from the Palomar Transient Factory (Law et al. 2009) to search for variable AGN in 35,000 galaxies with stellar mass M * < 10 10 M , identifying variability in 102 galaxies with masses M * < 3×10 9 M (Baldassare et al. 2020b).
As some low mass AGN vary on hourly timescales, very high cadence surveys have proved effective for the discovery of variability from dwarf galaxy AGN. For example, Burke et al. (2020) produced a 30 minute cadence, one month long light curve of the AGN in the archetypal dwarf galaxy NGC 4395 using data from the Transiting Exoplanet Survey Satellite (TESS). The ∼ 10 5 M black hole was variable with a characteristic timescale of 1.4 +1.9 −0.5 days. Martínez-Palomera et al. (2020) used data from the HiTS imaging survey, which undertook one week of high cadence (4 times per night) and high coverage (120-150 deg 2 ) observations with the Dark Energy Camera each year for three years, to search for short timescale IMBH variability. They identified 500 galaxies with hourly, small amplitude variation in their week-long light curves. They estimated that 4% of dwarf galaxies contained an IMBH based on their results. Shaya et al. (2015) monitored ∼500 galaxies with the Kepler telescope (Borucki et al. 2010) and found that while 4% showed bright AGN activity, many other galaxies exhibited faint (down to 0.001 magnitude) variability which may also have been due to the presence of a low mass AGN.
Recently, Secrest & Satyapal (2020) undertook a search for variable AGN in the mid-IR. This was motivated by the sensitivity of mid-IR studies to low luminosity AGN which are frequently optically obscured and Compton-thick (Ricci et al. 2016;Annuar et al. 2017) and frequently unobservable in the soft X-rays (Polimera et al. 2018), along with the low contamination rates of supernovae due to weak mid-IR emission (Smitka 2016). To produce mid-IR light curves, Secrest & Satyapal (2020) used multi-epoch photometry from the Wide-field Infrared Survey Explorer (Wright et al. 2010, WISE ). WISE mapped the sky in the W1 (3.4µm) and W2 (4.6µm) bands with 6 month cadence over an 8 year baseline between the initial observations in 2010, the Near-Earth Object Wide-field Infrared Survey Explorer mission from 2010-2011 (NE-OWISE ; Mainzer et al. 2011) and the reactivation mission beginning in 2013 (NEOWISE -R; Mainzer et al. 2014). Secrest & Satyapal (2020) found a sample of 2,199 dwarf galaxies of stellar mass M * < 2 × 10 9 M and redshift 0.02 < z < 0.03 from the NASA-Sloan Atlas which had corresponding sources in the AllWISE catalog. Amongst this sample, only 2 (0.09%) showed significant variability in light curves produced by the AllWISE Multiepoch Photometry Table and the NEO-WISE -R Single Exposure Source Table. Variability-based search strategies have been particularly good at finding AGN candidates in dwarf galaxies which are optically obscured or unidentifiable by their spectroscopic signatures. For example, only 25% of the optically variable AGN candidates in galaxies of mass M * < 10 10 M found in PTF were classified as AGN or Composite galaxies based on their narrow emission lines (Baldassare et al. 2020b). This is likely due to a combination of star formation dilution of AGN emission lines and the hardening of the accretion disk SED around lower mass BHs which extends the partially ionized zone and reduces the emission line ratios normally used for AGN classification (Cann et al. 2019). Variability-based strategies therefore have an important place for finding AGN which would otherwise be missed due to biases in other selection techniques.
Previous searches for active IMBH candidates in dwarf galaxies with spectroscopic, radio, X-ray and mid-IR observations have not provided a complete, unbiased census of the IMBH occupation fraction because different strategies are hindered by high star formation rates, obscuration and low accretion luminosities. Therefore, despite the discoveries of these large samples of dwarf AGN candidates, there are very few confirmed IMBHs with well-sampled SEDs and measured virial black hole masses which we can use to occupy the sparsely populated low mass ends of a number of black hole-galaxy scaling relations. This motivates the development of effective search strategies for identifying substantially larger samples of IMBH candidates. This will provide more opportunities for careful confirmation and multiwavelength characterization of the best candidates, especially those with broad Balmer lines.
In this paper we present a comprehensive search for AGN-like variability from a large sample of dwarf galaxies in the optical and mid-IR in order to build upon the previous successes of Baldassare et al. (2018Baldassare et al. ( , 2020b and Secrest & Satyapal (2020) with SDSS, PTF and WISE. The Northern Sky Survey of ZTF Phase I (Mar 2018-Sep 2020) had an average cadence of 3 days and this was supplemented by higher cadence sub-surveys with hourly to 1 day cadences (Bellm et al. 2019b). The ongoing ZTF Phase II (Oct 2020-present) Northern Sky Survey has a 2 day cadence. By comparison, PTF had a footprint of ∼ 8000 deg 2 and a ∼ 5 day cadence. To expand upon the mid-IR variability search of Secrest & Satyapal (2020), we have made use of new forward modeled photometry catalogs of time-resolved WISE coadds to produce more sensitive photometry of a larger dwarf galaxy sample. We do this work in preparation for the upcoming Legacy Survey of Space and Time (LSST) at Vera C. Rubin Observatory (Ivezić et al. 2019) which will provide a more complete census of variable IMBHs over the next decade due to its expected single visit limiting magnitude of g ∼25 at a 3 day cadence spanning ∼10 years.
In Section 2 we describe our dwarf galaxy sample selection process. In Section 3 we describe our procedure for ZTF forced photometry of the dwarf galaxy sample.
In Section 4 we describe our selection strategy for finding the optically variable AGN. In Section 5 we present our sample of IMBH candidates and supernovae from ZTF and describe their multiwavelength and spectroscopic properties. In Section 6 we describe our selection of variable AGN in WISE based on forward modeled photometry of time resolved coadds and in Section 7 we discuss the multi-wavelength and spectroscopic properties of the WISE -selected IMBH candidates. In Section 8 we discuss the merits of the two selection strategies and the properties of the IMBH candidates in further depth.

DWARF GALAXY SAMPLE SELECTION
We selected a sample of dwarf galaxies using the NASA-Sloan Atlas (NSA) version v 1 0 1 1 . The NSA produced stellar mass estimates for galaxies by fitting 5 templates derived from stellar population synthesis models (Bruzual & Charlot 2003;Blanton & Roweis 2007) to SDSS images of galaxies after subtracting the sky background (Blanton et al. 2011). We used this catalog of stellar masses to allow for more direct comparison to other AGN variability searches which classify dwarf galaxies using stellar masses from this catalog (e.g. Reines et al. 2013;Baldassare et al. 2018Baldassare et al. , 2020bSecrest & Satyapal 2020).
When compiling our list of dwarf galaxies from the NSA we required redshifts of 0.02 < z < 0.35 and elliptical Petrosian masses of M * < 3 × h −2 10 9 M . We selected an h value of 0.73 for consistency with Reines et al. (2013) and Secrest & Satyapal (2020) such that our mass cutoff corresponds to M * < 10 9.75 M . After finding some high redshift quasars listed in the catalog with erroneous redshifts and underestimated stellar masses, we required the NSA redshift to be derived from an SDSS spectrum (described by the ZSRC column in the NSA table). This resulted in a final sample of 81,462 dwarf galaxies.
For the ZTF photometry, we also required that the objects overlap with > 100 high quality g and r band ZTF images over a > 1 year baseline. Due to computational limitations, we selected a random subset of the parent light curve sample consisting of 25,714 dwarf galaxies for our ZTF variability search. For our WISE variability search we required there to be > 5 epochs over a > 3 year baseline and this resulted in a final sample of 79,879 objects for the mid-IR search.
We produced a control sample of optically variable AGN with host galaxies of mass M * > 10 9.75 M from  Figure 1. Pearson r correlation coefficient and χ 2 /N in g-and r-band calculated from ZTF light curves for 3 different binning timescales. We show the entire dwarf galaxy population with ZTF photometry in blue and the AGN control sample (from host galaxies of stellar mass M * > 3 × 10 9.75 M ) in green. The cutoffs used for AGN candidate selection are shown in black dotted lines for each statistic. We required the 3 statistics to satisfy the cutoffs for all 3 binning timescales in order for a candidate to be selected. a parent sample of 5,493 variable ZTF AGN found by Ward et al. (2021). This AGN sample was obtained from the ZTF alert stream (Patterson et al. 2019) using the AMPEL alert broker AMPEL (Alert Management, Photometry and Evaluation of Lightcurves; Nordin et al. 2019) to crossmatch ZTF transients to AGN catalogs using catsHTM (Soumagnac & Ofek 2018) and check for AGN-like variability with the Butler & Bloom (2011) quasar modeling routine. We matched the AGN from this sample to the NASA-Sloan Atlas with a 5" radius and found 1,053 objects with measured galaxy masses. Both the control sample and the dwarf galaxy sample were processed by the following photometry pipeline.

ZTF PHOTOMETRY OF DWARF GALAXIES
Difference imaging to detect variability requires the production of reference images, which are co-added stacks of exposure images to support image differencing downstream. To produce new ZTF forced photometry with custom, deeper references than the main ZTF alert pipeline we implemented the ZUDS photometry pipeline 2 on our dwarf galaxy and AGN control samples. We generated reference images for each field, CCD, quadrant and filter from ZTF single epoch science images by selecting up to 60 images as close as possible in time which met the following criteria established by Masci et al. (2019) for the main ZTF alert pipeline: 1. Seeing within the range 2."0 ≤ FWHM ≤ 5."0, with priority given to images with lower seeing values.

Color coefficients given by
.15 or −0.05 ≤ CLRCOEFF(r) ≤ 0.22 for filters g and r respectively.
6. Global pixel median ≤ 1900 DN or ≤ 1600 for filters g and r respectively.
When there were multiple field, CCD and quadrants containing the source for a particular filter, we selected the reference image containing the largest number of high quality ZTF science images in the coadd to be the reference image. We produced 1000"x1000" cutouts of each single epoch science image and reference image, produced subtractions, then undertook aperture photometry with a 3".0 radius. We then applied the aperture to PSF correction factors produced by the main ZTF image pipeline to produce PSF photometry light curves. We measured the baseline flux of each object in the reference image and added this to the fluxes measured from the image subtractions.
In order to improve our sensitivity to low S/N variability from faintly varying AGN and prepare the data for calculation of the Pearson correlation coefficient between g-and r-band photometry, we binned the data in temporal bins. We first applied zeropoint corrections, then undertook error-weighted binning in flux space using bins of 5, 10 and 20 day increments. Therefore, each galaxy had 3 light curves with different bin sizes for the calculation of variability statistics. This binning procedure may have reduced our sensitivity to objects with optical variability only on timescales <20 days and therefore may have introduced biases against particularly low mass AGN. However, we determined that binning was necessary to both allow for the calculation of the Pearson r correlation coefficient whilst confidently identifying light curves with correlated variability over a range of timescales, as expected from AGN power spectra (Kelly et al. 2009;MacLeod et al. 2011;Burke et al. 2021) With this procedure we generated light curves of 25,714 dwarf galaxies and the 1,053 AGN from the high mass galaxy control sample.

SELECTION OF VARIABLE IMBH CANDIDATES
In order to detect real variability amongst the sample of dwarf galaxy light curves we trialed a number of variability statistics including the Pearson correlation coefficient (as used by Secrest & Satyapal (2020)), the goodness of fit of the light curves to the quasar structure function (using qsofit (Butler & Bloom 2011), as used in Baldassare et al. (2018Baldassare et al. ( , 2020b; Ward et al. (2021)), the excess variance (Sánchez et al. 2017;Martínez-Palomera et al. 2020) and the χ 2 /N in g and r bands. We found that a combination of the Pearson correlation coefficient and χ 2 produced the best separation between the AGN control sample and the majority of the non-variable dwarf galaxy population. The Pearson correlation coefficient r between the binned g and r band fluxes was calculated as: where C fg,fr is the covariance between g and r bands given by and σ fg and σ fr are the g and r band variability amplitudes given by: and the expectation value f was given by the median flux.
We also calculated the χ 2 /N of the light curves in g and r bands: The distribution of the Pearson correlation coefficient r and the χ 2 /N values for g and r bands for the dwarf galaxy and control AGN light curves with 5, 10 and 20 day bin sizes is shown in Figure 1. We applied cutoffs of r > 0.2, r > 0.3, and r > 0.4 for 5, 10 and 20 day bin sizes respectively, and χ 2 /N > 3 in both filters for all bin sizes, to classify dwarf galaxies as variable. Each dwarf galaxy was required to meet this cutoff in all 3 binning timescales to ensure that correlations and variance was not produced as artifacts of bin phase and size. After applying these cutoffs, we found 130 dwarf galaxies with statistically significant variability. We then manually inspected the light curves and difference images to remove light curves with high variance due to poor quality photometry and to determine which light curves contained single flares with reddening consistent with supernovae. 36 supernovae-like flares were found, corresponding to a nuclear SN rate of 0.14±0.02% (0.05 ± 0.01% year −1 ) for dwarf galaxies in ZTF Phase I.
This process produced a final sample of 44 dwarf galaxies (0.17 ± 0.03%) with variability consistent with AGN activity. These constitute our set of optically variable IMBH candidates. The properties of this sample are summarized in Table 1. Examples of ZTF light curves of 3 dwarf galaxies with AGN-like variability are shown in Figure 2 and examples of 3 supernova light curves are shown in Figure 3. The distributions of the redshifts and host galaxy stellar masses of the optically variable dwarf galaxies are shown in Figure 4. All but 6 of the IMBH candidates were in galaxies of mass M * > 10 9 M . By comparison, a larger fraction of supernovae were observed in low mass galaxies and most were found in a narrow redshift range of 0.02 < z < 0.055.
We checked both the supernova and AGN candidate samples for spectroscopic classification on the Transient Name Server 3 . 25 of the 36 supernova candidates had published spectroscopic classifications. One of these was classified as a supernova but was not typed. 17 out of Note-Properties of the 36 supernova candidates with significant and correlated g and r band variability found in ZTF. The IDs, positions, redshifts and host galaxy stellar masses are those from the NSA catalog version v1 0 1. The presence of Balmer broad lines is indicated in the BLR column. The last column shows the spectroscopic classification made published by ZTF on the Transient Name Server.
the remaining 24 (71%) were classified as SNIa, SNIc or SNIc-BL. 7 out of 24 (29%) were classified as SNII or SNIIn. These classifications are shown in Table 2. We note that the remaining 11 SN without spectroscopic classifications can only be classified as SN candidates, as they may be AGN outbursts (Drake et al. 2019). None of the AGN candidates had spectroscopic classifications published on TNS.
In order to estimate the black hole masses of our IMBH candidates based on their host galaxy stellar mass we used the updated relationship between black hole mass and bulge mass derived by Schutte et al. (2019). This study found a linear relationship based on a sample of galaxies with carefully measured black hole and bulge masses including 8 dwarf galaxies with stellar masses log(M * ) < 8.5M : The relation has a scatter of 0.68 dex. The black hole mass estimates for the IMBH candidates are shown in Table 1.

SPECTROSCOPIC AND MULTI-WAVELENGTH PROPERTIES OF THE ZTF-SELECTED IMBH CANDIDATES
In order to determine the spectroscopic class of the IMBH candidates on the Baldwin, Phillips and Terlevich NSA164967 NSA136891 NSA33420 Figure 3. Three examples of ZTF light curves which passed the variability criteria and were classified as supernova candidates.  diagram (Baldwin et al. 1981;Veilleux & Osterbrock 1987) we fit the narrow emission line ratios with Penalized Pixel Fitting (pPXF) (Cappellari & Emsellem 2004;Cappellari 2016). This method finds the best fit stellar continuum and absorption model using a large sample of high resolution templates of single stellar populations adjusted to match the spectral resolution of the input spectrum. We simultaneously fit the narrow Hα, Hβ Figure 5. 35 objects (81%) are classified as star forming while 4 objects (9%) are in the composite region and 4 (9%) are classified as Seyferts. For comparison, we also show the density of the original 81,462 dwarf galaxy parent sample on the BPT diagram using emission line ratios quoted in the NSA catalog. From the grey contours it can be seen that the majority of the dwarf galaxy population is star forming but a small population extends out to AGN and LINER regions of the BPT classification scheme. Our typical star forming IMBH candidates lie at higher emission line ratios than this population, in part because our pPXF modeling considers stellar absorption of those lines.
The classification of the supernova host galaxies on the [O iii] λ5007/Hβ -[N ii] λ6583 /Hα BPT diagram is also shown in Figure 5. 34 host galaxies (94%) are classified as star forming while 1 object (3%) was in the composite region and 1 (3%) in the Seyfert region.
We crossmatched our ZTF-selected IMBH candidates to the active dwarf galaxies from the Secrest & Satyapal (2020) mid-IR variability search, the Baldassare et al. (2018Baldassare et al. ( , 2020b  NSA32653, NSA35747, NSA49405, NSA164884, NSA464884, NSA181600 and NSA451469 were the only objects in our IMBH sample to exhibit broad Balmer lines in archival SDSS spectra. These objects previously had their virial black hole masses estimated using the width of the Hα broad lines by Ho & Kim (2015) and Liu et al. (2019). These virial masses ranged between 10 6.3 M and 10 7.6 M and are shown in Table 1.

ZTF VARIABILITY OF PREVIOUSLY REPORTED PTF-SELECTED IMBH CANDIDATES
After discovering that 5 objects from the variabilityselected AGN sample from PTF (Baldassare et al. 2020b) were also variable in ZTF according to our selection thresholds, we decided to determine if the remaining AGN candidates in dwarf galaxies from PTF had ZTF variability which was missed by our selection criteria.
We first found that our parent dwarf galaxy sample had 152 objects which overlapped with the variable PTF sample. We then visually inspected the zuds-pipeline light curves we made of the 152 common dwarf galaxies and confirmed that no other dwarf galaxy had apparent variability. We then used the ZTF forced photometry service (Masci et al. 2019) to obtain alternative photometry of the sources with the original ZTF reference images. After removing poor quality images by requiring the procstatus flag be = 0, we measured the baseline flux from the reference images, applied zeropoints and combined the baseline and single epoch fluxes to produce the light curves of the 152 objects. We then visually inspected the candidates to look for any signs of variability. The alternative pipeline confirmed that only 5 candidates from our overlapping samples showed statistically significant variability in ZTF.

WISE SINGLE EPOCH FORCED PHOTOMETRY
As we aimed to search for mid-IR variability from a large sample of dwarf galaxies which may have been too faint to appear in the search of the original AllWISE catalog by Secrest & Satyapal (2020), we made use of forward modelled photometry of time-resolved WISE coadds (Meisner et al. 2018) made available through Data Release 9 4 of the DESI imaging Legacy Survey (Dey et al. 2019). This approach was previously implemented by Lang et al. (2016b) to produce forced photometry of 400 million SDSS sources with the deep unWISE coadds (Lang 2014;Meisner et al. 2017). Lang et al. (2016b) used The Tractor (Lang et al. 2016a) to use deeper, higher resolution source models from SDSS to produce photometry of blended and faint objects in the unWISE coadds. They were able to report fluxes and uncertainties from 3σ and 4σ detections which were not included in the original WISE catalog. More recently, they implemented this technique to produce time-resolved WISE photometry (Meisner et al. 2018). As WISE revisits each field at a ∼6 month ca-dence, with an increased cadence towards the poles, and takes ∼ 10 exposures each visit, they coadded the exposures from each visit to produce forced photometry on each coadd. This therefore provided light curves with approximately 15 fluxes measured over ∼ 8 years.
We crossmatched the SDSS galaxy positions of our dwarf galaxy sample to the unWISE source catalog and pulled the single epoch WISE photometry for the closest unWISE source. Some light curves showed an oscillation behaviour when multiple sources were contained within the large WISE PSF and the WISE flux was distributed across the multiple sources differently in each epoch. To overcome this, we combined the total flux of the sources within a radius of 3 WISE pixels (3x2.75") to ensure that all dwarf galaxy variability was captured within the combined fluxes. We removed light curves where an erroneously high flux ( f − f > 5σ) in both W1 and W2 bands produced a r and χ 2 /N above the cutoffs.
We also removed light curves where source confusion within the WISE PSF size still resulted in an oscillating behaviour using the following criteria. If the summed difference between every adjacent W1 flux measurement was greater than twice the summed difference between every W1 flux offset by two epochs, the light curve was flagged as bad quality. We produced WISE light curves of the AGN control sample with the same procedure.
For each source we calculated the Pearson correlation coefficient r between the W1 and W2 bands and the χ 2 /N for each band. The distribution of these statistics for the two samples is shown in Figure 6. We required that dwarf galaxies have r > 0.75 and χ 2 /N > 1.0 in W1 and W2 to be considered a variable AGN candidate.
Of the 79,879 dwarf galaxies for which WISE single epoch forced photometry was available, 124 were removed due to light curve quality flags. Of the remaining 79,755 light curves, we found that 165 had fluxes which met our variability criteria. One dwarf galaxy, NSA253466, was detected by our variability criteria due to the well-studied Type IIn supernova SN 2010jl (Stoll et al. 2011) which has had detailed follow-up in the near-IR, X-ray and radio (e.g. Fransson et al. 2014;Chandra et al. 2015) showing interaction with the dense circumstellar medium. We removed this object from our AGN candidate sample. We removed 3 other objects with supernovae visible in both WISE and ZTF data: NSA559938 (ZTF18aamftst: an SNIIn), NSA143427 (ZTF18acwyvet, an SNIIL), and NSA230430 (ZTF20aaupkac: an SNIa). We removed 14 other galaxies with light curves showing single flares, often with color changes characteristic of SNe: NSA20892, NSA32356, NSA143207, NSA236644, NSA250558, NSA253466, NSA274965, NSA340533, We show the entire dwarf galaxy population with WISE photometry in blue and the AGN control sample (from host galaxies of stellar mass M * > 3 × 10 9.75 M ) in green. The cutoffs used for AGN candidate selection are shown in black dotted lines for each statistic. We required the 3 statistics to satisfy the cutoffs in order for a candidate to be selected.
NSA355173, NSA379733, NSA475418, NSA502699, NSA528212, NSA548379. The properties of these supernova host galaxies are summarized in Table 3. The properties of the remaining final sample of 148 AGN candidates (corresponding to 0.19 ± 0.02% of the dwarf galaxy sample) are summarized in Table 4 and 6 examples of variable WISE light curves are shown in Figure 7.

PROPERTIES OF THE WISE -SELECTED VARIABLE IMBH CANDIDATES
The distributions of the redshifts and host stellar masses of the mid-IR variable dwarf galaxies are shown in Figure 8. We exclude NSA64525 from the histogram because the estimated mass provided by the NASA-Sloan Atlas (10 5.39 M ) is inconsistent with the stellar dispersion velocity of 193 ± 8 km/s measured by the SDSS spectroscopic pipeline. This dispersion velocity indicates that the BH itself likely has a mass of ∼ 10 7 M according to the M −σ * relation (Kormendy & Ho 2013) and the host galaxy therefore has mass comparable to ∼ 10 9.5 M based on the M BH − M * relation (Schutte et al. 2019). The mid-IR variability selection finds a higher fraction of variable AGN at low redshifts compared to the overall dwarf galaxy sample and shows a slight preference for higher mass galaxies.
The classification of the mid-IR IMBH candidates on the [O iii] λ5007/Hβ -[N ii] λ6583 /Hα BPT diagram is shown in Figure 9. 100 objects (69%) are classified as star forming while 32 objects (22%) are in the composite region, 3 (2.1%) were classified as LINERs and 10 (6.9%) as Seyferts. 12 of the IMBH candidates (8.1%) have broad Balmer lines. These objects also had their virial black hole masses estimated using the width of the Hα broad lines by Ho & Kim (2015) and Liu et al. (2019). These virial masses range between 10 5.5 M and 10 8.2 M respectively and are shown in Table 4.
In order to determine how many WISE -selected IMBH candidates were also optically variable, we produced ZTF photometry of the remainder of the 148 dwarf host galaxies which were not included in the original ZTF search. We found that 15 out of 148 (10%) met the optical variability criteria for ZTF outlined in Section 3.5. 7 of the 15 AGN which were variable at both wavelengths had visible broad lines in their spectra, and amongst these broad line AGN, the virial masses ranged between M BH = 10 6.3 M and M BH = 10 8.2 M . Only two of these sources (NSA164884 and NSA451469) were found in the original ZTF sample, due to the smaller sample of dwarf galaxies for which we produced optical photometry. The AGN with both mid-IR and optical variability are indicated in the last column of Table 4. NSA451469 had also been found in PTF by Baldassare et al. (2020b). We crossmatched our mid-IR selected IMBH candidates to the active dwarf galaxies from the Secrest & Satyapal (2020) mid-IR variability search, the Baldassare et al. (2018Baldassare et al. ( , 2020b   One object, NSA638093, had previously been found in the WISE variability search by Secrest & Satyapal (2020) and in the mid-IR color selection box by Latimer et al. (2021) (object number 11, listed with NSAID 151888 from version v 1 1 2 of the NSA), where they used Chandra to find X-ray emission which may have been consistent with X-ray binaries instead of AGN activity. NSA386591 appears in the Reines et al. (2013) spectroscopic search for BPT AGN and Composites and the Latimer et al. (2021) mid-IR color selection search (ID number 6) where it was found to have an X-ray luminosity consistent with AGN activity and too large to be produced by X-ray binaries.
We undertook a search for radio emission from the IMBH candidates in the Karl G. Jansky Very Large Array Sky Survey (VLASS; Lacy et al. 2020). This survey covers a total of 33,885 deg 2 in the 2-4 GHz range with an angular resolution of ∼ 2".5 and will obtain a coadd 1σ sensitivity of 1 µJy/beam by survey end in 2024. We searched for crossmatches within 10" in Table 2 of the VLASS Epoch 1 Quick Look Catalogue which contains ≈ 700, 000 compact radio sources with > 1 mJy/beam detections associated with mid-IR hosts from the unWISE catalog (Gordon et al. 2020). We found that one IMBH candidate, NSA87109, had a corresponding radio point source at a separation of 0.0038" from the optical NASA-Sloan Atlas galaxy position with a flux of 17.55 ± 0.27 Jy. This object is a known BL Lac which has also been detected in 0.1-2.4 keV X-rays (Massaro et al. 2009).
We show cutouts from the DESI Legacy Imaging Surveys (Dey et al. 2019) of 3 WISE candidates and 2 ZTF candidates with host masses M * < 10 8.2 M and the object with a VLASS radio detection (NSA87109, with stellar mass log 10 M * = 9.66) in Figure 10. The radio source and one other low mass IMBH candidate are in compact, blue galaxies while the other 4 low mass IMBH candidates reside in galaxies with complex morphologies and multiple stellar overdensities.

DISCUSSION
The number of dwarf galaxies which were variable in WISE corresponded to a 0.19 ± 0.02% variability fraction, while the optical variability fraction that we find from the ZTF AGN candidates is 0.17 ± 0.03%. Our results therefore suggest that the two methods are similarly effective for identifying IMBH candidates with the current sensitivities and baselines available for mid-IR and optical photometry of dwarf galaxies. This, however, will change with the improved optical sensitivities and baselines offered by the LSST survey over the next decade.
Our mid-IR active fraction was larger but within the uncertainty range of the active fraction found by Secrest & Satyapal (2020) (0.09 +0.20 −0.07 %) in a much smaller sample of 2197 dwarf galaxies with the main WISE photometry catalog. The optical variability fraction that we find is also consistent with the active fraction of 0.15±0.07% found for dwarf galaxies of the same mass in PTF (Baldassare et al. 2020b). The 90% fraction of WISE -selected AGN candidates which are variable in the mid-IR but not the optical likely arises from a combination of line-of-sight obscuration of nuclear optical emission due to the dusty torus and global obscuration from the host galaxy. Obscuration of nuclear optical emission from AGN with detectable mid-IR signatures has been observed for many Seyfert 2 galaxies (e.g. Goulding et al. 2011;Annuar et al. 2015;Ricci et al. 2016;Annuar et al. 2017). It is possible that the majority of our mid-IR variable AGN in dwarf galaxies are Seyfert 2s with obscured optical variability which are also not picked up by the BPT diagram due to the effects of their lower masses on the optical emission line ratios. Line-of-sight obscuration of optical emission of WISE candidates is supported by the fact that 7 of the 12 WISE -selected IMBH candidates with bright BLRs were variable in both the mid-IR and the optical.
We note that only 5 of the 152 dwarf galaxies in common with the Baldassare et al. (2020b) AGN candidate sample from PTF were variable in ZTF. This may be because of the use of differing statistical criteria for variability classification. It may also be the case that the longer 7 year baseline of the combined PTF and iPTF light curves in Baldassare et al. (2020b) improved their sensitivity to AGN which vary on longer timescales compared to shorter timescales. They indeed find that longer baseline data has a much higher detectable variability fraction, increasing by a factor of 4 from 0.25% for light curves with < 2 year baselines to 1% for light curves with > 2 year baselines. By comparison, our procedure of using deep references and stacking to detect fainter variability in more evenly-sampled, ∼ 3 year baseline light curves may make us more sensitive to variability on month to year long timescales, but may miss AGN with flux changes over longer 2-7 year timescales. An alternative explanation may be that a large fraction of low mass AGN are state-changing AGN (e.g. Frederick et al. 2019) which can switch their optical variability on and off over the decade-long timescale of the PTF and ZTF surveys. Indeed, Martínez-Palomera et al. (2020) found that the majority of IMBH candidates with optical variability on hourly to daily timescales from the SIBLING survey were not variable when observed again the following year.
81% of the ZTF-selected IMBH candidates were starforming on the BPT diagram and only 7 have been identified as AGN via their Balmer broad lines. Only 7 had been identified in previous dwarf galaxy AGN searches and these were all via their optical variability in SDSS or PTF. The non-AGN spectroscopic classifications of the majority of the sources indicates that optical variability selection can find AGN in dwarf galaxies which would be missed by other selection strategies.
Similarly, 100 (69%) of WISE -selected IMBH candidates were star-forming on the BPT diagram and therefore would have been missed by classic spectroscopic selection methods. 12 (8.1%) can be identified as AGN due to the presence of broad lines, 1 (0.63%) could have been found via radio emission alone, another 2 (1.27%) via the mid-IR selection box and 1 (0.63%) via X-ray emission. We therefore see that ∼ 70% of candidates from mid-IR variability selection could likely not have been found through other selection techniques. A higher fraction of the WISE -selected candidates are BPT-AGN compared to optically-selected candidates, perhaps indicating that the AGN emission lines of galaxies with mid-IR variability are less likely to be diluted by star formation.
Both the ZTF-selected and WISE -selected AGN candidate host galaxies tend to have higher masses and lower redshifts compared to the overall distributions of the host galaxy sample, likely due to the higher luminosities of these AGN. However, our selection method is still capable of detecting AGN variability from dwarf galaxies with redshifts up to z ∼ 0.15 and stellar masses down to M * = 10 7.52 M . The virial masses for the AGN candidates with broad lines go down to M BH = 10 5.485 M and in most cases are lower than the BH masses estimated from the M BH − M * relation, indicating that our search may have found MBHs which are undermassive for their hosts. We therefore conclude that our variability-selection approach is useful for selecting AGN which can populate the poorly sampled lower end of the M BH − σ * relation. Future work should take high resolution spectra of our IMBH candidates for fitting of the velocity dispersion from the stellar absorption lines to provide another independent estimate of the BH mass.
The discovery of 36 nuclear supernova candidates shows the usefulness of applying simple statistics to ZTF forced photometry of large galaxy samples to find supernovae candidates in dwarf galaxies. There are many motivations to study rates of supernovae in dwarfs such as explaining the increase in the rate ratio of superluminous supernovae to core collapse supernovae in low mass galaxies and whether the increased metallicity or specific star formation rates in dwarfs are the driving factor for this trend (Taggart & Perley 2021). It has also been found that the emerging class of fast blue optical transients like AT2018cow are preferentially hosted by dwarf galaxies . Our finding that 0.14% of dwarf galaxies contained nuclear supernovae during ZTF Phase-I (corresponding to 0.05% per year) and that most SN were within redshifts of 0.02 < z < 0.055 may provide insights into these questions.

CONCLUSIONS
In this paper we have presented a search for IMBH candidates by looking for variable AGN in dwarf galaxies of stellar mass M * < 10 9.75 M in the optical and mid-IR. We applied a new ZTF forced photometry pipeline to produce deep, high quality reference images for image subtraction and made g-and r-band light curves of 25,714 dwarf galaxies. These light curves were stacked in a range of time bins to improve sensitivity to faint variability. We applied statistical cutoffs to find significant and correlated variability between the two bands and found 36 supernova candidates and 44 AGN candidates. The supernova fraction was 0.05 ± 0.01% year −1 and the optically variable AGN fraction was 0.17 ± 0.03%.
To search for mid-IR variability we used Tractor forward modelled photometry of time-resolved WISE coadds. We found 148 dwarf galaxies with significant and correlated variability in the W1 and W2 bands after removing 14 supernovae. The mid-IR variable AGN fraction was 0.19 ± 0.02%.
We found that 81% of our ZTF AGN candidates would have been missed with classical spectroscopic classification on the BPT diagram. Of our WISE candidates, 69% would have been missed with spectroscopic classification and only 4 would have been detected via radio or X-ray detection or the mid-IR color selection box. While our candidates were slightly biased to low redshifts and high galaxy masses compared to the parent dwarf galaxy sample, they were effective in identifying AGN with virial masses as low as M BH = 10 5.485 M and in dwarf galaxies with stellar masses 10 7.5 < M * < 10 8.5 M . We therefore conclude, in accordance with previous variability searches, that optical and mid-IR variability selection is effective for finding low mass AGN in dwarf galaxies which would be missed by other spectroscopic selection techniques.
After checking the ZTF photometry of the 152 dwarf galaxies in common between our parent sample and the variability-selected AGN candidates from PTF (Baldassare et al. 2020b) we found that only 5 continue to show variability in ZTF. The lack of variability in our ZTF light curves may be due to the different baselines and sensitivities of the two search strategies and the differing methods for finding AGN-like variability.
With more detailed imaging and spectroscopic analysis, our variability-selected AGN candidates could help to populate the sparsely sampled end of the M * −σ relation and provide insights into black hole seed formation mechanisms, dwarf galaxy-black hole co-evolution and the accretion states of low mass AGN. Future work on these candidates will include forward modeling of ZTF and DECam images to determine their positions relative to their host galaxy nuclei and determine the off-nuclear fraction. The potential for forward modeling of ZTF images to determine the position of a variable point source relative to its host galaxy was demonstrated in Ward et al. (2021) for recoiling SMBH candidates and will provide a way to confirm the positions of IMBH candidates without X-ray or radio detections.
These candidates are just the tip of the iceberg in the search for optically variable AGN in dwarf galaxies, which will be greatly enhanced by the capabilities of the Legacy Survey of Space and Time (LSST) at Vera C. Rubin Observatory (Ivezić et al. 2019) over the next decade. The capacity of LSST to find fainter and more distant IMBHs in dwarf galaxies via their variability will tighten the constraints we can place on black hole seeding channels and the efficiency of massive BH growth. servatory, the University of Nottingham, the Ohio State University, the University of Pennsylvania, the University of Portsmouth, SLAC National Accelerator Laboratory, Stanford University, the University of Sussex, and Texas A&M University.