TREASUREHUNT: Transients and Variability Discovered with HST in the JWST North Ecliptic Pole Time Domain Field

The JWST North Ecliptic Pole (NEP) Time Domain Field (TDF) is a > 14 ′ diameter field optimized for multi-wavelength time-domain science with JWST. It has been observed across the electromagnetic spectrum both from the ground and from space, including with the Hubble Space Telescope (HST). As part of HST observations over 3 cycles (the “TREASUREHUNT” program), deep images were obtained with ACS/WFC in F435W and F606W that cover almost the entire JWST NEP TDF. Many of the individual pointings of these programs partially overlap, allowing an initial assessment of the potential of this field for time-domain science with HST and JWST. The cumulative area of overlapping pointings is ∼ 88 arcmin 2 , with time intervals between individual epochs that range between 1 day and 4+ years. To a depth of m AB ≃ 29.5 mag (F606W), we present the discovery of 12 transients and 190 variable candidates. For the variable candidates, we demonstrate that Gaussian statistics are applicable, and estimate that ∼ 80 are false positives. The majority of the transients will be supernovae, although at least two are likely quasars. Most variable candidates are AGN, where we find 0.42% of the general z ≲ 6 field galaxy population to vary at the ∼ 3 σ level. Based on a 5-year timeframe, this translates into a random supernova areal density of up to ∼ 0.07 transients per arcmin 2 ( ∼ 245 deg − 2 ) per epoch, and a variable AGN areal density of ∼ 1.25 variables per arcmin 2 ( ∼ 4500 deg − 2 ) to these depths.


INTRODUCTION
With the successful launch, commissioning, and first year of science observations, the James Webb Space Telescope (JWST) has opened new opportunities for the study of the faint (potentially m AB ≳ 30 mag) variable Universe.After discovering time-varying phenomena (i.e., transients and objects that vary in brightness or position), one would need the capacity to monitor these objects at their astrophysically relevant timescales, to determine the nature of their variability.For example, rapid follow-up of supernovae (SNe) is crucial to determine their types and distances.Flexible follow-up of active galactic nuclei (AGN) allows the timescale of AGN variability (which can range from several days to several years) to be measured, which can be linked to properties of the supermassive black hole (SMBH) at the center of its host galaxy (e.g., Vanden Berk et al. 2004).However, transient follow-up and variability monitoring with JWST become difficult when Sun angle restrictions, power generation, and micrometeoroid mitigation limit when any particular area of the sky can be observed.The Continuous Viewing Zones (CVZs) centered on the eclip-tic poles are the only locations where JWST can observe at any time of the year and avoid these constraints.Jansen & Windhorst (2018) selected a particularly suitable area within JWST's northern CVZ, the JWST North Ecliptic Pole (NEP) Time Domain Field (TDF), to enable the exploration of time-varying phenomena with JWST at high redshifts, but also within the halo of our own Galaxy, and in the extreme outer Solar System at high ecliptic latitudes.The JWST NEP TDF is a >14 ′ diameter field centered on (RA, Dec) J2000 = (17:22:47.896,+65:49:21.54),that was carefully chosen to minimize foregrounds, to be devoid of stars brighter than m K ∼ 15.5 mag, and maximize the observing efficiency for time-domain science with JWST.Its location near the NEP allows for follow-up of transients and variable sources at any cadence and at any time of the year.In addition, the NEP suffers minimal Zodiacal foregrounds, with minimal variation in the course of a year, naturally allowing for more sensitive observations per unit time.While the ecliptic latitude of the JWST NEP TDF is very high (b ecl ≃ 86.2 • ), it is located at intermediate Galactic latitude (b II ≃ 33.6 • ), providing a clear but relatively long sight-line through the halo of our Galaxy.This makes the field useful for deep Galactic time-domain science with JWST as well.
This field serves as an important testbed for advancing our understanding of the time-varying universe at fainter limits than can be reached from the ground at optical and near-IR wavelengths, where transients and variability can hold cosmological significance (e.g., dark energy, cosmic star formation rate, and the evolution of supermassive blackholes).Type Ia SNe are often used as standard candles to provide critical constraints on the Hubble constant, the mass density of the universe, the cosmological constant, the deceleration parameter, and the age of the Universe (e.g., Riess et al. 1998).Core Collapse (CC) SNe are equally important, as they not only provide most of the heavy elements in our universe (e.g., Matteucci & Greggio 1986), but also should trace the rate of instantaneous massive star formation.CC SNe rates and their connection to cosmic star formation rates remains poorly understood, as there seems to be a significant mismatch between observed CC SNe rates and what is expected (see, e.g., Cappellaro 2014).
In addition to studying transient events like SNe, we can uncover new insights into the nature of astrophysical processes that lead to AGN variability.AGN are powered by matter accreting onto a SMBH (Lynden-Bell 1969).All massive galaxies host central SMBHs (e.g., Magorrian et al. 1998;Kormendy & Ho 2013), and irregular or varying rates of infall, as well as turbulence and temperature fluctuations, can cause variations in brightness over time.The study of these variable AGN can provide a better understanding of the complex, unstable processes occurring around supermassive black holes (e.g., Shakura & Sunyaev 1976;Ulrich et al. 1997;Kawaguchi & Mineshige 1999).As AGN are linked to the evolution of their host galaxy (e.g., Gebhardt et al. 2000;Ferrarese & Merritt 2000;Kormendy & Ho 2013), this would also allow a better understanding of the co-evolution of galaxies and their central black holes.
AGN are often identified via X-ray or mid-IR emission, colors, or spectroscopic signatures, however, these methods can miss AGN that are intrinsically faint or lack X-ray emission.Variability in particular provides a unique way to identify AGN that might be missed via these other methods (e.g., Boutsia et al. 2009;Pouliasis et al. 2019).Lyu et al. (2022) explore various AGN selection methods in GOODS-S, including those based on X-ray emission, UV to mid-IR spectral energy distributions (SEDs), mid-IR colors, radio emission, and variability, and find ∼10% of the AGN in this field to exhibit optical variability 1 .We refer to Lyu et al. (2022) and Lyu et al. (2023) for a broader treatment of the role of various AGN selection methods.
Variability proves especially promising for identifying faint AGN, as there appears to exist an anti-correlation between AGN brightness and variability amplitude (Hook et al. 1994;Trevese et al. 1994;Cristiani et al. 1997;Giveon et al. 1999;Vanden Berk et al. 2004;Wilhite et al. 2008;Zuo et al. 2012).In addition, the amplitude or timescale of AGN variability is correlated accretion disk size, and therefore black hole mass (e.g., Xie et al. 2005;Wold et al. 2007;Li & Cao 2008;Wilhite et al. 2008;Burke et al. 2021) and redshift (Cristiani et al. 1990;Hook et al. 1994;Trevese et al. 1994;Vanden Berk et al. 2004).A variability survey can thus provide a more complete census of the AGN population and its evolution over cosmic time.In addition, AGN identified through variability may be spectroscopically confirmed at any later time.
A large number of studies focusing on AGN variability have been published in the last few years (for a recent review, see Lyu & Rieke 2022).At optical wavelengths, AGN have been identified via their variability on scales from weeks to years (Cohen et al. 2006;Morokuma et al. 2008;Villforth et al. 2010;Sarajedini et al. 2003Sarajedini et al. , 2011)).
As part of Hubble Space Telescope (HST) programs GO 15278 (PI: R. Jansen) and GO 16252/16793 (TREA-SUREHUNT; PIs: R. Jansen & N. Grogin) UV-Visible imaging of the JWST NEP TDF was secured with the F275W (0.272 µm), F435W (0.433 µm), and F606W (0.592 µm) filters.Areas of partial overlap between individual visits of the observations allow a first look at 2to 4-epoch object variability and transients in this field.For a detailed description of the observational design of these HST surveys, data reduction, stacked images, and source photometry, we refer the reader to Jansen et al. (2024a,b;in prep.).
While previous HST variability studies like Sarajedini et al. (2003) (Hubble Deep Field North) and Cohen et al. (2006) (Hubble Ultra Deep Field) achieved impressive depths (29.0 and 30.5 mag, respectively), they are limited to approximately ∼6.5 arcmin 2 .Villforth et al. (2010), on the other hand, leveraged the expansive area (≳500 arcmin2 ) of the Great Observatories Origins Deep Survey (GOODS), but to a relatively shallow depth (∼26.0 mag).At the same exquisite angular resolution (full width half max ≲ 0. ′′ 09), TREASUREHUNT observations provide a unique combination of depth (m AB ≲ 29.5 mag) and area (∼88 arcmin 2 ) for a variability analysis.
In this work, we provide the first time-domain study of the JWST NEP TDF at visible wavelengths, presenting 12 transients (supernova candidates) and ∼100 variable candidates (AGN candidates) identified using TREA-SUREHUNT data.In § 2, we briefly summarize the TREASUREHUNT HST data, and data processing specific to the detection of transients and variable objects.
2. DATA AND PROCESSING

HST Observations and Data Processing
Observations from HST GO 15278 were taken between 2017 October 1 and 2019 February 9, and those from GO 16252+16793 (TREASUREHUNT3 ) between 2020 September 25 and 2022 October 31.The former program targeted the central (r ≲ 5 ′ ) portion of the JWST NEP TDF, and the latter an outer annulus to r ∼ 7. ′ 8. Together, these programs provide near-contiguous coverage with the Wide Field Camera 3 (WFC3/UVIS) in the F275W filter (λ c ∼ 0.272 µm) and with the Advanced Camera for Surveys (ACS/WFC) in the F435W and F606W filters (λ c ∼ 0.433 and ∼0.592 µm), with a total area of ∼194 arcmin 2 (ACS/WFC).The full HST coverage of the JWST NEP TDF is shown in Figure 1.
As described in more detail in Jansen et al. (2024a,b;in prep.), both programs used 4-orbit CVZ and pseudo-CVZ4 visits to reach 2-σ limiting depths of m AB ≃ 28.0, 28.6, and 29.5 mag in F275W, F435W and F606W, respectively.These limiting depths were determined using the SourceExtractor (Bertin & Arnouts 1996) MAG_AUTO aperture, within areas with single-epoch coverage (lighter grey in Figure 1), and are the adopted nominal survey depths.
The TREASUREHUNT data were retrieved from MAST and post-processed (after the standard ACS and WFC3 calibration pipeline steps of flatfielding and correction for Charge Transfer Efficiency (CTE) effects) to mask satellite trails, remove residual CTE trails, and fit and subtract the background (Jansen et al. 2024a,b;in prep.).The latter was needed because the observing strategy deliberately used the full duration of CVZ orbits, as well as pseudo-CVZ orbits where HST was allowed to dip closer to the Earth limb than the nominal limit for CVZ orbits.This resulted in excess background light that helped fill charge traps in the F275W exposures with their exceedingly dark natural sky background (e.g., O'Brien et al. 2023), but also rendered the observed background level meaningless.
The background-subtracted images were corrected for geometric distortion and aligned to the Gaia/DR3 (Gaia Collaboration et al. 2016, 2023) astrometric reference grid using the DrizzlePac tasks TweakReg and AstroDrizzle (see Hoffmann et al. 2021).The AstroDrizzle software also flags cosmic-rays, a crucial step when identifying faint variability, where cosmic ray induced signals could cause spurious variability if not properly removed.We expect, and tested (as discussed in more detail in the Appendix), that all cosmic ray and other image artifacts were largely removed through the AstroDrizzle step, as each single 4-orbit visit resulted in 8-9 exposures per filter that were stacked together.
All 8-9 exposures per visit were taken at the same pointing with small ≤0.2 ′′ dither offsets in RA and Dec.For each visit, we drizzle-combined the individual post-processed, rectified and aligned exposures per filter (Jansen et al. 2024a,b;in prep.).For ACS/WFC F606W and F435W, this resulted in 24 rectified and aligned science images (and 24 associated weight maps), and for WFC3/UVIS F275W in 23 images and associated weight maps, each at a pixel scale of 0. ′′ 030 pix −1 and centered on the nominal field center of the JWST NEP TDF.We also use the full mosaics of all drizzle-combined images per filter, generated with the same image dimensions and pixel scale, and aligned to Gaia/DR3 to within 0. ′′ 008 in RA and Dec.
The footprints of individual ACS/WFC (and to a lesser extent WFC3/UVIS) visits partially overlap one another, with a cumulative area of overlap of ∼88 arcmin 2 .These overlapping regions are key to searching for transients and variability within this field.Figure 1 showcases the overlapping regions in progressively darker shades of grey for areas of 2-, 3-and 4-epoch overlap.The time interval between observations in these regions of overlap range between 1 day and 4+ years.Regions of overlap between visits in the mosaic are highlighted in progressively darker shades of grey for areas of 2-, 3-, and 4-epoch overlap.This paper focuses on these regions of overlap, with a total area of 87.94 arcmin 2 , to investigate transients and variable sources.Red circles mark the positions of the 12 identified transients, while blue squares and diamonds mark the 190 unique variable candidates discovered in this study.We expect ∼80 of these sources to be false positives (see § 3.2).Symbol shapes and colors identify the samples from which individual variables were selected, as indicated by the legend at the top of the figure.For comparison, the JWST/NIRCam coverage of this field is outlined in orange.We estimate photometric redshifts ( § 4.3) for objects with both HST and JWST photometry.One visit (jeex02; at bottom right) with a bad astrometric solution was ignored in this work.
We note that while the very small number of epochs and sparse time sampling can catch transients like SNe and variable sources like AGN, they cannot provide a light curve to trace their evolution.Further epochs would provide a more robust identification of variable sources, where we may characterize the scatter in the measured brightness (e.g., Pouliasis et al. 2019;Sokolovsky et al. 2017).
While our primary focus is the identification of time variable objects within the TREASUREHUNT HST data, we opt to incorporate ancillary data to better understand certain objects of interest.Specifically, we incorporate JWST/NIRCam (0.8-5 µm), NuSTAR and XMM-Newton, Chandra, and MMT/Binospec data.

JWST Observations and Data Processing
JWST/NIRCam observations within the JWST NEP TDF were obtained between 2022 August 25 and 2023 May 30 as part of JWST GTO 2738 (PIs: R. Windhorst & H. Hammel) in 7 broadband filters (F090W through F444W) and 1 mediumband filter (F410M) as part of the PEARLS program (Windhorst et al. 2023).Figure 1 shows the JWST observations outlined in orange, highlighting the 4 epochs, corresponding to the four orthogonal "spokes".Each spoke consists of two NIRCam pointings that partially overlap in the center.Each pointing includes a total exposure time of 3.5 hours across all 8 filters, or about 3000 sec per filter.Each spoke has an area of 13. ′ 7, totaling 54. ′ 7 of JWST observations.The data were retrieved from MAST and postprocessed by the PEARLS team using their custom pipeline to mitigate 1/f noise, alleviate wisps in the NIRCam/SW filters, mask snowball artifacts, and flatten the background across read-out amplifier boundaries.The individual post-processed images where rectified and aligned to Gaia/DR3 independent of the HST images, but using nearly identical methods.Full mosaics of the field were created for each filter with a 0. ′′ 03 platescale.The 5 sigma point-source limit is typically between 28.0 and 29.1 mag, depending on the filter.For more details on the reduction, calibration and post-processing, we refer to Windhorst et al. (2023), Robotham et al. (2023), andJansen et al. 2024c (in prep.).In § 4.3, we will use PEARLS aperture photometry of transient and varying objects discovered in areas with both HST and JWST coverage to fit SEDs.

X-Ray Observations and Data Processing
The JWST NEP TDF has also been observed extensively in X-rays with NuSTAR, XMM-Newton, and Chandra.NuSTAR has monitored the field since Cycle 5 (PID 5192,6218,8180,and 9267;PI: F. Civano).For this work we use the Cycle 5+6 data and catalog published in Zhao et al. (2021aZhao et al. ( ,b, 2024)).The 21 NuSTAR observations took place over 28 months (between 2019 September and 2022 January), while the three XMM-Newton observations (part of the NuSTAR cycle 6 program) were taken over 15 months (between 2020 October and 2022 January).The total NuSTAR exposure from Cycle 5+6 is 1.5 Ms, reaching a flux limit of 3.3×10 −14 erg cm −2 s −1 in the 3-24 keV band.We note that this is the deepest NuSTAR field taken thus far.The total XMM-Newton exposure time is 62 ks, reaching a limit of 8.7×10 −16 erg cm −2 s −1 in the 0.5-2 keV band.
We also use the first 1.3 Ms of Chandra observations of the JWST NEP TDF (PID 19900666, 20900658, 21900294 and 22900038; PI: W.P. Maksym) spanning cycles 19-23.Chandra observed the field 46 times between 2018 April 12 and 2022 September 25 (a span of 4.6 years), with exposures ranging from 9 to 59.9 ks.Beginning 2021 August 29, these occurred in 90 ks epochs every 3 months, typically broken into 2 or more individual visits (OBSIDs) spanning no more than 30 days per epoch.We reprocessed and reduced the data using ciao (Fruscione et al. 2006), corrected the cross-observations using source catalogs generated by wavdetect, and astrometrically registered them to Pan-STARRS (Chambers et al. 2016).We then used merge_obs to generate combined event files and exposure maps.The faintest 3σ detection in the catalog is 1.35 × 10 −16 erg s −1 cm −2 in the full 0.5-7 keV band assuming galactic absorption and a Γ = 1.4 power law.

Ground Based Spectroscopy
MMT/Binospec (Fabricant et al. 2019) observations of the JWST NEP TDF in 2019 September (PI: C. Willmer) used a 270 lines/mm grating and a 1 ′′ slit width to yield spectra spanning ∼3900-9240Å sampled at 1.30Å/pix and with a resolution of ∼4.9Å (full width half max) at 6500Å.The exposure time was 1 hr (4×900 s), reaching S/N ∼ 5 at m AB ∼ 22.5 mag for quiescent galaxies, and a magnitude deeper for galaxies with significant emission lines.

METHODS TO IDENTIFY TRANSIENTS AND VARIABILITY
In this section, we detail our methods for identifying transients and variability in the 55 areas of overlap (shown as darker shades of grey in Figure 1) with a combined area of ∼88 arcmin 2 .
Transients are usually defined as objects that only appear for a short duration.For this work, where we expect most transients to be SNe, we consider relatively bright point sources that are present in one epoch, yet do not appear in another.We searched for transients using difference images, and considered any source detected therein a transient candidate.In contrast, variable sources tend to show less extreme differences in brightness and are not as evident in difference images.Variable sources are therefore identified after measuring the source flux and its associated uncertainty within a small aperture.Consequently, some sources we identify as transients may be variable sources with extreme amplitudes of variability (e.g., quasars), and some variable sources may be dim transients or transients caught well away from peak brightness (e.g., SNe).Nonetheless, we assume most transients we discover will be SNe, while most variable sources will be weak AGN.However, the methods presented here may identify a variety of objects, including variable stars, tidal disruption events, gamma ray bursts, or quasars.We did not identify any objects that showed appreciable movement between visits.For all objects discovered here, further analysis is required to determine their true nature.

Transient Search
Transients were identified by visual inspection of difference images, each of which was created as follows.First, we bring the drizzled images of individual visits to the same pixel grid using their Gaia-aligned WCS and the Astropy (Astropy Collaboration et al. 2013, 2018, 2022) affiliated package reproject.Then, we generate two distinct difference images for each overlapping region in the deeper F606W images, i.e., 110 difference images for the 55 overlapping areas.For instance, if Visit 1 overlaps with Visit 2, we generate a difference image by subtracting Visit 2 from Visit 1, and another difference image by subtracting Visit 1 from Visit 2. In each difference image we then visually identify sources with positive flux that are consistent in appearance with bearing the imprint of the point spread function (PSF).We did not use the F435W images for our transient search, because they are ∼0.9 mag shallower.
Although rare, bright cosmic ray (CR) induced signals can sometimes survive in final drizzled products and subsequently appear in difference images.We inspected the individual pipeline-calibrated, flat-fielded exposures at the corresponding position of each transient candidate to ensure they appear in each exposure, bear the imprint of the PSF, and are not (residuals of) CRs or other spurious noise.We originally found 20 transient candidates, of which 8 were confirmed to be due to CRs and were removed from further consideration.
We measured the apparent F606W magnitudes of each transient using the difference images.To measure transient magnitudes, we centered a circular aperture with a radius of 8 pixels on the transient in the difference image, and measured its magnitude within that aperture.These apertures are ≳5× the full width half max (FWHM) of the PSF and encompass the transient signal.We used the Python photometry package Photutils (Bradley et al. 2022) and specified the exact center of each aperture.We created a difference error map by adding the weight (WHT) maps from the individual visits in quadrature, then taking ERR = 1/ √ WHT.We direct this error map into Photutils to derive an uncertainty on the reported flux.Magnitude errors are then calculated as ∆m = (2.5/ln 10) • σ F /F , with F the measured flux, and σ F the uncertainty thereon, both as reported by Photutils.
For each transient, we also performed matchedaperture photometry on a subset of F435W difference images, generated for those regions of overlap that contained a transient.

Variability Search
To detect variability in galaxies, we first constructed source catalogs using SourceExtractor, with object positions and aperture photometry for each visit.For initial object detection, as summarized in Table 1, we require an object to be detected at the 1.5σ level above the local sky level (given a gain value corresponding to the median exposure time of the individual exposures in that visit and filter) after smoothing with a Gaussian kernel with a FWHM of 3 pixels, and require a minimum number of 5 contiguous pixels to minimize detection of spurious sources.For the same reason, we also generously masked the pixels along the detector borders, where fewer contributing exposures result in excess noise and imperfect rejection of CR signal.
We compared brightness measurements across distinct visits to identify potential variables.We specifically focus on detecting variability in the F435W and F606W filter observations, as the F275W filter magnitude limit is substantially brighter, the total area of overlap between individual F275W visits is substantially smaller, and the WFC3 F275W observations are not contemporaneous with the ACS/WFC F435W and F606W ones.

Photometry
Our goal is to leverage HST's high resolution to isolate the potentially varying nuclei within galaxies and SNe near their core while ignoring the non-varying extended component.To achieve this, we first set the minimum contrast for deblending in SourceExtractor (DEBLEND_MINCONT) to 0.001.This value is lower than is typical, but allows isolating the nuclei of galaxies and potentially locate variability in bright regions that may not be in the very center of a galaxy.In addition to a low Note.For each SourceExtractor parameter in the first column, the second an third columns list its value for identifying variable sources and for full aperture photometry, respectively.We used a circular aperture with a diameter of 0. ′′ 24 diameter for variability identification, specified using the PHOT_APERTURES parameter, and the FLUX_AUTO aperture for full aperture photometry.We set the GAIN value to the median exposure time of the individual exposures combined per visit.
deblending threshold, we opted to use apertures close to the FWHM of stellar objects.For our resolution of ≲0.′′ 09, we used a circular aperture with a diameter of 0. ′′ 24 (8 pixels).For apertures much larger than this the surrounding galaxy tends to dilute any variability of a compact central source.We henceforth refer to this small aperture as the 0. ′′ 24 aperture, and to the corresponding magnitudes as m 0.24 .The radius of this aperture is shown as the dashed blue vertical line in Figure 2 to facilitate comparison with detected stellar objects.
With DEBLEND_MINCONT = 0.001, the measured variability sometimes appears offset from the centers of galaxies.This could be due to off-center SNe or SMBHs that have yet to settle onto the centers of their host galaxies.To measure host galaxy properties for these offset events, we also ran SourceExtractor with DEBLEND_MINCONT = 0.06, which reduces deblending of galaxies into individual clumps.We thus have two types of aperture photometry for each source: 1) within the small circular 0. ′′ 24 diameter aperture used to identify variability, and 2) for the full Kron aperture, as identified by SourceExtractor.Table 1 summarizes the relevant SourceExtractor parameters used for these two types of measurements.For our matched aperture F435W photometry, we used SourceExtractor in dual image mode, employing the full F606W mosaic as the detection image.
To identify variability, we need to first carefully match the objects from one visit catalog to another.We did so using the Astropy match_coordinates_sky function, which finds the nearest object in the catalog from one visit to a given object in a catalog from another.To ensure slight differences in position do not contribute to the measured variability, we only considered two sources to be a match across two different visits if they are within one pixel of each other (0. ′′ 03).We increased this threshold to 0. ′′ 08 when matching the F606W catalog to the F435W catalog, to account for the fact that objects may have different PSF and surface brightness distributions at different wavelengths.Matches between the DEBLEND_MINCONT = 0.001 and DEBLEND_MINCONT = 0.06 catalogs were identified using a search radius of 3 arcsec, since the highly deblended catalog may find variability in the outskirts of host galaxies, and thus the positions may be offset by much more than our astrometric uncertainties and pixel size.
We opted to ignore stars in our variability analysis due to uncertainties in the ACS PSF (e.g., Villforth et al. 2010).We flagged as stars all objects with a FWHMt less than 0. ′′ 15 or m AB < 17 − 4 log 10 (FWHM) mag.Grey shading in Figure 2 indicates the stars and faint point-like sources that were removed from our analysis.In addition, we ignored objects with a SourceExtractor FLAGS5 > 2, including objects with saturated pixels, truncated footprints, corrupted apertures/footprints, or other issues.Last, we also ignored all objects within the jeex02 visit in the F606W filter, as it has a poor astrometric solution that hampers reliable alignment and variability measurements (Jansen et al. 2024a,b;in prep.).Some variable candidates, like variable stars and ultrabright AGN like quasars, will be point-like and excluded by our methods.Missing quasars in particular may bias our final AGN estimates.Nonetheless, these objects are relatively rare.In addition, our methods would miss any AGN that does not vary on the timescales we probe.Variability identified in this paper is therefore the minimum amount of variability expected in this field.

Variable Source Selection
Sources were identified as variable if they varied by more than 3× their combined photometric uncertainty.Therefore, careful consideration of uncertainties is crucial to ensure variable sources are robustly identified.Combining and resampling multiple lower-resolution images into a higher-resolution composite image (using, e.g., the drizzle algorithm) causes both the signal and the noise in adjacent pixels to become correlated.Photometric uncertainties are known to be underestimated when correlated pixel noise is not taken into account (e.g., Casertano et al. 2000;Labbé et al. 2003;Blanc et al. 2008;Whitaker et al. 2011;Papovich et al. 2016).In addition, the ACS PSF is known to vary due to telescope breathing (e.g., Villforth et al. 2010), which could potentially contaminate variability identification performed within small apertures.Finally, cosmic rays may cause changes in brightness if not filtered out properly.We employed an empirical approach to statistically account for correlated pixel noise and PSF changes, as well as other sources of error, regardless of their origin.
Figure 3 highlights our method to best calibrate photometric uncertainties, which is based on the variability search methods of Cohen et al. (2006).Measured magnitudes were compared for all objects that appear in more than one visit.For each such object, Figure 3 shows the magnitude difference, ∆m 0.24 , versus the mean magnitude, m 0.24 .We used the distribution of these magnitude differences to calibrate the empirical (true) photometric scatter, where the inner brown curves contain ∼68.26% of the data points in each magnitude bin (shown as black circles and curves).To account for correlated pixel noise, a varying PSF, and other unaccounted-for sources of error, we multiplied all SourceExtractor magnitude er-rors by a fixed scale factor (1.15 for F606W and 1.35 for F435W) to best match the observed distribution of photometric scatter.The scaled result is represented by the large brown circles.Magnitude difference uncertainties (for F606W) were calculated as: with N the number of sigma, F 1 (F 2 ) the measured flux for Visit 1 (Visit 2), and σ 1 (σ 2 ) the corresponding scaled flux error.
We emphasize that Figure 3 is used to calibrate photometric uncertainties and not to identify variable candidates.The blue squares represent objects flagged as variable at ≥3σ ∆m in F606W.These objects are determined to be those which vary by more than 3× their individual photometric uncertainty.In general, these sources fall outside the brown 3σ curves in Figure 3, but this isn't always the case.Some objects outside the brown 3σ curves are not flagged as variable if their individual photometric uncertainties are larger than the median.Similarly, some objects inside the brown 3σ curves are flagged as variable if their individual photometric uncertainties are smaller than the median.
In general, we do not expect zeropoint offsets between visits to be a concern, as changes in zeropoint over the duration of the TREASUREHUNT program are <1% for ACS/WFC (Bohlin et al. 2020;O'Brien et al. 2023).The histogram on the right side of Figure 3 is centered at ∆m 0.24 = 0.02 mag (with a standard deviation of 0.25 mag that is dominated by the photometric uncertainties of faint sources), and therefore shows that the zeropoints between detectors and different regions on the same detector as calibrated by the standard HST pipeline are correct to ≲0.02 mag.
We only considered objects brighter than m 0.24 ≃ 28.6 mag in F435W or m 0.24 ≃ 29.5 mag in F606W (corresponding to the predicted 2σ magnitude limits for sources in the field) in at least one epoch.We manually sorted through variable candidates, and removed 63 sources that are either contaminated by diffraction spikes from a bright star, or are close to the edge of a detector.In a vast majority of these cases a diffraction spike overlapped the object in one visit but not in another, observed at a different orientation.
We created a variability catalog for objects that vary more than 3σ ∆m in F435W, and similarly a variability catalog for objects that vary more than 3σ ∆m and 5σ ∆m in F606W.In addition, we created catalogs with sources that vary (in the same direction) in both filters with 2σ ∆m and 3σ ∆m significance.The significance of the catalogs for objects that vary in both filters to 2σ ∆m is √ 2 (2σ ∆m ) ≃ 2.83σ ∆m , and to 3σ ∆m is We plot the difference in measured F606W magnitude, ∆m0.24, versus the average measured F606W magnitude, m0.24, for objects observed more than once in regions of overlapping HST coverage, where ∆m0.24 is defined in the sense Visit 2 − Visit 1.Each small grey point is an individual object, as detected and measured with SourceExtractor.The large brown circles depict SourceExtractor photometric errors, scaled by a factor 1.15 to match the observed distribution of magnitude differences (shown as small black circles) measured in 0.5 mag wide bins for 25.5 < m0.24 ≤ 28.0 mag such that ∼68.3% of all data points fall between them.The dashed black curves represent the inner ∼68.3% range of the distribution at each magnitude and are polynomial fits to the black points.The brown curves therefore represent median 1, 2 and 3σ ranges as fit to the brown circles.Blue squares identify objects that vary at ≥3σ ∆m in F606W, where σ ∆m is defined in Equation 1.The histogram on the right-hand side illustrates the distribution of magnitude differences for all objects detected in the field.In grey, we list mean, median, standard deviation, and median absolute deviation of the y-axis distribution.
√ 2 (3σ ∆m ) ≃ 4.24σ ∆m .We consider the latter sample (variability detected in both filters to 3σ ∆m ) to be our most conservative, robust sample, as objects that vary in both filters are most likely to be true variables, although the shallower F435W images severely limit their number.

Accounting for False Detections
Since variable sources are identified using 2σ ∆m or 3σ ∆m significance levels, we can expect false detections (e.g., noise peaks).We herein refer to the individual variable sources identified in this work as variable candidates, since we cannot determine which individual sources exhibit genuine variability and which are spurious detections.
We predict the number of false detections using Gaussian statistics.As demonstrated in Figure 4, we observe a Gaussian distribution in ∆m 0.24 for all 4,059 galaxies with 28.0 < m AB < 28.5 mag.We confirm that brightness bins of 25.0 < m AB < 27.0 (2,550 galaxies), 27.0 < m AB < 28.0 (4,892 galaxies), and 28.5 < m AB < 29.0 mag (5,336 galaxies) similarly follow a Gaussian distribution.
Given a Gaussian distribution, we expect a false detection rate of 0.27% (assuming a two-tailed test, where the difference in magnitude can either be positive or negative) for variable candidates identified at 3σ ∆m significance.That rate falls rapidly for levels higher than 3σ ∆m .We also identify variable candidates that vary in both filters (in the same direction) with 2σ ∆m and 3σ ∆m significance.Since we implement the requirement that the variability must be of the same sign in both filters, we estimate the amount of false detections using a one-tailed test (i.e., 0.135% false detections at 3σ ∆m ) for these samples.
To demonstrate the methodology by which we assessed the number of false positives as a function of the significance of the detected variability, Figure 5 shows the number of detected variable candidates, the expected number of false positives (and percentage of the total number of objects in the parent sample), and the number of genuine variable sources in 0.5σ ∆m wide bins for the ) for all galaxies (4,059 total) between 28.0 < m AB < 28.5 mag.The black circles show the number of galaxies (N ) for each ∆m0.24 bin, where the error bar is √ N .The solid black curve is a Gaussian fit to the data.The vertical blue lines represent the 3σ thresholds based on the Gaussian fit.This proves that the distribution of ∆m0.24 follows a Gaussian distribution, such that the false detections can be assumed to follow Gaussian statistics (see Section 3.2.3).
F606W sample that varies by more than 3σ ∆m and the F435W+F606W sample that varies by more than 2σ ∆m .
In addition to false detections due to Gaussian statistics, image artifacts like cosmic rays, hot pixels, bias stripping, amplifier offsets, CTE trails, dust motes, optical ghosts, and cross-talk, may cause artificial variability to be detected.Of these, cosmic rays are the brightest and most prevalent, yet will not affect our results due to robust flagging and having 8 exposures (per pointing, filter and epoch) to stack using the AstroDrizzle software.We explore this in more detail in the Appendix.The other image artifacts are largely removed by the standard HST pipeline.Therefore, we do not expect and find no evidence for image artifacts to bias our results.Nonetheless, confirmation of which individual sources genuinely vary will require follow-up observations.

RESULTS
In this section, we present our findings regarding the identification and characterization of transients and variable objects in our data set.We detected a total of 12 transients, including two quasar candidates, and 190 unique galaxies exhibiting variability (where ∼80 are expected to be false positives).Figure 1  The number of genuine variables and false positives add up to the total number of candidates in each bin.The red labels above each bin show the percentage of all objects in the field that are expected to be false positives for that σ ∆m -bin.The top panel is for sources that varied by more than 3σ ∆m in F606W, and the bottom panel for sources that varied by more than 2σ ∆m in both F606W and F435W, and did so in the same sense.Hence, assuming Gaussian statistics, the upper panel reflects a two-sided test, and the lower panel a one-sided test.In total, we estimate ∼80 false positives in the sample of sources that varied by more than 3σ ∆m in F606W, the sample of sources that varied by more than 3σ ∆m in F435W, and the sample of sources that varied by more than 2σ ∆m in both F606W and F435W.
diamonds).Variable candidates are identified in almost all areas of overlap, except those with insufficient astrometric precision or excess noise.We note that the techniques used to identify transients and variable candidates (as described in §3) cause some objects detected as variable candidates to also be classified as transient candidates (e.g., the quasars).Rather than appearing and then disappearing, these actually vary in brightness over time with a sufficiently large amplitude that they were detectable in our difference images.

Transients in the JWST NEP Time Domain Field
We identified a total of 12 transients, for each of which we show 3 ′′ ×3 ′′ cutouts of the F606W images from the relevant visits and of the corresponding difference image in Figure 6.For each of these transients, Table 2 lists the coordinates, dataset IDs and date of observation of each visit, time interval between epochs, as well as the measured apparent brightness (difference in magnitude between the two epochs) in F606W.We also report the F435W magnitude, if the matched-aperture flux exceeded the F435W detection limit.We call the epoch in which the transient signal appears "Visit 1".The transients range in brightness from m AB ∼ 24.8 mag to ∼27.5 mag in F606W.This may not correspond to the peak brightness of the transient, as the measured brightness depends on when the transient event image was taken during its respective light-curve.At magnitudes fainter than m AB ∼ 27.5 mag, the difference images become dominated by image noise, preventing reliable identification of transients.
In all, we detected ∼0.14 transients per arcmin 2 of overlapping area, or ∼491 deg −2 .Since two epochs are needed to detect a transient, the areal density per epoch is therefore ∼0.07 transients per arcmin 2 (∼245 deg −2 ), or ∼0.77 per ACS/WFC footprint per epoch.For comparison, Dahlen et al. (2012) searched for SNe in the GOODS fields (using the ACS/F850LP filter), and identified 118 SNe in a total area of 1 deg 2 to m F 850LP ∼ 26 mag.If we only consider the 4 transients with m F 606W ≲ 26.0 mag, we find an areal density of ∼164 deg −2 , consistent with Dahlen et al. (2012) when factoring in differences in detection filter, SNe colors and redshift distribution, and small number statistics.We also note that it is likely that not all our transients are SNe.

Transients with Significant X-ray Detections
Next, we checked whether any of the transients were detected in X-rays within XMM-Newton, NuSTAR, or Chandra observations of the field.We matched the transient positions with published and preliminary X-ray source catalogs, and found that 3 of the transients are noticeable X-ray emitters.While an in-depth analysis of the X-ray observations of these sources is beyond the scope of the present paper and is deferred to future work, we present the measured XMM-Newton, NuSTAR, and Chandra fluxes in Table 3.
For the XMM-Newton and NuSTAR catalogs, we used a matching radius of 5 ′′ and 30 ′′ , respectively, and found a likely match for transients T7, T10, and T11.We inspected the Chandra 0.5-7 keV event-file images for obvious sources, and here also identified sources consistent with T7 and T10.For the other transients, we extracted source counts from a region with a radius of 2 ′′ and background counts from a nearby source-free region with a radius of 10 ′′ , but found no evidence for any X-ray counterparts above 1σ significance.Chandra upper limits vary across the field due to detector geometry, as well as vignetting and PSF spreading for off-axis sources, but we expect typical 3σ upper limits of ∼1.4×10 −5 ct s −1 for the full 1.3 Ms exposure.For T7 and T10, we then extracted source fluxes using the ciao tool srcflux on the cross-registered event files for each OBSID.To convert count rates to fluxes, we used a model with assumed Galactic absorption and a Γ = 1.4 power law.At r    ≃ 6.5 ′ from the field center (characteristic of T7 and T10), this corresponds to 3σ upper limits of 3.2×10 −16 erg s −1 cm −2 .Note, however, that the epochs of individual X-ray observations do not necessarily match or overlap with specific HST visits in which a transient was observed or absent.

Notes on Individual Transients
Transient T1 is detected only in F606W (although it is faintly discernible also in F435W) and appears in or superimposed on the outskirts of a diffuse, possibly clumpy disk galaxy, ∼1.′′ 08 S and 0. ′′ 16 W of the estimated galaxy center.It has a red (F435W −F606W ) color compared to the galaxy in a color composite image, and its point-like morphology and complete absence in Visit 2 (1.06 yr earlier) suggest it is a supernova (SN) within this m F606W ∼ 23.9 mag host galaxy.
Transient T2 is detected in both F606W and F435W and appears off-center in a centrally concentrated galaxy (possibly an elliptical galaxy) or in its fainter and partially overlapping apparent companion galaxy.It is located ∼0. ′′ 14 E and 0. ′′ 09 N of the core of the brighter (m F606W ∼ 24.5 mag) galaxy and closer in projected distance (∼0.′′ 07 W and 0. ′′ 10 S) to the center of the fainter one (not separately photometered by SourceExtractor).Its point-like morphology and complete absence in Visit 2 (0.30 yr earlier) suggest it is a SN within either galaxy, with a visible color similar to that of its host.
Transient T3 is detected only in F606W.Its very red color is similar to that of a nearby extended disk galaxy or grouping of galaxies that is barely visible in F435W.Its point-like morphology and complete absence in Visit 2 (0.31 yr earlier) suggest it is a SN associated with this m F606W ∼24.0 mag host galaxy, appearing ∼1.′′ 11 W and 0. ′′ 32 S of the center of the brightest region, although the relatively large separation from the main body of this potential host may allow different interpretations.
Transient T4 is detected in both F606W and F435W within a m F606W ∼ 23.0 mag disk galaxy with knots of apparent star formation and a secure MMT/Binospec spectroscopic redshift of z = 0.615 .It appears ∼0.′′ 10 N and 0. ′′ 10 E (a projected distance of ∼1.0 kpc) of the estimated galaxy center, and its color is similar to that of the host.T4 appears point-like and is completely absent in Visit 2 (0.90 yr earlier), but is much fainter than its host galaxy, suggesting a SN caught either before or significantly past its peak brightness.T4 is presently the only transient with a spectroscopic redshift.
Transient T5 is detected in both F606W and F435W and appears well outside a nearby, highly inclined disk galaxy.It is located ∼0. ′′ 10 E and 0. ′′ 01 N of the core of that bright m F606W ∼ 22.5 mag galaxy, with a color that is similar or slightly bluer.No hint of this point-like source is discernible in Visit 2 (1.18 yr later), suggesting that T5 may be a SN associated with the disk galaxy, although the relatively large separation from this potential host allows other interpretations.
Transient T6 is detected only in F606W (with perhaps a hint discernible in F435W) within a red and diffuse disk galaxy, ∼0.′′ 09 E and 0. ′′ 09 N of the estimated galaxy center.It has a red color consistent with that of the galaxy, and its point-like morphology and complete absence in Visit 2 (1.18 yr earlier) suggest it is a SN within this m F606W ∼ 24.4 mag host galaxy.
Transient T7 is detected in both F606W and F435W, appears point-like and isolated.It lacks any nearby potential host galaxy, although there may be a hint of signal ∼0.′′ 59 due E and due W, and ∼0.′′ 23 due S, as well as some "fuzz" ∼0.′′ 42 to the NW in both Visits that could be due to a very faint potential host or galactic companions.Unlike most transients reported here, T7 is also clearly detected in Visit 2 (2.24 yr later) and has significant X-ray detections by both Chandra and XMM-Newton (see Table 3).This transient is thus unlikely to be a SN, but rather a quasar that exhibits significant variation in brightness.
Transient T8 is detected only in F606W and is our faintest transient at m F606W = 27.57mag.It appears ∼0.′′ 27 E and 0. ′′ 04 S of a faint (m F606W ∼ 27.3 mag) and small nearby galaxy that may be similar in color, although possibly not quite as red.This point source is absent in Visit 2 (1.02 yr earlier), and the proximity to a potential host galaxy suggests that it could be a SN.
Transient T9 is detected only in F606W (although it is faintly discernible also in F435W) and appears in what looks like either a spiral arm of a relatively face-on disk galaxy or in an interacting galaxy pair6 , ∼0. ′′ 20 E and 0. ′′ 08 N of the estimated galaxy center.Its color is redder than that of an adjacent clump in the same spiral arm, but is consistent with the color of most of the disk.Its point-like morphology, projected location within a spiral arm, and absence in Visit 2 (3.87 yr later) suggest that T9 is a Type II SN associated with massive star formation within this m F606W ∼ 23.4 mag host galaxy.
Transient T10 is detected in both F606W and F435W, appears point-like and is isolated, lacking any nearby potential host galaxy.Unlike most transients reported here, T10 is also clearly detected in Visit 2 (0.30 yr earlier) and has significant X-ray detections by both Chandra and XMM-Newton (Table 3).This transient is thus unlikely to be a SNe, but rather a quasar that exhibits significant variation in brightness.
Transient T11 is detected in both F606W and F435W, appears point-like and is isolated.The nearest galaxies are ∼2.′′ 3 to the S and ∼2.′′ 6 to the N, both >3 galaxy diameters away.Its color is relatively blue, allowing it to be clearly visible in F435W though m F606W = 27.55 mag.This point source is absent in Visit 2 (0.30 yr later).Interestingly, T11 was detected with NuSTAR in the 3-8 keV band, and thus is unlikely to be a SNe.
Transient T12 is detected only in F606W and appears adjacent to a stellar bar across the center of an extended face-on spiral disk galaxy, ∼0.′′ 12 S and 0. ′′ 02 E of the estimated galaxy center.Its color is redder than that bar, but consistent with the color of other portions of the faint disk.Its point-like morphology, projected location within the disk close to the central portions, and absence in Visit 2 (0.89 yr earlier) suggest that T12 is a SN within this m F606W ∼ 22.5 mag host galaxy.
The four brightest transients, T1, T2, T7, and T10, all of which are brighter than 25.6 mag, may be bright SNe or quasars.The significant X-ray detections of T7 and T10 specifically argue for the latter.In contrast, one of the faintest transients in this study, T11 with m AB ∼ 27.5 mag, is especially interesting as it is the only transient without a potential host detected in F606W, yet with noticeable X-ray emission.

Variability in the JWST NEP Time Domain Field
We identified 190 unique candidate variable candidates that meet the selection criteria of various samples discussed in § 3.2, of which we estimate ∼80 to be false positives (the exact number depends on how many false detections appeared in multiple variability samples).Fig- ure 1 shows the locations of all 190 variable candidates.Symbol shapes (squares and diamonds) and hues identify the sample in which each variable was identified.Table 5 (at the end of this paper) provides a full list of IDs, celestial coordinates, dataset IDs and date of observation of each visit, as well as brightest measured apparent magnitude for each variable candidate.We also list the change in brightness between visits (in the sense Visit 2 − Visit 1), and the significance of variability in units of σ ∆m , for both the F435W and F606W filters.
We emphasize that we can not claim any individual variable candidate varying at < 5σ ∆m significance to be a genuine variable.We can only place strong constraints on the overall number of variables in the area samples in a statistical sense.Conversely, we also do not exclude any specific sources as false detections from this final variability catalog, because we can only account for false positives in a statistical sense.We acknowledge that when an individual source is of interest, the source should be carefully analyzed and may require additional observations to ensure it is a genuine variable.
More relevant for population statistics, Table 4 lists the inferred number of genuine variable candidates, factoring in false detections, for different significance levels in either the F606W or F435W filters, or both.Figure 7 shows 3-color composite cutouts for a representative sample of 12 variable candidates.The bottom row shows our most conservative sample of 3 galaxies that vary in both F435W and F606W at ≥3σ ∆m significance.Examples of variable sources drawn from the less restrictive samples are shown in the first three rows.
We estimate the total number of variable sources using the 1st, 3rd, and 4th rows of Table 4, since variables in the 2nd and 5th rows have more stringent selection criteria and are already included within these less restrictive samples.There are 206 objects listed in these specific rows of Table 4, where 16 are duplicates (i.e., matching more than one of these sample selection criteria), leaving 190 unique variable candidates.With an estimated ∼80 false positives, this results in ∼110 unique variable sources, or 0.42% of all 26,468 detected sources in the field.This is the maximum variable source fraction inferred in this work, as these samples are not necessarily independent.Our most conservative sample of ≥3σ ∆m variable objects in both F435W and F606W only makes up 0.02% of the parent population.
We also estimate the areal density of variable sources using the total number of inferred genuine variables.These ∼110 unique variable candidates translate to ∼1.25 variables per arcmin 2 (∼4500 deg −2 ) of overlapping area.We can therefore expect ∼14 new variable sources per additional fully overlapping ACS/WFC footprint (∼11.3  1) the sample of variable candidates, with the significance limit given between parentheses, (2) the total number of sources in the parent sample, limited to the F435W parent sample in the case of a combined F606W+F435W selection, (3) the number of objects satisfying the criteria at face value, (4) the number of false positive detections expected, (5) the number of genuine variable sources after statistical correction for false positives, and ( 6) the inferred percentage of sources in the parent sample that are variable.
arcmin 2 ), assuming F606W imaging to similar depths as the TREASUREHUNT observations.The 3σ ∆m threshold, which determines by what magnitude a galaxy must vary to be flagged as variable, is key to putting our findings in context with other work.Higher noise levels naturally result in fewer detected galaxies, so the threshold helps to normalize these variations.The solid brown curves in Figure 3 show the 1σ ∆m , 2σ ∆m , and 3σ ∆m thresholds (for the F606W filter) as a function of magnitude.For all sources with m 0.24 ∼ 26.0 mag (a 1.0 mag wide bin centered on 26.0 mag), the median 3σ ∆m threshold is ∆m 0.24 = 0.14 mag.In other words, a typical 26.0 mag source would need to vary by at least 0.14 mag to be flagged as variable for this specific data set.For bins centered on 27.0 and 28.0 mag, the 3σ ∆m thresholds are 0.24 and 0.41 mag, respectively.For the F435W filter, the 3σ ∆m thresholds for bins centered on 26.0, 27.0, and 28.0 mag are 0.26, 0.44, and 0.82 mag, respectively.
Although our variability identification pipeline was developed to identify variable AGN, we still detect 6 out of the 12 transients presented in Section 4.1.Specifically, it identifies transients T2, T3, T5, T6, T9, and T11.Transients T7 and T10 were originally identified to exhibit variability, but were flagged and removed as star candidates due to their small FWHM.Transients T1, T4, T8 and T12 did not satisfy the 3σ ∆m threshold due their relatively large photometric uncertainties.

Photometric Redshift Estimation and SED Fits
To gain a better understanding of the transients and variable candidates in the field, we estimated photometric redshifts using EAZY (Brammer et al. 2008).This specifically allows us to understand the distances of the variable candidates, and to explore various properties of them, such as their masses, ages, dust extinction, and specific star formation rates (sSFRs).
Since that would be impossible with just the three HST filters, we also used PEARLS JWST/NIRCam images (see § 2).Consequently, we estimate redshifts only for the subset of galaxies that fall within the overlapping coverage of both HST and JWST.We opted to not use the HST F275W filter for photometric redshift estimations, because most objects remain undetected at sufficient significance.We thus employed a total of 10 filters for estimating photometric redshifts: HST/ACS F435W and F606W, and JWST/NIRCam F090W, F115W, F150W, F200W, F277W, F356W, F410M, and F444W.To ensure the HST and JWST images are on the same pixel scale, we reprojected each HST image onto the PEARLS F444W image pixel grid using reproject.In each filter, we set all pixels without coverage in all 10 filters to NaN values, facilitating flagging of objects at image edges.
We created photometric catalogs for input into EAZY using SourceExtractor.Recognizing that some galaxies exhibit variability that may bias a redshift estimation, we ran SourceExtractor on HST drizzled images of each individual visit instead of the full HST mosaic, and use the measurement with the fainter flux.We assume the fainter measurement represents the photometry if no variability were present -galaxies may temporarily brighten, but they cannot get dimmer on human time scales.For the JWST data, we used the full mosaics, since only a single epoch yet exists at any given location within the field.The NIRCam F444W image served as the SourceExtractor detection image in all cases.For HST images, we set MAG_ZEROPOINT to 26.49 mag for F606W and to 25.65 mag for F435W.For our JWST images, MAG_ZEROPOINT was set to 28.0865 mag for all filters, as appropriate for their 0. ′′ 030 platescale, to convert from MJy sr −1 (see Windhorst et al. 2023).Other relevant parameters are the same for HST and JWST, and are as listed in the right-most column of Table 1.
EAZY(v1.0) was run on the resulting 10-band photometry.Following the approach outlined in Figure 3, we scale SourceExtractor magnitude uncertainties by a factor of 3 upon input into EAZY.This factor of 3 is larger Color composites of 12 representative variable candidates, with F606W, F435W, and F275W shown in red, green, and blue hues, respectively.Each row of 100×100 pixel (3 ′′ ×3 ′′ ) cutouts correspond to a different sample: 1) F606W-only 3σ ∆m , 2) F435W-only 3σ ∆m , 3) F435W and F606W 2σ ∆m , and 4) F435W and F606W 3σ ∆m variability.A cyan or black circle at the center of each cutout represents the 0. ′′ 24 aperture in which variability was detected.Above each cutout, we note the visit identifier and decimal year of observation.At the top of each cutout, we include the measured m0.24 magnitude for that visit.In the left panel of each pair of cutouts, at the bottom, we list the variable candidate identification number and whether it is variable in the F606W filter, the F435W filter, or both.Most of the variable candidates will be AGN, although some will be faint or obscured SNe.
than the factor of 1.15 used for scaling the 0. ′′ 24 aperture magnitudes in §3.2, primarily due to the larger aperture sizes used here.We used the eazy_v1.0templates on a redshift grid spanning 0.01 ≤ z ≤ 15, allowing only single templates to be fit.We adopted z_ml as the redshift estimate and did not include any priors.
We used the EAZY redshifts as input to a custom IDL spectral energy distribution (SED) fitting code that was also recently used by Meinke et al. (2021Meinke et al. ( , 2023)).This code starts with Bruzual & Charlot (2003;BC03) SED models assuming a Salpeter IMF and an exponentially decaying star-formation history specified by the decay scale, τ , which ranges from 0.01 to 100 Gyr in 16 logarithmically spaced steps.Extinction by dust is specified on a grid of 0.0 ≤ A V ≤ 4.0 mag in steps of 0.2 mag assuming a Calzetti et al. (2000) extinction law.For a given redshift, the age, T , is not allowed to exceed the age of the Universe (assuming a Planck Collaboration et al. (2020) cosmology).The best-fit model is chosen by minimizing the χ 2 and scaling for the stellar mass, M .The (unweighted) present-day value of the star formation rate, Ψ(T ), is computed using the best-fit stellar mass, The y-axis labels on the left-hand side of each plot represents the number of galaxies in the general population, and those on the right-hand side (in blue) the number of galaxies with variability.Above each histogram panel, we plot the fraction of galaxies with variability (100% × Nvar/N general ) for each bin.The uncertainties assume Poisson statistics with σN ∝ √ N .This comparison only includes the subset of galaxies already flagged as reliable in § 3.2.1, with both HST and JWST imaging (where photometric redshifts and SEDs can be fit), with z ph < 5.5, χ 2 < 10, and where the 2σ confidence interval in p(z) spans <1 redshift unit.
With 10-band, 0.4-5 µm photometry, we required χ 2 < 10 and z < 5.5.Since all objects of interest are detected in F606W, they must be at z < 5.5.However, these cuts alone will not eliminate all bad fits, especially if the photometric uncertainties in some crucial filters are large.We added an additional condition that the limits on the 2σ confidence interval in redshift probability distribution, p(z), as computed by EAZY (namely l95 and u95) must not differ by more than 1 redshift unit.These conditions result in a sample with reliable photometric redshift estimates that includes 22 galaxies exhibiting variability (28% of the original sample of 79 variables that fall within the HST and JWST footprints) and 3296 normal galaxies (31% of the 10,664 galaxies that fall within the HST and JWST footprints).We only considered galaxies that were already marked as reliable in § 3.2.1 (the 26,468 sources in Table 4).
Due to the necessarily limited grid of SED model parameters, any galaxy with an inferred sSFR < 10 −5 yr −1 is categorized as "quiescent" and its sSFR is set to 10 −5 yr −1 .Our templates also cannot distinguish stellar population ages younger than 10 Myr, so we set all apparent age solutions <0.01 Gyr to 0.01 Gyr.Given that age is inversely correlated with the inferred sSFR, where log 10 (Age) ∼ −2 corresponds to log 10 (sSFR) ∼ 2, we also set all log 10 (sSFR)> 2 to 2. This adjustment is solely made to assess collective trends in galaxies with variability.
In Figure 8, we compare photometric redshift (z ph ), stellar mass, stellar population age, specific star formation rate (sSFR), and dust extinction (A V ) of galaxies that exhibit variability to typical galaxies in the general population.Overall, we sample a broad range of redshifts, masses, ages, star formation histories, and dust extinction.The most prominent peaks in the fraction of galaxies with variability with respect to the general galaxy population (expressed in %) occur at a redshift of z ph ≃ 2.5, a stellar mass of ∼10 8 and ∼10 10 M ⊙ , an age of ∼10 −0.5 Gyr, a sSFR of ∼10 −5 Gyr −1 , and an attenuation A V ≃ 1.0 mag.The median redshift of our variability sample is z ph = 1.3, but we sample variability up to z ph = 2.9 .Similarly, our galaxies with variability sample masses M = (0.02-20)×10 9 M ⊙ , ages T = (0.01-2)Gyr, sSFR = (10 −5 -10 2 ) Gyr −1 , and A V = (0.0-2.4) mag.
Three transients fell within an area that also has PEARLS 8-filter JWST/NIRCam coverage, yet all three had EAZY photometric redshifts with a χ 2 > 10.This could be due to a significant non-thermal component from an AGN, or due to the presence of strong emission lines that could mimic a continuum break.As transient T4 has a secure MMT/Binospec spectroscopic redshift of z = 0.615, we used the Code Investigating GALaxy Emission (CIGALE) (Boquien et al. 2019) to better fit the SED of its host galaxy.CIGALE has the benefit of incorporating emission lines into its models, while the previously mentioned SED fitting algorithm does not and resulted in poorer fits for this transient host.We used a delayed burst-quench star-formation history model (Ciesla et al. 2017), since the standard exponentially declining model can struggle to reproduce higher SFRs, and its morphology suggests that the transient host likely underwent a recent star-formation episode.The e-folding times tested were 1, 3, 5, and 8 Gyr, the main stellar population spanned 1-8 Gyr in increments of 1 Gyr, the burst/quench age ranged from 0.1, 0.5, and 1 Gyr, and the possible SFR ratios before and after these burst ages were 0.1, 0.3, and 0.6.Our stellar component used the BC03 library using the Salpeter IMF, testing metallicities of 0.004, 0.02, and 0.05 Z ⊙ with young stellar ages being either 10 or 100 Myr.We included a nebular component to account for emission lines using all the default parameters, where the nebular component is selfconsistent with the stellar population in the sense that its strength is determined by the Lyman continuum of the total SED.To account for dust attenuation, we used Calzetti et al. (2000) extinction with possible modified power-law slopes of −0.4,0.2, and 0. For dust emission, we used the Dale et al. (2014) model with α slopes of 1.0, 2.0, and 3.0 .
The best fit SED for the T4 host galaxy is shown in Figure 9 and corresponds to a stellar population age of 3.6 Gyr, A V = 0.5 mag of attenuation by dust, and a stellar mass of 1.2×10 10 M ⊙ .Its sSFR is 0.5 M ⊙ Gyr −1 .The reported best-fit value of 0.05 of the reduced χ 2 indicates that the (conservative) photometric uncertainties upon input to CIGALE were overestimating the actual ones somewhat.

DISCUSSION
An analysis of the detailed properties of individual transients and variable candidates is beyond the scope of this paper.We can, however, make some general statements about the properties of these source populations, using their observed optical light distributions (morphology) and EAZY SED fits.From the point-like morphology and off-center locations at small projected distances to likely host galaxies, we infer that the majority of the discovered transients are SNe.The majority of the variable candidates, on the other hand, are found either in or very near the center of galaxies or are isolated without a discernible host, and are surmised to be AGN.Most of these AGN are normal SMBHs in galaxies, with stochastically fluctuating accretion at lower luminosities than seen in quasars.Nonetheless, it is important to consider that other phenomena, such as SNe from central starbursts, stochastic microlensing events, or tidal disruption events, could also contribute to the observed central variability.SNe that are heavily obscured and SNe that are caught at very early or at late times can contribute to variability detected in local peaks in surface brightness well away from host galaxy centers.

Classification of Transients
Based on their positions relative to likely host galaxies and the PSF-like appearance in the difference images, we propose that nine of the transients (T1, T2, T3, T4, T5, T6, T8, T9, and T12) are SNe.It is reasonable to expect that ∼40-50% of detected SNe are Type Ia, as suggested by previous studies (Dahlen et al. 2012;Graur et al. 2014;Rodney et al. 2014;Cappellaro et al. 2015).If we assume that we have a minimum of nine SNe, we can anticipate a minimum of three Type Ia SNe with the remainder being CC SNe.
The TREASUREHUNT observing strategy provided piecemeal UV-Visible imaging rather than time-domain monitoring.As a result, these SN candidates will have long faded and can no longer be followed up for spectroscopic identification.Nonetheless, now that this initial imaging exists, future monitoring observations could efficiently expand the present sample and would allow such spectroscopic follow-up.The sum of the areas of overlap used in this work samples ∼88 arcmin 2 , which is about half of the total HST coverage of this field (∼194 arcmin 2 ).Therefore, we can expect about 2× as many transients if this field were observed again.Full overlapping coverage of this field could yield at least six new Type Ia supernovae, and coupled with spectra and rapid follow-up, can provide essential constraints on fundamental cosmological parameters like the Hubble constant, the mass density of the universe, the cosmological constant, the deceleration parameter, and the age of the Universe (e.g., Riess et al. 1998Riess et al. , 2023)).In addition to Type Ia SNe, comprehensive HST coverage of this field could yield at least six CC SNe observations.Particularly when combined with ancillary data, this could contribute significantly to resolving the observed discrepancies in CC SN rates compared to theoretical models (Cappellaro 2014).With a sufficiently large sample, the SN rate as a function of redshift could also be constrained.Furthermore, the increased volume of observations enhances the likelihood of capturing rare or unexpected events, potentially leading to groundbreaking discoveries in the field.
At least three of the transients are unlikely to be SNe.They either were detected in both epochs or they are isolated with no apparent host galaxy.When considering only the optical emission, transients T7 and T10 exhibit substantial changes in brightness between the two epochs, of ∆m = 0.83 ± 0.06 mag and 0.82 ± 0.09 mag, respectively.That would be equivalent to the presence of an additional source with m AB ∼ 24.82 and 25.52 mag.The observed centroid positions in both epochs match to within 0. ′′ 001, such that the additional signal is indistinguishable from being due to a nuclear source that changed brightness.We therefore suggest that T7 and T10 may be quasars.The fact that they both were detected in soft X-rays at flux levels >2×10 −15 erg s −1 cm −2 lends additional evidence to this interpretation.
T11 is a unique transient with no visible host, yet was detected in hard X-rays by NuSTAR.If that X-ray emission is indeed associated with T11 then it could be a faint (m AB ≳ 29 mag), high-redshift quasar that briefly flared up to be detectable in F606W.If the X-ray emission is unrelated, then we speculate that T11 may represent a SN located within an undetected host galaxy.That host could be a faint dwarf galaxy at z ≲ 6, or a more massive host at z ≳ 6 to explain the non-detection in F606W in the second visit.Unfortunately, T11 falls outside the area with JWST coverage, so we cannot distinguish between these two scenarios at this time.Chakrabarti et al. (2018) have demonstrated that the rates of SNe in dwarf galaxies can be three times higher (per unit volume) than in typical spiral galaxies.Given the limited understanding of star formation in dwarf galaxies, a substantial number of SN candidates like T11 could offer a unique opportunity to identify otherwise undetected dwarf galaxies and gain more insights into their star formation rates.
Overall, determining the nature of each transient remains challenging without spectroscopic data secured close in time.Nonetheless, these detected events showcase the considerable potential and effectiveness of transient science within the JWST NEP TDF.

Properties of Variable Candidates
Our methods were developed to isolate variable candidates that are coincident with local maxima in surface brightness that can be isolated and photometered using SourceExtractor.This naturally optimizes detection of variable AGN, but faint SNe are also detected when located offset from the center of the host galaxy.In addition, some offset variability may also be due to accretion disks of SMBHs that have not yet settled in the center of the host galaxy (e.g., Reines et al. 2020).Examples of offset variables are shown in Figure 7 (objects V001 and V169).
In general, the number of variable sources inferred here is expected to be fewer than the true number in the field, as some sources will vary on timescales that were not probed in this work.Constructing a sample of genuine variable sources in the field requires follow-up observations at various time sampling intervals, possibly extending over decades.
Figure 8 shows that galaxies exhibiting variability span a broad range of stellar mass, stellar population age, star formation history, and dust extinction.This underscores that rest-frame UV-Visible variability can be a powerful tactic to identify accreting AGN in a wide cross section of the galaxy population.Especially notable is that AGN were detected through their variability even in dusty galaxies with A V ≳ 2 mag.
The top portion of each panel in Figure 8, of the fraction of galaxies with variability within each histogram bin for the quantity plotted, shows peaks in most of the panels, although the formal significance of these peaks (assuming Poisson counting statistics) is generally ≲2σ due to small number statistics in individual bins.
First, in Figure 8a, we find a relatively high fraction of galaxies with variability around z ∼ 2.5, while at both lower and higher redshifts variability appears to be more rare.This is in contrast with both Villforth et al. (2012), who find optical-near-infrared variable AGN predominantly toward lower redshifts, and Sarajedini et al. (2011), who find an increase in optically varying AGN with increasing redshift.Zhong et al. (2022) find a similar number of optically variable AGN to peak at z ∼ 1.We suspect that the number of variable sources as a function of redshift is influenced by cosmic variance of a relatively rare source population, and thus differs depending on the field observed.The fraction of variable sources may be flat over the 0 ≲ z ≲ 3 redshift range when averaged over an area of sufficient size.
The observed stellar mass distribution (Figure 8b) of galaxies exhibiting variability shows an apparent excess at both ∼10 8 M ⊙ and ∼10 10 M ⊙ .We do not detect variability in the most massive galaxies, either due to small number statistics (most likely) or because their SMBHs are less likely to accrete at a sufficient rate to be detectable or have a longer period of relative quiescence between major accretion events.We also do not detect variability in galaxies with M ≲ 10 7 M ⊙ , likely due to a combination of both small number statistics and larger photometric uncertainties in these faint systems, such that our 3σ ∆m threshold is less likely to be met for a given amplitude of variability.If the relative distribution is taken at face values, then it could also hint at the presence of two separate populations of sources with variability: one associated with lower mass ∼10 8 M ⊙ hosts, and another associated with higher mass hosts ∼10 10 M ⊙ .Such could be the case if variability associated with CC SNe favors lower mass host galaxies, and variability associated with SN Ia and AGN favors higher mass hosts.
In Figure 8c and d we see that a significant fraction of the galaxies with variability are best fit with stellar population models characterized by active star formation (young ages ≲30 Myr and high sSFR ≳30 Gyr −1 ).Conversely, there is also a sizable fraction that is quiescent (sSFR ≲ 0.01 Gyr −1 ) with stellar populations older than 1 Gyr.A possible excess appears for intermediate population ages (60-600 Myr).
Last, Figure 8e shows that while the galaxies exhibiting variability largely track the distribution of extinction seen in the general population, there may be an excess in the variability sample for A V ∼ 1 mag.

Constraints on SMBH Mass, Accretion and Radiation Lifetimes
The sparse and random time sampling of this field does not allow us to trace the evolution of faint AGN, but may still teach us various properties of them.AGN brightness varies on timescales corresponding to the light crossing time in the accretion disks and gas clouds surrounding SMBHs.For example, for a SMBH with a mass M ∼ 10 8 M ⊙ , the timescale for variability cannot be shorter than ∼1 week in the rest-frame (e.g., Xie et al. 2005).For the majority of sources where the variability is due to AGN activity, that variability originates close to the central SMBH.If the timescale we sample is the minimum timescale of variation, we can provide rough constraints on the SMBH masses following Xie et al. (2005).According to Schwarzschild Black Hole Theory (SBHT), the mass of a black hole (M • in M ⊙ ) is equal to the minimum timescale of variation (∆t min in yr) as M • = 4.29×10 11 ∆t min .Time intervals between observations in various areas of overlap in the TREA-SUREHUNT data range from 1 day (0.0027 yr) to 4.78 yr.At the median redshift of our variable candidate sample, z ph = 1.3, the range in rest-frame timescales becomes 0.001-2.08yr, giving SMBH masses of 5.1×10 8 ≲ M • ≲ 8.9×10 11 M ⊙ .SMBH masses ≳10 11 M ⊙ are unreasonable: objects within the region of overlap with a 4.78 yr interval between observations likely vary on much shorter timescales and we did not sample the minimum timescale of variation.Even a more typical time interval of ∼1 yr between epochs will yield a ∼10 11 M ⊙ SMBH mass.Taking the full range in redshifts sampled into account, we will therefore only claim that in select areas of overlap we can probe SMBHs with masses of M • ≳ (3-9)×10 8 M ⊙ , but with no firm upper limit.
We can also frame our finding that a maximum of 0.42% of the general field galaxy population shows significant variability in the context of the accretion and radiation lifetimes of AGN.If we assume that: 1) all galaxies seen in this field with HST have a central SMBH (e.g., Kormendy & Ho 2013), 2) the central SMBH is always visible when accreting, and 3) that a currently accreting AGN will be variable, then we would deduce that the AGN in our field are actively accreting for 0.42% of the time.Assumptions 2) and 3) are simplified, as some AGN will be hidden at rest-frame UV-Visible wavelengths by surrounding dust.This fraction could be as large as 2/3 of all AGN, if we use the ∼30% average Lyman continuum escape fraction of weak AGN (Smith et al. 2018(Smith et al. , 2020(Smith et al. , 2024) ) as proxy for the fraction of AGN with direct unobscured sight lines to the observer.Figure 8c implies an average SED age of ∼10 8 years (∼100 Myr) for the parent population, although the spread in age is wide (0.01-10 Gyr).A 0.42% fraction of visibly variable AGN would correspond to an implied average AGN activity lifetime of 0.42%×10 8 ∼ 4×10 5 yr (and a range of 4×10 4 -4×10 7 yr).For comparison, Rawes et al. (2015) estimated the optical synchrotron electron life times of AGN with visible jets observed with HST and Chandra at ∼10 4 yr.However, these authors state that their estimated synchrotron lifetimes may be too short by at least a factor of two, and may be longer if synchrotron electrons are re-accelerated in the ambient magnetic field.On the other end of the activity scale, Jakobsen et al. (2003) show that at least some quasars can remain active more or less continuously for ≳10 Myr.With these assumptions and significant uncertainties, a visible AGN accretion time of ∼10 4 to 10 7 yr could indeed result in a 0.42% variability fraction in a galaxy population characterized by an average SED age of ∼100 Myr (and a range of 0.01-10 Gyr) whose stars and gas feed that central engine.
Future work will need to secure and analyze larger samples of optically variable sources with deep X-ray imaging to better constrain UV-Visible synchrotron lifetimes, and wider and deeper HST+JWST images with spectroscopic (NIRISS or NIRSpec) redshifts to improve redshifts and the characterization of their stellar population.

Caveats and Estimated Population Size
It is important to note that with HST data alone, distinguishing between variability caused by variable AGN and faint SNe is challenging.Variability at locations offset from the cores of galaxies could potentially indicate SN events or SMBHs that have not yet reached the galaxy center.In all four catalogs, we identify a total of 25 cases showing variability at locations other than the core, as identified through visual inspection.Identifying AGN will require complementing this work with ancillary observations of the field to determine whether galaxies exhibiting variability also emit mid-IR emission.Prolonged monitoring of variable objects within this field will contribute to a deeper understanding of SMBHs within galaxies, including their evolution over time and their influence on galactic physics.
We also note that the variable candidates presented here are a strict lower bound to the number present in the field: many genuine variable sources, even if at z ≲ 6 and detectable in F606W, were certainly missed because they did not vary sufficiently between specific individual visits to meet our 3σ ∆m threshold, or if the timescale of their variability exceeded the duration of the TREASUREHUNT program.In the present study, the photometric uncertainties in each epoch of observation contributed equally to that 3σ ∆m threshold.Any future additional observations to similar depths would reduce the photometric uncertainty in the reference magnitudes, thus lowering the variability detection threshold.The magnitude and level of variability for all sources in this work is collected in Table 5, so that they can be compared in future studies within the JWST NEP TDF.
Although significant (i.e., >3σ ∆m ) direct detections of variability to such exceeding faint limits through systematic monitoring campaigns with HST, JWST, or Roman would require a very large community investment, it is not outside the realm of what is possible.Even at slightly shallower depths, and employing similar extrapolations, variability could account for a large fraction of all accreting SMBHs.

CONCLUSION
This study highlights the JWST NEP TDF as an exceptional site for investigating time-variable astronomical phenomena, including SNe and (weak) variable AGN.We identified 12 transients and ∼100 variable sources in ACS/WFC F606W and F435W images from the HST TREASUREHUNT program.We provide positional information, dates of observation, intervening time interval, and magnitudes for each of the 12 transients, which range in brightness from ∼24.8 to ∼27.5 mag.Their areal density is ∼0.07 transients per arcmin 2 (∼245 per deg 2 ) per epoch.We argue that the vast majority of these transients are SNe.Three transients (T7, T10, and T11) were detected in X-rays, of which two (T7 and T10) appear isolated and point-like and are likely newly identified quasars.Transient T11 has no visible host galaxy, and its nature remains uncertain.We suspect it to be a faint quasar if the X-ray detection is indeed associated, or otherwise either part of a very faint dwarf galaxy at z ≲ 6, or a transient in a more massive host at z ≳ 6.One transient (T4) has a spectroscopic redshift z = 0.615 from MMT/Binospec.
Our variability search revealed that 0.42% of the general z ≲ 6 field galaxy population exhibits variability at ≳3σ ∆m significance to depths of 29.5 mag in F606W and 28.6 mag in F435W, with an areal density of ∼1.25 variables per arcmin 2 (∼4500 per deg 2 ).We carefully identified variable candidates, using the measured distribution of magnitude differences for objects observed more than once to calibrate our photometric uncertainties, where the scaled photometric uncertainty is σ ∆m , as defined in Equation 1. Sources are flagged as variable if they varied by more than 3σ ∆m in F435W or F606W, or varied by more than 2σ ∆m in both filters in the same sense.This revealed 190 unique variable candidates, of which we estimated ∼80 to be false positives (using Gaussian statistics).Most of the variable candidates are coincident with the cores of galaxies, indicating potential AGN variability, while a smaller fraction appears associated with faint SNe.We also estimate photometric redshifts for our sample, from which we estimate masses, ages, dust extinction, and sSFRs.We briefly explore ways to estimate SMBH mass and AGN timescales using the time interval between overlapping visits and the observed rates of variability.
In conclusion, this work firmly establishes the JWST NEP TDF as a pivotal field for time-domain science with an initial harvest of transients and galaxies exhibiting variability.We emphasize the ability to identify AGN through their variability in the JWST NEP TDF.Future follow-up observations with HST or JWST could greatly increase the numbers of directly detected variables and to lower variability amplitudes.Figure 10.We analyze the potential effect of cosmic rays on the measured flux difference of variable objects in an 8×8 pixel (0. ′′ 24×0.′′ 24) area centered on each object.A pixel is considered at risk of cosmic ray contamination if it encounters a cosmic ray in four or more contributing exposures.Plotted is a random selection of 10,000 galaxies (represented by black points), with median values indicated by green boxes.Left: The absolute value of the difference in measured F606W flux versus the total number of pixels (in Visit 1 + Visit 2) possibly impacted by cosmic rays.The red shaded region shows objects that are removed for containing too many cosmic rays.The y-axis is logarithmically scaled except between -0.1 and 0.1, which is scaled linearly.The x-axis is linearly scaled up to 20, where it is then logarithmically scaled to 50.Right: The absolute value of the difference in measured F606W flux versus the difference in the number of pixels potentially affected by cosmic rays.
One object with very high measured differences in flux also have a large number of pixels with cosmic rays.We opt to reject and remove two objects with more than 25 pixels (∼25% of the pixels in the aperture) affected by cosmic rays.There were no variable candidates that met this criteria.
We also confirm the distribution of ∆m 0.24 is Gaussian for sources potentially affected by cosmic rays, shown in Figure 11.More importantly, since a majority (∼ 80%) of sources in our dataset are potentially affected by cosmic rays, our empirical uncertainties in Section 3.2.2already account for cosmic ray affects.This is specifically represented by the fact that the Gaussian fit to the CR-affected sources (red line) is within 0.01 mag to the Gaussian fit for the full data (black line).The latter (black line) was used to calibrate our empirical uncertainties in Section 3.2.2.We note that, by eye, there appears to be a slight divergence from the Gaussian fit at σ ∆m < −0.5 and σ ∆m > 0.5.However, only the outer-most points seem to contribute to this trend, corresponding to less than 5 sources (or less than 0.2% of sources).
In conclusion, this evidence suggests that cosmic rays rarely appear in the final drizzled images.The spurious transients (in Section 3.1) were likely outliers, and any other spurious cosmic rays to this degree would have been identified during the transient search.In addition, if cosmic rays affect the brightness of sources regularly, our empirical estimation of uncertainties in Section 3.2.2already accounts for this.Nonetheless, we do not completely eliminate the possibility that some (albeit very few) of our variable candidates may be affected by cosmic rays that cause artificial variability.In other words, although we assume a majority of our sources are not affected by cosmic rays, follow-up observations remain necessary to confirm which sources are truly variable, as well as reveal other sources that were missed in this work.Histogram of the distribution of magnitude differences measured between two epochs (∆m0.24) for 4,059 galaxies between 28.0 < m AB < 28.5.We compare the sample of CR-affected galaxies (red) to the full sample (black).The points show the fraction of CR-affected galaxies (N ) for each ∆m0.24 bin, where the error bar is √ N .The solid lines are Gaussian fits to the data.The black solid line represents the Gaussian fit to the full sample.Our empirical uncertainties estimated in Figure 3 already account for cosmic ray effects, as shown by the black line closely following the red line.

Figure 1 .
Figure1.Full extent of HST ACS/WFC coverage of the JWST NEP TDF.Regions of overlap between visits in the mosaic are highlighted in progressively darker shades of grey for areas of 2-, 3-, and 4-epoch overlap.This paper focuses on these regions of overlap, with a total area of 87.94 arcmin 2 , to investigate transients and variable sources.Red circles mark the positions of the 12 identified transients, while blue squares and diamonds mark the 190 unique variable candidates discovered in this study.We expect ∼80 of these sources to be false positives (see § 3.2).Symbol shapes and colors identify the samples from which individual variables were selected, as indicated by the legend at the top of the figure.For comparison, the JWST/NIRCam coverage of this field is outlined in orange.We estimate photometric redshifts ( § 4.3) for objects with both HST and JWST photometry.One visit (jeex02; at bottom right) with a bad astrometric solution was ignored in this work.

Figure 2 .
Figure 2. Star-galaxy separation, following Windhorst et al. (2011), using the SourceExtractor MAG_AUTO magnitudes in F606W versus full width at half maximum (FWHM).Objects within grey regions, with FWHM < 0. ′′ 15 or m AB < 17 − 4 log 10 (FWHM), are flagged as stars and are not considered in this study.We also show the radius (0. ′′ 12) of the aperture used to identify variability as a vertical blue dashed line.

Figure 3 .
Figure3.Calibration of photometric uncertainties to identify variable sources.We plot the difference in measured F606W magnitude, ∆m0.24, versus the average measured F606W magnitude, m0.24, for objects observed more than once in regions of overlapping HST coverage, where ∆m0.24 is defined in the sense Visit 2 − Visit 1.Each small grey point is an individual object, as detected and measured with SourceExtractor.The large brown circles depict SourceExtractor photometric errors, scaled by a factor 1.15 to match the observed distribution of magnitude differences (shown as small black circles) measured in 0.5 mag wide bins for 25.5 < m0.24 ≤ 28.0 mag such that ∼68.3% of all data points fall between them.The dashed black curves represent the inner ∼68.3% range of the distribution at each magnitude and are polynomial fits to the black points.The brown curves therefore represent median 1, 2 and 3σ ranges as fit to the brown circles.Blue squares identify objects that vary at ≥3σ ∆m in F606W, where σ ∆m is defined in Equation 1.The histogram on the right-hand side illustrates the distribution of magnitude differences for all objects detected in the field.In grey, we list mean, median, standard deviation, and median absolute deviation of the y-axis distribution.

Figure 4 .
Figure 4. Histogram of the distribution of magnitude differences measured between two epochs (∆m0.24) for all galaxies (4,059 total) between 28.0 < m AB < 28.5 mag.The black circles show the number of galaxies (N ) for each ∆m0.24 bin, where the error bar is √ N .The solid black curve is a Gaussian fit to the data.The vertical blue lines represent the 3σ thresholds based on the Gaussian fit.This proves that the distribution of ∆m0.24 follows a Gaussian distribution, such that the false detections can be assumed to follow Gaussian statistics (see Section 3.2.3).

Figure 5 .
Figure5.Histograms of the number of variable candidates (light grey) and estimated numbers of false positives (red) and genuine variable sources (blue) as a function of the significance (in units of σ ∆m ) of the variability detected.The number of genuine variables and false positives add up to the total number of candidates in each bin.The red labels above each bin show the percentage of all objects in the field that are expected to be false positives for that σ ∆m -bin.The top panel is for sources that varied by more than 3σ ∆m in F606W, and the bottom panel for sources that varied by more than 2σ ∆m in both F606W and F435W, and did so in the same sense.Hence, assuming Gaussian statistics, the upper panel reflects a two-sided test, and the lower panel a one-sided test.In total, we estimate ∼80 false positives in the sample of sources that varied by more than 3σ ∆m in F606W, the sample of sources that varied by more than 3σ ∆m in F435W, and the sample of sources that varied by more than 2σ ∆m in both F606W and F435W.

Figure 6 .
Figure 6.F606W images of each of the 12 transients identified within the JWST NEP TDF.For each transient we show a 3-panel composite of 100×100 pixel (3 ′′ ×3 ′′ ) cutouts from the image in which the transient is present (left), the image in which the transient is absent (middle), and the difference image (right).Labels give the transient ID, the visit IDs and dates of observation, the time interval between epochs, and the apparent F606W magnitude of the transient.Transients T7, T10 and T11 were also detected in X-rays (see Table 3), with T7 and T10 likely quasars.Grey-scale bars to the right of each 3-panel set are in units of e − /s.
Table columns list (1) the transient ID, (2) the X-ray observatory that also detected the transient source, (3) the measured X-ray flux, and (4) the energy band in which that flux was measured.The XMM-Newton 2-10 keV upper limits correspond to the 90% confidence level.

Figure 7 .
Figure 7.Color composites of 12 representative variable candidates, with F606W, F435W, and F275W shown in red, green, and blue hues, respectively.Each row of 100×100 pixel (3 ′′ ×3 ′′ ) cutouts correspond to a different sample: 1) F606W-only 3σ ∆m , 2) F435W-only 3σ ∆m , 3) F435W and F606W 2σ ∆m , and 4) F435W and F606W 3σ ∆m variability.A cyan or black circle at the center of each cutout represents the 0. ′′ 24 aperture in which variability was detected.Above each cutout, we note the visit identifier and decimal year of observation.At the top of each cutout, we include the measured m0.24 magnitude for that visit.In the left panel of each pair of cutouts, at the bottom, we list the variable candidate identification number and whether it is variable in the F606W filter, the F435W filter, or both.Most of the variable candidates will be AGN, although some will be faint or obscured SNe.

Figure 8 .
Figure 8.Comparison of galaxy properties in the general galaxy population and in galaxies that exhibit variability.Shown are distributions of (a) photometric redshifts (z ph ), (b) stellar mass, (c) stellar population age, (d) specific star formation rate (sSFR), and (e) extinction by dust (AV ), for the general population (solid grey) and galaxies with variability (blue hatched).The y-axis labels on the left-hand side of each plot represents the number of galaxies in the general population, and those on the right-hand side (in blue) the number of galaxies with variability.Above each histogram panel, we plot the fraction of galaxies with variability (100% × Nvar/N general ) for each bin.The uncertainties assume Poisson statistics with σN ∝ √ N .This comparison only includes the subset of galaxies already flagged as reliable in § 3.2.1, with both HST and JWST imaging (where photometric redshifts and SEDs can be fit), with z ph < 5.5, χ 2 < 10, and where the 2σ confidence interval in p(z) spans <1 redshift unit.

Figure 9 .
Figure 9. SED fit of transient T4 with CIGALE, adopting the MMT/Binospec spectroscopic redshift of 0.615.The red points with vertical error bars show the measured fluxes of the transient host and their associated uncertainties.The horizontal bars on each data point represent the width of the bandpass.The light grey spectrum represents the best CIGALE fit, with the reduced χ 2 value shown.
Figure11.Histogram of the distribution of magnitude differences measured between two epochs (∆m0.24) for 4,059 galaxies between 28.0 < m AB < 28.5.We compare the sample of CR-affected galaxies (red) to the full sample (black).The points show the fraction of CR-affected galaxies (N ) for each ∆m0.24 bin, where the error bar is √ N .The solid lines are Gaussian fits to the data.The black solid line represents the Gaussian fit to the full sample.Our empirical uncertainties estimated in Figure3already account for cosmic ray effects, as shown by the black line closely following the red line.

Table 1 .
Relevant SourceExtractor parameters for variability identification and aperture photometry.

Table 2 .
Transients identified in HST TREASUREHUNT observations of the JWST NEP TDF Tablecolumns(1) through (8) list the transient ID, the celestial coordinates (identified by eye) of the transient in decimal degrees, the root names of the visit where the transient is present (Visit 1) and where it is absent (Visit 2), the corresponding dates of observation, and the time interval ∆t between Visits 1 and 2. The two final columns list the measured F435W and F606W magnitudes of the transient (m AB at t1 relative to t2, measured within a circular aperture centered on the transient with a radius of 8 pixels), where the magnitude uncertainties are given in parentheses.

Table 3 .
Transients identified in HST observations of the JWST NEP TDF with significant X-ray detections

Table 4 .
Number of variable candidates identified and inferred in F606W and F435W.