Neutrino follow-up with the Zwicky Transient Facility: Results from the first 24 campaigns

The Zwicky Transient Facility (ZTF) performs a systematic neutrino follow-up program, searching for optical counterparts to high-energy neutrinos with dedicated Target-of-Opportunity (ToO) observations. Since first light in March 2018, ZTF has taken prompt observations for 24 high-quality neutrino alerts from the IceCube Neutrino Observatory, with a median latency of 12.2 hours from initial neutrino detection. From two of these campaigns, we have already reported tidal disruption event (TDE) AT 2019dsg and likely TDE AT 2019fdr as probable counterparts, suggesting that TDEs contribute >7.8% of the astrophysical neutrino flux. We here present the full results of our program through to December 2021. No additional candidate neutrino sources were identified by our program, allowing us to place the first constraints on the underlying optical luminosity function of astrophysical neutrino sources. Transients with optical absolutes magnitudes brighter that − 21 can contribute no more than 87% of the total, while transients brighter than − 22 can contribute no more than 58% of the total, neglecting the effect of extinction and assuming they follow the star formation rate. These are the first observational constraints on the neutrino emission of bright populations such as superluminous supernovae. None of the neutrinos were coincident with bright optical AGN flares comparable to that observed for TXS 0506+056/IC170922A, with such optical blazar flares producing no more than 26% of the total neutrino flux. We highlight the outlook for electromagnetic neutrino follow-up programs, including the expected potential for the Rubin Observatory.


INTRODUCTION
Astrophysical neutrinos are produced through the interaction of accelerated hadrons with matter or photons.A flux of astrophysical neutrinos with energies in the TeV-PeV range, was first discovered by IceCube in 2013 (IceCube Collaboration 2013).Recent results suggest that a substantial fraction of these high-energy neutrinos are produced in the cores of Active Galactic Nuclei (AGN) (Abbasi et al. 2022), with additional evidence for neutrino emission from the nearby AGN NGC 1068 (Aartsen et al. 2020).Beyond this static component, various transient or variable source classes have been proposed as possible contributors to the neutrino flux, including gamma-ray bursts (GRBs) (Waxman & Bahcall 1997), core-collapse supernovae (CCSNe) (Murase et al. 2011), TDEs (Farrar & Gruzinov 2009) and blazars (Mannheim 1993).All of these proposed neutrino source classes have electromagnetic signatures at optical wavelengths.
The Zwicky Transient Facility (ZTF) is an optical telescope with a 47 sq.deg field of view (Bellm et al. 2019a;Dekany et al. 2020).Since ★ E-mail: rdstein@caltech.edu† LSSTC Data Science Fellow first light in 2018, ZTF has operated a dedicated neutrino follow-up program, in which the arrival directions of IceCube neutrino alerts are observed with Target-of-Opportunity (ToO) observations (Graham et al. 2019).This program has led to the identification of two further likely high-energy neutrino sources, the TDE AT 2019dsg (Stein et al. 2021b) and the probable TDE AT 2019fdr (Reusch et al. 2022, though see Pitik et al. 2022 for an alternative interpretation).Accounting for the contribution of higher-redshift sources, these results suggest that at least 7.8% of neutrino alerts arise from the broader TDE population (Reusch et al. 2022).Archival analysis of ZTF data revealed further evidence of a correlation between such flares and high-energy neutrinos (van Velzen et al. 2021).
In this paper we outline the full results of the ZTF neutrino followup program, which has to date included 24 dedicated neutrino followup campaigns.This sample enables novel constraints to be set on the neutrino emission of a broad range of optical transient and variable populations.
The paper is organised as follows: Section 2 outlines the program itself, including trigger criteria and optical candidate selection.Section 3 outlines transient candidates identified by the program, and subsequent electromagnetic observations to determine their nature.Section 4 outlines optical AGN flares found coincident with neutrinos, and Section 5 provides data on two candidate neutrino sources identified in the literature.Section 6 considers the various constraints that can be placed on different possible neutrino source populations from our program.Section 7 summarises the main results, and outlines how such follow-up programs may improve with future observatories.

NEUTRINO FOLLOW-UP WITH ZTF
Neutrino alerts are generally published by IceCube in the form of automated Gamma-ray Coordination Network (GCN) Notices1 , with initial estimates of the statistical uncertainty on the neutrino position.These positions are then superseded after a few hours by a GCN Circular with an updated localisation that also incorporates systematic uncertainties (Lagunas Gualda et al. 2021).Given the substantial increase in localisation area once systematic effects are accounted for, with increases of factor 5 not being uncommon, we rely on the latter category to perform our search for neutrino counterparts.
With ZTF, we aim to observe all accessible high-quality neutrino alerts from IceCube.We define high-quality alerts as those with a high probability to be of astrophysical origin ('signalness' > 50%), or those which are well-localised (a 90% localisation area < 10 sq.deg.).Though IceCube labels alerts as Gold or Bronze based on average quality, individual Bronze alerts have been reported with signalness values greater than 50% (e.g.IC211208A) and Gold alerts have been reported with signalness values less than 15% (e.g.IC201130A).We therefore ignore the labelling of these streams, and select exclusively based on the signalness and localisation.
We have followed up 24 neutrinos in the period from survey start on 2018 March 20 to 2021 December 31, out of a total of 79 neutrino alerts published by IceCube during that time.Table 1 summarises each neutrino alert observed by ZTF.From 2019 June 17, IceCube published neutrino alerts with improved selection criteria (V2) to provide an elevated alert rate (Blaufuss et al. 2019).In addition to 1 of the 12 alerts under the old selection, ZTF followed up 23 of the 67 alerts published under the V2 selection.Midway through the ZTF program, an additional cut on neutrino alert galactic latitude (|b| > 10 deg) was introduced to avoid crowded fields with many stars.
Each neutrino localisation region can typically be covered by one or two observations of fields in a predefined ZTF 'grid' tiling of the sky.Multiple observations are scheduled for each field, with both  and  filters, and a separation of at least 15 minutes between images.These observations typically last for 300 s, with a typical limiting magnitude of 21.5.ToO observations are typically conducted on the first two nights following a neutrino alert, before swapping to serendipitous coverage with shorter 30 s exposures and a 2-day cadence as part of the public survey (Bellm et al. 2019b).As can be seen in Figure 1, our first coverage of events has a median latency of 12.2 hours from neutrino detection.Some latency is unavoidable because the neutrino localisation itself is typically only released with a delay of ≳2 hours, but additional latency arises primarily due to observability constraints.Poor weather can prevent observations on the first night after neutrino detection, leading to 20% of alerts observed with a latency >24 hours.Serendipitous coverage from the public survey, with a median latency of 24 hours after neutrino detection, reduces the latency for some campaigns.
As for all ZTF data, these observations are first processed by the Infra-red Processing and Analysis Centre (IPAC) to identify detec-tions in difference images (Masci et al. 2019).These detections are then packaged as 'alerts' (Patterson et al. 2019), and processed by our dedicated data analysis pipeline, NuZTF (Stein et al. 2021a), which searches for extragalactic ZTF detections coincident with external triggers.For neutrinos followed-up by ZTF, we define spatial coincidence as requiring that an object lies within the reported 90% localisation rectangle from IceCube, and define temporal coincidence as requiring that an object is detected at least once following the neutrino arrival time.
NuZTF is built using the AMPEL software framework (Nordin et al. 2019), based on a search algorithm for extragalactic transients.Cuts are applied to reject spurious detections, stars and solar system objects (see Stein et al. 2021b for more details).Searching for detections in the window from neutrino arrival time to 14 days postneutrino, these cuts typically yield 1 good candidate per ∼3 sq.deg. of observed sky.
Promising candidates are prioritised for spectroscopic classification, to confirm or rule out a possible association with a given neutrino.Once classified, an object can then be cross-referenced to relevant neutrino emission scenarios for that population.In particular, optical signatures we look for include: • Supernovae with evidence of CSM interaction.High-energy neutrinos are thought to be produced when CCSNe occur within a dense circumstellar medium (CSM), with the resultant shock collisions then generating neutrino emission simultaneously with the optical lightcurve (Murase et al. 2011).The presence of such CSM interaction also results in characteristic narrow lines in the optical spectrum, so these models generally apply to the Type IIn supernova population which exhibits these lines.The neutrino emission is expected to be highest close to optical peak, and to then decay over time.In this case, the expected optical signature would be any supernova with evidence of ongoing CSM interaction.
• Supernovae with relativistic jets.Some supernovae have been observed to launch relativistic jets as part of the core-collapse process (Galama et al. 1998).Those jets which proceed to escape the surrounding stellar envelope and CSM can be observed as long GRBs if they are oriented towards Earth.Analogously, where an on-axis supernova jet does not escape the stellar envelope, there would instead be a so-called 'choked jet' (Nakar 2015).For both scenarios, neutrino emission would primarily be expected during the 'prompt phase', in the ∼100s after supernova explosion (Waxman & Bahcall 1997;Senno et al. 2016).This scenario would then lead to a young supernova, typically of Type Ic-BL, appearing at the location of the neutrino.The supernova would have an explosion time compatible with the neutrino detection time, and since SNe brighten over a period of days, this optical signature would be delayed relative to the neutrino itself.
• GRB Afterglows.Another signature of the supernova jet scenario would be the direct detection of a long-GRB afterglow.Models have also predicted neutrino emission for short GRBs, so a short-GRB afterglow could also be a potential counterpart (Waxman & Bahcall 1997).These GRB afterglows would not be detected before the neutrino detection, and would fade rapidly over the next few hours before falling below the ZTF detection threshold.
• AGN Flares.AGN flares, and especially blazar flares, have been suggested as neutrino sources (Bednarek & Protheroe 1999), though the neutrino emission itself would not necessarily be directly correlated to the optical emission.For example, for the standard two-hump Spectral Energy Distribution (SED) model, the optical emission could serve primarily as a tracer for photon target density but not necessarily PeV proton luminosity.We restrict ourselves to searches for AGN undergoing significant optical flaring coincident with a neutrino.Neutrinos could also be produced in AGN without coincident optical flares, but such neutrino emission scenarios are not best probed with an optical follow-up program such as ours.
• Tidal Disruption Events.TDEs have been suggested as neutrino sources, through multiple emission channels such as jets, outflows or in coronae (see Hayasaki 2021 for a recent review).The timescale for neutrino production remains unclear, but would not be expected prior to the TDE itself.Non-thermal emission from TDEs can last several hundred days, so the signature in this case would be any 'ongoing' TDE coincident with a neutrino.
We do not enforce any additional cut on candidate distance, because IceCube detects neutrinos emitted throughout the universe rather than being restricted only to the local universe (see Strotjohann et al. 2019 for a more detailed explanation of this effect).However, given the limiting magnitude of ZTF, the candidates we find will nonetheless be biased towards lower redshifts.
We do not explicitly reject objects with a history of variability, because variable objects have been proposed as possible neutrino sources.However, our program is intended to identify increased optical flux that is contemporaneous with a neutrino's detection, so only variable objects with significantly enhanced flux relative to reference images are selected by our pipeline.The blazar flare of TXS 0506+056 fell into this category (IceCube Collaboration et al. 2018), and we would be capable of identifying similar examples.
To date, the NuZTF pipeline has identified 172 candidates for visual inspection out of an observed area of 154.33 sq.deg across 24 neutrinos, using a search window of 14 days after each neutrino detection.This corresponds to an initial density of 1.05 candidates per sq.deg. of sky.The full list of candidates for each neutrino is given in the Appendix.
Visual inspection then enables us to further classify objects and reject background detections.Viewing difference images directly enables us to identify additional image artefacts.We select likely stars through cross-matches to Gaia (Gaia Collaboration et al. 2018), where we reject sources with significant (3) evidence for parallax, and to SDSS star/galaxy morphology classifications (Stoughton et al. 2002).
We then flag AGN through matches to catalogued sources in the Milliquas catalogue (Flesch 2021), or via WISE colour cuts (Wright et al. 2010;Stern et al. 2012).We further cross-match to NASA's NED database, to flag any missing catalogued sources (Helou et al. 1991).We seek to distinguish between 'routine' AGN variability and extreme AGN flares.We search for evidence of flaring activity at the time of neutrino detection using the data provided in the ZTF alert packets (Patterson et al. 2019), which are based on difference images.For cases where a source appears to be significantly variable, or may have been flaring at the time of neutrino detection, we run dedicated forced photometry on the science images to produce a source lightcurve (Masci et al. 2019).We reject AGN with no evidence for contemporaneous flaring as 'AGN variability'.After removing those sources flagged as stars ( 17), image artefacts (17) or AGN variability (84), we are left with 54 'interesting candidates' (5 AGN flares, 12 confirmed transients and 37 unclassified sources).The full breakdown in classification is shown in Figure 2.
These interesting candidates include potential transients, which    we seek to classify spectroscopically.Some objects will have already been classified serendipitously, in particular those brighter than 19.0 mag selected by the ZTF Bright Transient Survey (Fremling et al. 2020;Perley et al. 2020).The efficiency with which candidates were classified can be seen in Figure 3. Above a peak apparent magnitude of 19.5, almost all candidates are classified.There were 106 fainter candidates in total, of which 68% were classified.The spectroscopic programs which supported our program are listed in Table 2.
The transients are further broken down by subclass in Figure 4.  Four could be immediately excluded as candidates based on their classification as SNe Ia, a population not predicted to emit high-energy neutrinos.Of the remainder, beyond the two TDEs, no further sources exhibited electromagnetic signatures from the theoretically-predicted neutrino emission scenarios listed above.However, we cannot exclude the possibility that future theoretical work proposes additional neutrino emission scenarios not considered here, and therefore we cannot definitively rule out these transients as neutrino sources.We can state that none of them are consistent with existing models.
A selection of highlighted results is given in the following sections.We specifically list transients which were unclassified at the time of our ToO observations, and the additional follow-up that we took to confirm their nature.ZTF data for two other candidate neutrino sources from the literature, PKS 1502+106 and BZB J0955+3551, are also outlined in Section 5. We omit ZTF data for the probable neutrino-TDEs AT 2019dsg and AT 2019fdr, as these have already been released in dedicated publications (Stein et al. 2021b;Reusch et al. 2022).All other classified sources are listed in Appendix tables B1-B20.
For four neutrino campaigns (IC200107A, IC201007A, IC201222A and IC210922A), no candidates were identified, and there are no corresponding lists in the appendix.

SN 2019pqh and IC190922B
Follow-up of IC190922B by ZTF identified the candidate supernova SN 2019pqh/ZTF19abxtupj (Stein et al. 2019d).The lightcurve is shown in Figure 5, where upper limits are illustrated with triangles.The arrival time of the neutrino on 2019 September 22 is marked with a dotted line, and the supernova is detected in the subsequent ToO observations.The neutrino arrival time was close to optical peak, consistent with a CSM-interaction scenario.
However, a spectrum was taken by the NUTS2 collaboration (Holmbo et al. 2019), and the supernova was classified as a Type II supernova without spectroscopic signatures of CSM interaction (Reguitti et al. 2019).A higher-resolution spectrum of the object was also obtained on 2019 September 28, shown in Figure 6, using the Low Resolution Imaging Spectrometer (LRIS) spectrograph at the Keck observatory (PI: Yan) (Oke et al. 1995).A historical spectrum of the host galaxy, taken by the Sloan Digital Sky Survey (SDSS; Abolfathi et al. (2018)), is also shown in Figure 6  supernova classification, with the best match being a Type IIb supernova (SN 1993J, Barbon et al. 1995) 2 days before peak, also shown in Figure 6.With this redshift, a peak absolute magnitude of −18.6 was derived, atypically bright for such a Type II supernova (see e.g.Lyman et al. 2016).One explanation for this enhanced luminosity could be CSM interaction, through which additional kinetic energy is converted to electromagnetic emission.However, the lack of corresponding narrow line spectroscopic signatures generally disfavours the existence of CSM-interaction, and thus any associated neutrino emission from this object.It is therefore likely that SN 2019pqh is instead unrelated to the neutrino IC190922B.

SN 2020lam and IC200530A
ZTF serendipitously observed the localisation of neutrino alert IC200530A on 2020 May 30, just 10 minutes after detection (Stein 2020e), as part of routine survey operations (Reusch et al. 2020b).Additional ToO observations were then conducted on 2020 May 31 in  and  band, and again on 2020 June 1.During ZTF followup of IC200530A, SN 2020lam/ZTF20abbpkpa was identified as a candidate supernova and potential optical counterpart (Reusch et al. 2020b).Spectroscopic observations were triggered using the NOT/ALFOSC spectrograph on 2020 June 6 (PI: Sollerman), which confirmed SN 2020lam as a Type II supernova (Reusch et al. 2020c) using SNID.This spectrum is shown in Figure 7, alongside the match-  As seen in the lightcurve in Figure 8, the supernova was close to peak at neutrino detection time.The object then rapidly cooled, and thus reddened, as is typical for supernovae.Given the neutrino arrival time, CSM-interaction would be the only viable neutrino production mechanism.However, the spectrum shown in Figure 7 had no narrow lines, and therefore did not provide any evidence supporting such CSM interaction.SN 2020lam was therefore likely unrelated to IC200530A.

SN 2020lls and IC200530A
SN 2020lls/ZTF20abdnpdo was also identified as a candidate supernova on 2020 May 30, during ZTF follow-up of IC200530A (Reusch et al. 2020b).Spectroscopic observations were again triggered using the NOT/ALFOSC spectrograph on 2020 June 12 (PI: Sollerman), which confirmed that SN 2020lls was a Type Ic supernova without broad-line features (Reusch et al. 2020d).This spectrum is illustrated in Figure 9, alongside a matching Type Ic supernova spectrum from SNID mapped to the same redshift (Taubenberger et al. 2006).Given that the supernova had not been detected in alert data prior to the neutrino arrival time, and that it belonged to the subpopulation associated with relativistic jets, SN 2020lls was a candidate for the choked-jet neutrino production model.
However, as can be seen in Figure 10, forced photometry analysis (Reusch 2020) revealed a lower-threshold -band ZTF detection pre- ceding the neutrino arrival.Additionally, modelling of the lightcurve using the MOSFIT software (Guillochon et al. 2018) revealed an estimated explosion date predating the neutrino by a week.In combination, these results disfavoured any supernova explosion origin for the neutrino, suggesting that SN 2020lls was instead unrelated to IC200530A (Reusch et al. 2020d).

AGN FLARE CANDIDATES
While the vast majority of AGN detections from our pipeline were categorised as 'AGN variability', visual inspection revealed five AGN which appeared to possibly undergo optical flaring at the time of neutrino detection.The forced photometry lightcurves of these five flares are shown in Figure 11.We attempt to quantify whether the optical lightcurves of these AGN identify them as candidate neutrino sources.
We can consider possible optical signatures associated with neutrino emission.One scenario is the optical flaring observed for TXS 0506+056 during the detection of neutrino IC170922A (IceCube Collaboration et al. 2018).In particular, the optical apparent V-band magnitude of TXS 0506+056 was observed to increase from 15.0 to 14.5 during the time of neutrino detection, corresponding to a flux increase of >50%, over a period of 50 days, relative to the pre-neutrino baseline.
AGN can also exhibit short-term variability for periods of hours or days, but we caution that the detection of a high-energy neutrino alert is a process that requires a substantial fluence at the IceCube detector, even after accounting for the significant Eddington bias associated with cosmic neutrino detection (Strotjohann et al. 2019).The corresponding neutrino flux that is required is inversely proportional to the duration of neutrino emission, and therefore associating a neutrino detection with a temporary electromagnetic signature lasting hours or days would imply an extremely high average neutrino flux for the duration of that signature.Such highly luminous rapid neutrino flares are not well motivated theoretically, it is therefore unlikely that short AGN flares are indicators of neutrino production.
In contrast, longer-term electromagnetic signatures can serve as tracers for neutrino emission.For example, month-long flaring periods of substantially elevated flux can dominate the neutrino emission of blazars (see e.g.Rodrigues et al. 2021).Very long flares, with durations of years, could also be relevant for neutrino production.However, given the relatively short baseline of ZTF observations, our neutrino follow-up program is not well-suited to identify them.We therefore restrict ourselves to searching for such month-long optical flares, as was observed for TXS 0506+056.
We calculate the median flux for each of the five AGN, and each ZTF filter, in a ±25 day window centered on the neutrino detection.We divide this instantaneous flux by the median flux of the source in that filter over the entire ∼4 year ZTF baseline, giving a proxy for relative optical flare strength.These values are given in Table 3.Of the five AGN, only two (ZTF18aavecmo and ZTF20aamoxyt) had a median instantaneous flux >50% above the baseline median flux.ZTF18aavecmo reached this threshold in both g and r band, while ZTF20aamoxyt reached this threshold only in r.We conclude that the remaining three AGN (ZTF18abrwqpr, ZTF18abxrpgu, ZTF19aasfvqm) do not exhibit substantial neutrino-coincident optical flares, and we therefore find no evidence to suggest they are counterparts to high-energy neutrinos.
ZTF18aavecmo (top panel of Figure 11), cross-matched to source WISEA J170539.32+273641.2, is classified as a likely QSO in the Milliquas catalogue.It underwent a single coherent flare lasting approximately one year, with a peak flux roughly triple the quiescent flux measured by ZTF.It was coincident with neutrino IC200530A, detected during the decay of the optical flare.However, this flare was extremely faint, with a median flux at the time of neutrino detection was   ≈ 5 × 10 −13 erg cm −2 s −1 .This is a factor of 20 lower than the flux observed for TXS 0506+056 during the detection of IC170922A (IceCube Collaboration et al. 2018).
ZTF20aamoxyt (middle panel of Figure 11) is a likely AGN flare, detected coincident with IC200929A.It appears to be spatially consistent with its host galaxy nucleus, and cross-matched to WISEA J015853.53+035126.6.Based on a WISE colour of W1−W2=0.7, it is a possible AGN.The neutrino IC200929A was detected during an extended year-long flare.Much like ZTF18aavecmo, ZTF20aamoxyt at the time of neutrino detection was substantially fainter than TXS 0506+056, with a median flux of   ≈ 3 × 10 −13 erg cm −2 s −1 .
We thus identify no optical AGN flares which resemble the multiwavelength flare of TXS 0506+056 in 2017, from any of our 24 neutrino follow-up campaigns.We emphasize that our results do not preclude a significant degree of neutrino emission from AGN more broadly, but they do disfavour scenarios where the vast majority of astrophysical neutrinos are produced by bright AGN optical flares.There is no tension with scenarios where AGN neutrino emission is not dominated by bright optical flares, for example the 'steady state' AGN neutrino models tested in Abbasi et al. (2022) or scenarios where AGN neutrino emission is correlated only to gamma-ray flares.Constraining such scenarios would require more comprehensive multi-wavelength analysis of AGN lightcurves, incorporating other wavelengths in addition to the optical data presented here.A more systematic study of correlations between ZTF-detected AGN flares and neutrinos, including corresponding chance coincidence probabilities, will be the subject of a future analysis.

CANDIDATE NEUTRINO SOURCES FROM THE LITERATURE
We here provide data on two candidate neutrino sources reported in the literature.However, we caution that none of the objects presented here were selected by our pipeline as ZTF candidates, and therefore are not considered part of our systematic search for neutrino counterparts.We would not claim any such object as a candidate neutrino sources in our neutrino follow-up program, because the chance coincidence probability would be unquantifiable.Any search for additional candidate neutrino sources, beyond those candidates found by our pipeline following ToO observations, would require an independent and unbiased systematic analysis procedure.

PKS 1502+106
The neutrino IC190730A was reported by IceCube in spatial coincidence with PKS 1502+106 (Stein 2019a), a particularly gammabright Flat Spectrum Radio Quasar (FSRQ) with a catalogued redshift of z=1.84 (Albareti et al. 2017) .The object was observed by ZTF as part of ToO observations, and was detected under the ZTF candidate name ZTF18aaqnqzx (Stein et al. 2019c).The blazar had already been repeatedly detected as part of the routine survey operations, with both positive and negative flux changes relative to survey reference images.The blazar lightcurve is shown in Figure 12, using data from science images with the ZTF forced photometry service (Masci et al. 2019).The neutrino arrival time is marked in blue.There was no significant flaring observed for this source coincident with the neutrino.The blazar at this point was dimmer than survey reference images, with the neutrino arriving during a year-long fading, and consequently was not selected by our follow-up pipeline as a possible counterpart.There is thus no evidence from the contemporaneous ZTF data to suggest a causal connection between IC190730A and PKS 1502+106, consistent with data from other observatories which did not see any evidence of short-term flaring (Franckowiak et al. 2020).
Data from the Owens Valley Radio Observatory (OVRO) did reveal that the radio flux was elevated in the months preceding the neutrino detection relative to the decade-long observation baseline, behaviour which has also been claimed for TXS 0506+056 and other neutrino-coincident blazars (Kiehlmann et al. 2019).Comprehensive time-dependent modelling has found that the detection of a neutrino alert from PKS 1502+106 is consistent with the multi-wavelength observations of this object, so a neutrino-blazar association is plausible but likely unrelated to the flaring activity (Rodrigues et al. 2021;Oikonomou et al. 2021).In any case, the new optical data presented here can be used to further constrain such neutrino emission scenarios.

BZB J0955+3551
IC200107A was a high-energy neutrino reported by IceCube (Stein 2020a) which was later identified to be in spatial and temporal coincidence with a blazar undergoing a dramatic simultaneous Xray flare (Krauss et al. 2020;Giommi et al. 2020c).The source BZB J0955+3551 (also known as 4FGL J0955.1+3551 and 3HSP J095507.9+355101),located at a redshift of 0.55703 (Paliya et al. 2020), belongs to the specific subclass of extreme blazars, which are characterised by synchrotron peaks at very high frequencies, which had been proposed as especially promising candidates of high-energy neutrinos (Padovani et al. 2016).
More comprehensive multi-frequency modelling has confirmed that the detection of a neutrino alert from an extreme blazar is plausible, though the simultaneous X-ray flare may not be directly related to the neutrino production (Paliya et al. 2020;Giommi et al. 2020b;Petropoulou et al. 2020).The ZTF lightcurve for BZB J0955+3551 is shown in Figure 13.There is no evidence of any optical flaring on short or long timescales coincident with the detection of IC200107A.

LIMITS ON NEUTRINO SOURCE POPULATIONS
With our program, we did not find any likely candidate counterparts from any population except TDEs.We can consider limits that can be placed on other potential sources of astrophysical neutrinos given the non-detections.These limits will clearly not apply for TDEs, because for this population probable counterparts were detected.If AT 2019fdr is ultimately found not to be a TDE, for example if it instead an SLSN (Pitik et al. 2022), the limits would not apply to that population either.It will however apply to all other populations which would be detected by ZTF, provided at least one detection occurred within the 14 day window after neutrino arrival time.Our search has no requirement that an optical lightcurve peaks after the neutrino detection, so these limits also apply to older/fading transients.
For each neutrino, we can consider the probability that an astrophysical counterpart would be detected.A counterpart could only be detected if a given IceCube neutrino was astrophysical, with this as P astro probability being reported by IceCube in GCN notices as the 'signalness' parameter.For each neutrino that was indeed astrophysical, the source could only then be detected if it lay within the area observed by ZTF.We can estimate this probability, P obs , by assuming that the 90% probability is uniformly distributed across the rectangle reported by IceCube, A IC , such that:  obs = 0.9 ×  ZTF  IC (1) where A ZTF is the area observed by ZTF after accounting for detector chip gaps.
The probability to find an optical counterpart is then given by the joint probability that the neutrino is astrophysical, P astro , that the astrophysical source lay in the observed ZTF area, P obs , and the probability that a given counterpart would be detectable with our program  detectable .The values of  astro and  obs for each alert are given in Table 4.
The detectable probability will depend on the selection efficiency,  det , of our program.This selection efficiency in turn depends on the apparent magnitude of the electromagnetic counterpart.Motivated by our classification efficiency in Figure 3, we assume completeness for objects brighter than 19.5 mag, and a classification efficiency of 68% for objects fainter than this (this assumption is illustrated with the red dashed line in Figure 3).We additionally assume a conservative 95% detection efficiency for sources to be found by our pipeline, if said source was imaged by the camera.Chip gaps in the detector are already accounted for in Equation 1.Because the detection efficiency will decrease as the objects approach the ZTF limiting magnitude of 21.5 for 300s exposures, we neglect objects fainter than 21 mag in our calculation: The fraction of astrophysical neutrino sources that are detected by our program will depend on the properties of a given population.For a power law neutrino spectrum, the differential neutrino particle flux at Earth for a transient population with a redshift-independent luminosity distribution is proportional to: where  is the intrinsic neutrino spectral index such that     ∝  − , R(z) =()/(0) is the normalised source redshift evolution for a population with rate (),   is the luminosity distance and   is the comoving volume (see e.g Murase ( 2007)).By normalising Equation 3, we can derive a probability density function (PDF) for the redshift of detected neutrinos: PDFs for  dist (), calculated using the flarestack code (Stein et al. 2020a) for an E −2 spectrum, are shown in Figure 14 for redshift evolutions from a 'GRB-like' population (Lien et al. 2014) and from a Star-Formation-Rate population ('SFR-like') (Strolger et al. 2015).We also show a blazar-like redshift evolution assuming a 'primarilydensity evolution' (PDE) following the formulation of Ajello et al. (2015), motivated by recent observations favouring such a distribution (Marcotulli et al. 2020).The blazar redshift distribution is expected to vary with underlying source luminosity, so we assume a blazar gamma-ray luminosity of   ≈ 10 46 erg s −1 , motivated by observations of TXS 0506+056 (Padovani et al. 2018), and the best fit parameters for a PDE blazar evolution from Ajello et al. (2015).
It can be seen in Figure 14 that GRB-like populations tend to be at greater distances than SFR-like ones, with GRB-like neutrinos being emitted from a median redshift of z = 1.34, whereas SFR-like neutrinos would have a median distance of z = 0.64 and blazarlike neutrinos would have a similar median of z = 0.57.This has a direct impact on the population properties compatible with our limits, because a neutrino population dominated by nearby sources will generally produce counterparts with brighter apparent magnitudes.
For a given source evolution, the probability of detecting a counterpart will then ultimately depend on the underlying luminosity function of the population.For an absolute magnitude, , the counterpart detection probability is equal to the integrated product of the probability that a counterpart has a given redshift,  dist (), and the detection efficiency of our program for the apparent magnitude, (, ), corresponding to that redshift: The impact of different evolutions and absolute magnitudes can be seen in Figure 15.For sources with an absolute magnitude of −21, our program would be sensitive to counterparts up to a redshift of z ≈ 0.45, beyond which  > 21 so  selection = 0.For an SFRlike evolution, this would correspond to  detectable (−21) = 26%, but for the higher-z GRB-like neutrino distribution, we would instead find  detectable = 16%.For a fainter absolute magnitude of −17, our program would probe a much smaller volume up to redshift z ≈ 0.1, so then  detectable would be 5% and 4% for SFR-like and GRB-like populations respectively.
Combining these values, the joint probability for us to find a counterpart during a follow-up campaign is given by: where  is the fraction of astrophysical neutrino sources with an absolute magnitude equal to or brighter than .The probability that no counterpart was detected in any of our 24 follow-up observations is then given by: The probability of no counterpart detection is given in Figure 16 as a function of .The results of our program strongly disfavour scenarios where all neutrino sources have bright absolute magnitudes.The horizontal dashed line in Figure 16 represents a 10% chance of nondetection, and thus a 90% confidence limit.We can use this threshold to set a limit on the luminosity function of neutrino sources, by choosing the appropriate fraction  such that  no counterpart (,  ) > 0.1 .
These constraints on  () at 90% CL are illustrated in Figure 17, for the two source evolutions.These are generic constraints on the underlying luminosity function of neutrino sources, and are agnostic to the actual nature of the neutrino sources which follow the redshift evolutions.They constrain the aggregate neutrino flux emitted by e.g. a SFR-like population, and thus apply equally well to a composite neutrino flux with e.g.multiple SFR-like neutrino populations.To the best knowledge of the authors, this is the first time generic constraints on the optical luminosity function of neutrino sources have been derived, though a similar procedure has already been used to derive limits from optical searches for counterparts to gravitational waves (Kasliwal et al. 2020).One novel consequence of these general limits are the first observational constraints on the contribution of the brightest superluminous supernova to the diffuse neutrino flux, under the assumption these trace the star formation rate.Objects brighter than an absolute magnitude of −22 can contribute no more than 58% of the total astrophysical neutrino alerts if SFR-like.
During the multi-wavelength flare of TXS 0506+056 which coincided with the detection of neutrino IC170922A, the source exhibited a g-band optical flux of ∼6 mJy/m=14.5(IceCube Collaboration et al. 2018), more than double the mean archival g-band magnitude of 2.7 mJy/m=15.3 in Pan-STARRS1 (Chambers et al. 2016).This would correspond to a transient flare of absolute magnitude M  ≈-26 in difference imaging.Blazar flares such as this can contribute no more than 26% of the total neutrino flux.
It should be noted that none of these values account for the impact of dust extinction, which would shift the curves in Figure 17 rightwards to higher luminosities.However, we do not expect that this would significantly impact the results presented here.Our limits are only valid for a given source evolutions, and would need to be adjusted for alternative ones.For example, if the neutrino source evolution were strongly negative (lower number density at higher redshift), a larger fraction of the flux would arise from local sources.Therefore, our limits would instead be more constraining.It should also be noted that these limits assume that a given transient could pass our selection criteria outlined in Section 2, and therefore do not apply to extremely rapid transients such as GRB afterglows, which peak and fade on timescales ≲1 day.Such objects are not well captured by the ZTF public survey cadence or our typical neutrino follow-up observation cadence, and are unlikely to be detected multiple times in order to pass our selection criteria, so our detection efficiency will be somewhat lower.

CONCLUSIONS
The ZTF neutrino follow-up program coincided with the introduction of the upgraded IceCube alert selection, yielding one unretracted alert every 2 weeks and one ZTF follow-up campaign every 4 weeks on average.The program resulted in the identification of two probable neutrino sources (Stein et al. 2021b;Reusch et al. 2022), and the first limits on the optical luminosity function of neutrino sources.
Though the limits presented here constrain only the very brightest transients such as superluminous supernovae (SLSNe), they will continue to become more stringent over time if no new counterparts are identified.As can be seen in Figure 18, extrapolating our analysis to a neutrino sample that was twice or four times as large would lead to substantially more constraining limits, and will be achieved on the present trajectory with 2 or 6 additional years of observations.
Although the data analysis presented considered candidates detected up to 14 days after neutrino detection, our early real-time counterpart searches generally focussed on counterparts detected in the ToO observations scheduled for the first two nights after neutrino detection.Motivated by the systematic analysis performed here, and to improve sensitivity to time-delayed optical signatures such as neutrino emission from choked jets, we have modified our ToO observation strategy to better cover a range of transient timescales.We now trigger deep 300s in  and  band on the first night of observations to obtain deep upper limits or faint detections, and to additionally yield colour information for any active transient.However, we replaced our second pair of 300s exposures with a series of 30s exposures spread over subsequent nights, to complement the public survey and ensure good coverage of the photometric evolution of candidates.Forced photometry is only possible for images from the public survey after they have been published as part of the regular ZTF Data Releases, but with this ToO monitoring we can perform forced photometry analysis in real time (Reusch 2020).We can also better prioritise spectroscopic follow-up with photometric classification.One shortcoming of the ZTF program thus far has been the relatively poor sensitivity to very rapid transients such as GRB afterglows, owing to the median latency of 12.2 hours to first coverage.We plan to implement automated triggering with ZTF, similar to that operated by other observatories such as ASAS-SN (Necker et al. 2022), enabling low-latency observations for at least some favourable neutrino alerts with appropriate accessibility.Dedicated analysis of low-latency follow-up campaigns would yield more stringent constraints on GRB afterglows as neutrino sources.
The results and analysis presented here can serve as a pathfinder for future triggered neutrino follow-up programs with wide-field instruments.In particular, ToO observations with the upcoming Vera C. Rubin Observatory would offer an unprecedented opportunity to probe neutrino sources to much higher redshifts (Ivezić et al. 2019).Multi-band observation coverage would enable photometric classification of many candidates, substantially extending the classification efficiency presented in Figure 3 to much greater depths.An illustration of this is presented in Figure 19, assuming that the same neutrino sample in Table 1 had instead been observed with the Rubin Observatory.For a comparable 60% classification efficiency down to 24th mag, the corresponding limits on the neutrino luminosity function would be much more constraining for lower magnitudes.However, for very luminous optical sources such as blazar flares, the performance of both surveys for such a neutrino sample would be comparable.Given that there are only expected to be ∼12 astrophysical neutrinos in our sample, observations will never be able to overcome the 90% limit from Poisson counting statistics even if they had a perfect 100% efficiency.Instead, as seen in Figure 18, only larger neutrino samples can enable stricter limits on bright sources.
Beyond optical observatories, similar electromagnetic neutrino follow-up programs are planned for telescopes currently under construction.These include near infra-red (NIR) wavelengths with WIN-TER (Lourie et al. 2020;Frostig et al. 2020), at ultra-violet (UV) wavelengths with ULTRASAT (Sagiv et al. 2014), andin gammarays with CTA (Cherenkov Telescope Array Consortium et al. 2019;Carosi et al. 2021).These new instruments, in concert with the continuation of existing follow-up programs, will enable us to study the dynamic neutrino sky across the entire electromagnetic spectrum.

APPENDIX A: NOT FOLLOWED UP
Those alerts not observed by ZTF are summarised in Table A1.Of those 57 alerts not followed up, the primary reasons were proximity to the Sun (18/55), alerts with poor localisation and low astrophysical probability (16/55) and alert retraction (10/55).The full breakdown of neutrino observations statistics can be seen in Figure A1.

Figure 1 .
Figure 1.Latency between neutrino detection and first ZTF coverage.The median latency time of 12.2 hours is indicated by the vertical dotted line.

Figure 2 .
Figure 2. Breakdown of the classification of 172 candidates selected by our program for visual inspection.

Figure 3 .
Figure 3. Top: Apparent magnitude distribution of candidates selected for visual inspection.Bottom: Classification efficiency as a function of peak apparent magnitude.The red dashed line indicates our step-function approximation of classification efficiency.

Figure 4 .
Figure 4. Breakdown of the 12 identified transients by subclass.
. Both the transient and host galaxy exhibit prominent Balmer lines, highlighted in orange in Figure 6, from which a redshift of 0.134 is derived.A template-matching classification using SNID (Blondin & Tonry 2007) confirms a Type II

Figure 5 .Figure 6 .
Figure 5. ZTF lightcurve of SN 2019pqh.The arrival time of neutrino IC190922B is marked by the dashed blue line.

Figure 7 .
Figure 7. Spectrum of SN 2020lam, taken on 2020 June 6.A similar spectrum, from Type IIP supernova SN 2005cs, is shown for comparison.

Figure 8 .
Figure 8. ZTF lightcurve of SN 2020lam.The arrival time of neutrino IC200530A is marked by the dashed blue line.

Figure 9 .Figure 10 .
Figure 9. Spectrum of SN 2019lls, taken on 2020 June 13.A similar spectrum, of Type Ic supernova SN 2004aw, is shown for comparison.

Figure 12 .
Figure 12.ZTF lightcurve of blazar PKS 1502+106.The arrival time of neutrino IC190730A is marked with the vertical dashed line.

Figure 13 .
Figure 13.ZTF lightcurve of blazar BZB J0955+3551.The arrival time of neutrino IC200107A is marked with the vertical dashed line.

Figure 14 .
Figure 14.PDF for neutrino sources as a function of redshift, for both GRBlike and SFR-like source evolutions.

Figure 15 .
Figure 15.Cumulative counterpart detection probability as a function of redshift.

Figure 16 .Figure 17 .
Figure 16.Probability of detecting no counterpart as a function of absolute magnitude for  =1.The dotted line corresponds to 90% confidence.

Figure 18 .
Figure18.Upper limits (90% CL) on the luminosity function of neutrino sources for an SFR-like evolution that would be derived for a ZTF neutrino sample that was twice (  =48) or four times (  =96) the size of the sample presented here.

Figure 19 .
Figure19.Upper limits (90% CL) on the luminosity function for an SFR-like population with our sample of 24 observed neutrino alert and our classification efficiency (ZTF ), and limits that would be obtained for a comparable neutrino follow-up program with the upcoming Rubin Observatory.

Figure A1 .
Figure A1.Breakdown of the neutrino follow-up program, as of 2021 Dec 31.

Table 1 .
Summary of the 24 neutrino alerts followed up by ZTF since survey start on 2018 March 20.

Table 2 .
Summary of dedicated spectroscopic programs for our neutrino follow-up program.

Table 4 .
Probability of finding a counterpart for each neutrino, assuming counterparts are sufficiently bright to be detected by our ZTF neutrino followup program.