Multiwavelength Spectral Analysis and Neural Network Classification of Counterparts to 4FGL Unassociated Sources

The Fermi-LAT unassociated sources represent some of the most enigmatic gamma-ray sources in the sky. Observations with the Swift-XRT and -UVOT telescopes have identified hundreds of likely X-ray and UV/optical counterparts in the uncertainty ellipses of the unassociated sources. In this work we present spectral fitting results for 205 possible X-ray/UV/optical counterparts to 4FGL unassociated targets. Assuming that the unassociated sources contain mostly pulsars and blazars, we develop a neural network classifier approach that applies gamma-ray, X-ray, and UV/optical spectral parameters to yield descriptive classification of unassociated spectra into pulsars and blazars. From our primary sample of 174 Fermi sources with a single X-ray/UV/optical counterpart, we present 132 P_bzr>0.99 likely blazars and 14 P_bzr<0.01 likely pulsars, with 28 remaining ambiguous. These subsets of the unassociated sources suggest a systematic expansion to catalogs of gamma-ray pulsars and blazars. Compared to previous classification approaches our neural network classifier achieves significantly higher validation accuracy and returns more bifurcated P_bzr values, suggesting that multiwavelength analysis is a valuable tool for confident classification of Fermi unassociated sources.


INTRODUCTION
The Fermi Gamma-ray Space Telescope -Large Area Telescope (Fermi -LAT) 4FGL catalog includes 5064 sources, 3376 of which are associated with extragalactic blazars or nearby pulsars (Abdollahi et al. 2020). 352 other 4FGL sources include supernova remnants, X-ray binaries, starburst galaxies, and other objects. 1336 4FGL sources are "unassociated", lacking confident astrophysical explanations or source counterparts at longer wavelengths. Given that the bulk of the 4FGL sources are blazars or pulsars it is feasible that many of the unassociated sources are also blazars or pulsars that failed to be classified. The unassociated sources therefore tease systematic expansions to blazar and pulsar catalogs by extending to less obvious sources.
Identifying blazars and pulsars among the unassociated 4FGL sources is an important step towards confident population studies of both classes. The Fermi -LAT unassociated sources might include blazars that are lower luminosity, higher redshift, or viewed at a larger angle with respect to the jet than their more easily detected and identified cousins in the established Fermi blazar catalog (Ferrara et al. 2015). A more complete blazar catalog will serve numerous scientific goals, including the verification and analysis of the blazar sequence (e.g., Fossati et al. 1998;Ghisellini et al. 2017) as a theoretical unifying scheme for blazars.
Similarly, there may be undiscovered pulsars among the unassociated 4FGL sources, valuable additions to the significantly shorter list of Fermi pulsars. The Fermi pulsar list includes canonical and millisecond pulsars but is also notable for dramatically expanding research into 'Black Widow' pulsars (e.g. Wu et al. 2018;, and the unassociated sources may contain several pulsars of these various types. Finally, some unassociated objects may defy classification as blazars or pulsars even after multiwavelength analysis. Identifying and classifying the "low-hanging fruit" of pulsars and blazars among the unassociated sources is a first step in identifying and studying other astronomical objects. While gamma-ray observations from Fermi -LAT are the foundation for the 4FGL catalog, more confident classification of pulsars and blazars can be achieved by extending spectral analysis to lower energies. To this end, the Neil Gehrels Swift Observatory (aka Swift) (Gehrels et al. 2004) is conducting a continuing survey of the Fermi unassociated sources. With the XRT (Burrows et al. 2005) and UVOT (Roming et al. 2005) instruments onboard Swift, X-ray observations from 0.3 − 10.0 keV and UV-visual observations from 450 − 900 nm systematically cover the uncertainty ellipse of each Fermi -LAT unassociated source.
Given the variety of spectral features in blazars and pulsars, this multiwavelength capability is a powerful tool for studying 4FGL unassociated sources at lower energies and for supporting identification and classification efforts. For example Kaur et al. (2019) sorted 217 high-S/N unassociated 3FGL sources into blazars and pulsars using machine learning on gamma-and a single X-ray spectral parameter. The sample was limited to only those sources with a single possible X-ray counterpart in the error ellipse. Training an ML routine with known gamma-ray pulsars and blazars, the authors identified 173 likely blazars with P bzr > 90% (134 with P bzr > 99%) and 13 likely pulsars with P bzr < 10% (7 with P bzr < 1%). From their initial list, 31 sources from the 3FGL unassociated list defied categorization and were labeled 'ambiguous'. Recent work by this collaboration (Kerby et al. 2021) continued that analysis by including more detailed X-ray spectral fitting and systematically adding spectral parameters to machine learning. Adding UV-visual observations from the Swift-UVOT telescope is particularly useful, as pulsars are usually extremely dim in the UV-visual range Saz Parkinson et al. 2016) while blazars emit at all wavelengths and at low redshift can be observed across the electromagnetic spectrum (Ghisellini & Tavecchio 2008).
The basis of expanding analysis of the unassociated sources to lower energies is the assumption that a gamma-ray source is likely to also be emitting X-rays due to the ubiquity of X-ray synchrotron emission in energetic gamma-ray systems. Given the spatial resolution and X-ray sensitivity of our observations with Swift-XRT and the ubiquity of pulsars and blazars emitting both in gamma-and X-rays, if only a single X-ray source is present in a gamma-ray ellipse then there is likely a relationship between the two. Because gammaray emitters like pulsars and blazars normally emit Xrays via mechanisms like synchrotron radiation, an association between an unassociated gamma-ray source and a solitary X-ray sources within its uncertainty ellipse is theoretically sound. The automated analysis after Swift observations calculates the probability of a coincident but unrelated X-ray source in the Fermi 4FGL 95% confidence region based on exposure time and confidence region size, and for a typical 4ks exposure in these ∼5 arcmin semi-major axis regions, the probability is < 0.01. For the handful of exposures with longer than 4ks and/or larger than typical 4FGL confidence regions sizes, the probability of additional spurious X-ray source detections increases. Finally, there are a large number of 4FGL target fields for which no X-ray detection was found in the typical 4 ksec exposure, which is consistent with the estimated low chance probability for spurious X-ray source detection.
In previous work on the Fermi unassociated sources (e.g. Kaur et al. 2019;Kerby et al. 2021) analysis was restricted to Fermi unassociated sources with only one high-S/N X-ray source in their gamma-ray confidence ellipse. Still, it may be the case that Swift observations reveal multiple X-ray sources in the uncertainty ellipse of a single Fermi unassociated gamma-ray source, as a 4FGL gamma-ray uncertainty ellipses spans several arcminutes while the positional uncertainty of a XRT detections is a few arcseconds. In this case, other factors must be considered before determining a likely association between high-and low-energy photons. It is less likely that there will be degeneracy between UV-visual and X-ray sources, as both the XRT and UVOT instruments on Swift have relatively small PSFs.
Our approach to classification of 4FGL unassociated sources has several unique factors compared to other classification efforts and builds on our work on the 3FGL catalog. While Lefaucheur et al. (2017) conducted IR observations of unassociated sources to look for lowenergy counterparts, skipping over X-ray/UV/optical observations would complicate attempts to use spectral information to link gamma-ray sources to counterparts. The sky is much more prolifically filled with IR sources than X-ray sources, and a gamma-ray uncertainty ellipse with just one possible X-ray counterpart might have dozens in the IR bands. Searching in the radio band is similarly difficult, requiring assumptions to link gamma-ray and radio emission without interpreting the intermediate wavelengths, but the detection of characteristic radio pulsations is a direct route to locating pulsars and some previous works have used radio properties to predict pulsar membership after efficient searches (for example, Frail et al. 2018). Recently, Zhu et al. (2021) conducted ML on the 4FGL unassociated sources, but restricted their analysis to only gamma-ray properties. Analysis of X-ray observations specifically can capture the synchrotron peaks of non-thermal emitters like blazars, a valuable region for discriminating the spectra of pulsars and blazars.
In this work, we extend investigations of the unassociated sources by combining gamma-ray data with Swift-XRT and -UVOT observations for 4FGL unassociated targets observed thus far by Swift's systematic program. With this multiwavelength set of parameters, we build a neural network classification (NNC) routine using samples of known pulsars and blazars to classify the unassociated sources. In section 2, we describe the gamma-ray, X-ray, and UV/optical observations and sources used in this work, plus reasons for excluding certain entries from our base sample. Next, in section 3 we describe our production of spectral parameters and present tabulated results. In section 4 we discuss the classifications of the unassociated gamma-ray sources via a NNC method. In section 5 we present and discuss spectral fitting and classification results, summarize our findings, and posit next steps.

Fermi-LAT Unassociated and Swift-XRT Sources
The unassociated sources of the Fermi -LAT 4FGL catalog are comprised of gamma-ray sources of unknown astrophysical nature and no known counterpart (for a discussion of the entire 4FGL catalog, see Fermi-LAT-Collaboration 2019). After ten continuous years of observations, the 4FGL-DR2 catalog (Ballet et al. 2020) contains 5064 sources in total, of which over 1300 are unassociated. 1410 4FGL sources are targets for the Swift-XRT survey of unassociated sources (slightly higher than the total number of currently unassociated sources due to sources gaining or losing associations between DR1 and DR2). By observing the unassociated gamma-ray sources, the Swift program provides highresolution X-ray, UV, and visual observations to find lower-energy counterparts to the unassociated sources. This survey is a further development after previously analyzing unassociated 3FGL sources (detailed in Kerby et al. 2021) and has covered approximately 500 4FGL unassociated targets with > 4 ks of observations each. So far Swift has detected possible X-ray counterparts within the Fermi -LAT uncertainty ellipse of 208 4FGL unassociated sources, with some unassociated sources having multiple possible X-ray counterparts.
Combining archival and new Swift observations, an automated analysis process described in Falcone et al. (2011) produces lists of X-ray sources. For our purposes, a source is deemed 'notable' if it is contained within the 95% uncertainty ellipse of the 4FGL source and if the signal-to-noise ratio calculated with the Ximage function SOSTA is greater than 4. We take the produced list of 238 notable X-ray sources as containing the possible Xray counterparts to the 4FGL unassociated sources, and focus on these detections as the sample for further X-ray analysis.
The HEASARC query interface allows for downloads of Swift-XRT observations within 8' of all 238 notable Xray sources spread across the 205 4FGL targets. Some X-ray sources were matched with Swift-XRT observations only upon expanding the HEASARC search radius to 10', as they sit close to the edge of the field of view. Expanding the search radius of the HEASARC query does not impinge too closely towards the edges of the 23.6 field of view of the Swift-XRT telescope.
We immediately eliminate two 4FGL targets from consideration, along with all related possible counterparts. The two targets, J1836.8-2354 and J1649.2-4513, have very uneven Swift-XRT coverage near the 4FGL centroid position due to a bright, highly-observed X-ray source just outside the field of view. The nine notable X-ray sources from those two targets appear to be related to the extremely long (> 100 ks) observations of the nearby X-ray source that is outside the uncertainty ellipse. Several other 4FGL targets are eliminated if their only possible X-ray counterpart is coincident with a catalogued star that contributes a false X-ray signal via optical loading, the pile-up of optical photons in the Swift-XRT detector, or by a star of such brightness that it would radically dominate the UV/optical emission from an unassociated source, putting it outside of the range of known pulsars and blazars.
Of the 205 4FGL targets with possible X-ray counterparts examined in this work, 14 have more than one notable X-ray source within their uncertainty ellipse. For the purposes of spectral analysis and ML classification, we fit all X-ray sources while making no judgements about which is most likely the counterpart to the gamma-ray source. We tabulate and discuss the unitary sources as our primary sample separately from the more "confused" unassociated sources with more than one X-ray/UV/optical possible counterpart.

Swift-UVOT Source Analysis
For each X-ray source, we also used Swift-UVOT data within 8-10 of the XRT centroid to search for X-ray/UV/optical counterparts using the Swift-UVOT telescope. For each XRT/UVOT source, we generated a fully integrated UVOT image by performing UVOTIMSUM on the UVOT FITS file. We placed circular extraction regions at the coordinates of these UVOT detections; if there was no UVOT detection within the ∼ 5 positional uncertainty of the XRT detection, we centered the extraction region on the X-ray source's centroid position to determine an upper limit to the brightness. The radius of the UVOT extraction region was set to 5 , the PSF of the Swift-UVOT telescope for a point source. UVOTSOURCE was then used to gather counts from the source region, and the background was measured using a source-free circular background region of radius 20 elsewhere in the image for each position. The background-subtracted count rate was converted to flux and magnitude following the process described in Breeveld et al. (2011) to obtain a UVOT magnitude in one of the bands in Table 2. UVOT counterparts imaged in multiple UVOT bands were analyzed once separately for each band.
We applied a 3σ detection threshold for UVOT counterparts using UVOTDETECT. Though the vast majority of XRT detections had unique UVOT counterparts as well, fifteen regions had clustered UVOT sources that could not be disentangled within the constraints of the PSF of the telescopes on Swift. Training, validation, and research datasets with complete column coverage are vital to the machine learning process described below, so we exclude those confused sources from the unassociated sample.

X-ray Spectral Analysis
In total the HEASARC query for Swift-XRT observations collected over 1000 individual observations. Here, we only use observations in the photon counting (PC) mode of the Swift-XRT, enabling two-dimensional imaging across the XRT field-of-view with reasonable energy resolution. Summed exposure time for the X-ray sources varies from a few kiloseconds to many dozens of kiloseconds.
Each level 1 event file was processed and cleaned using xrtpipeline v.0.13.5 from the HEASOFT software 1 . Only events graded 0 through 12 were used in this analysis. Merging the events and exposure files with other 1 https://heasarc.gsfc.nasa.gov/docs/software.html observations of each particular X-ray source was conducted using xselect v.2.4g and ximage v.4.5.1, resulting in a single summed event list for each source plus a summed exposure map and ancillary response file using xrtmkarf. This merging precludes any time series or variability analysis of the X-ray sources, but few X-ray sources in our sample have sufficient observation time or photon counts to enable detailed temporal analysis.
For each possible X-ray counterpart, xselect produced spectra for source and background regions. The source region was circular with radius 20 arcseconds, and the background region was annular with inner and outer radii of 50 and 150 arcseconds respectively. Both regions were centered on the coordinates of the examined X-ray source. The background region size was chosen to be far outside the PSF of a point source for the Swift-XRT detector (half-power diameter of 18 ).
If the count rate in the center of the source region exceeded 0.5 counts per second (with PC mode having frames of 2.5 s each, this is approximately greater than one photon per frame), drawing a new annular source region with an inner radius depending on the count rate would avoid photon pile-up and saturation on the detector. The possible X-ray counterparts to 4FGL unassociated sources are faint enough that none caused such pile-up.
We used Xspec v.12.10.1f (Arnaud 1996) to fit each spectrum. The fitting model included three nested functions: tbabs, cflux, and powerlaw. cflux calculated the total unabsorbed flux between 0.3 and 10 keV and tbabs modeled line-of-sight hydrogen absorption using galactic values from the nH lookup function described in Wilms et al. (2000). The galactic line-of-sight extinction is fixed at the catalog value for each spectrum analyzed. powerlaw is a simple power law with photon index Γ X . Uncertainties on the fitted photon index and X-ray flux were jointly measured using the iterative steppar routine, and the spectral fitting results are included in an accompanying machine-readable database.
Fitting was executed using the C-statistic as the optimization metric. The C-statistic is a useful fitting statistic for spectra with few counts, particularly in cases for which there are not sufficient photons to bin counts for a χ 2 fit (Cash 1976). Compared to the χ 2 statistic, which assumes Gaussian behavior in bins, the C-statistic operates under Poisson statistics. In this way, it is much more applicable to fitting spectra with very few counts in each energy bin. The C-statistic is given by with t the exposure time, m i the predicted count rate in any particular bin, and S i is the observed counts in each bin. For an X-ray source with only a few dozen detected X-ray photons, binning the events for a χ 2 approach would result in only one or two bins, which is not useful for detailed spectral fitting. The Cash statistic does not require any such binning by assuming the more appropriate Poisson distribution of events. The model incorporating hydrogen-absorbed power law spectra returned unusually high or low Γ X for nine 4FGL X-ray counterparts compared to the expected 0 < Γ X < 4 for pulsars and blazars. Given that some of these X-ray sources also coincided with catalogued stars, we viewed these sources as dubious for pulsar/blazar classification but worthy of additional investigation with other approaches. We do not include these spectra in the ML classification. These sources are listed in Table 1.

V-magnitude Conversions
The analysis of Swift-UVOT data produces magnitudes in the six UVOT bands (vv, bb, uu, w1, m2, w2). However, the training sample of known pulsars and blazars uses Johnson V magnitudes, requiring a conversion from the UVOT bands to the Johnson system. To facilitate this conversion, we convert the UVOT magnitudes to fluxes, then use a power-law scaling relationship to predict the V-band flux, converting back to magnitude. The power-law scaling relation uses the central wavelengths of the UVOT bands (given in Table 2) and the V band (540 nm) plus an assumed UV/optical spectral index of α c = 0.5.
The ratio of two fluxes at two different wavelengths is simply which inserted into the magnitude equation gives a predicted conversion between the magnitude in UVOT band m i and in the Johnson V band.
Clearly, a conversion from the UVOT vv magnitude to the Johnson V magnitude is the most direct and most desirable conversion, because the UVOT vv band has close to the same central wavelength as the Johnson V band. Unfortunately, most Swift-UVOT observations were in the higher-energy bands, requiring significant conversions and introducing uncertainties into the analysis. For Swift-XRT sources with observations in multiple UVOT bands, we only used the converted V magnitude originating from the UVOT band closest to the Johnson V band, regarding that observation as most faithfully approximating the Johnson V magnitude.
Conversion between one spectral band and another with an assumed slope is fraught with uncertainty. However, because the average V magnitudes of the known pulsar and blazar samples are so distinct (the medians separated by five magnitudes, at least a factor of 100x in flux), these uncertainties are tolerable as pulsars are typically much dimmer than blazars in visual and UV photons. The V magnitudes in our training and research samples are not corrected for extinction, which would mainly impact pulsars within the galactic plane. However, using the hydrogen column densities of the pulsars in our sample and the n H vs A V scaling relation given in Güver &Özel (2009), it is unlikely that an entire sample of pulsars is dimmed enough to create a fictitious separation in V magnitudes from the blazar population that would bias our classification efforts. Figure 1 shows the converted V magnitude distribution (using V magnitude upper limits for those pulsars without solid estimates) of the known Fermi pulsars and blazars compared to the UV/optical counterparts in our unassociated sample.

Known and Unknown Samples
To train ML classification routines, we gathered a sample of 74 known gamma-ray pulsars (including radioloud, radio-quiet, and millisecond pulsars) and 635 known gamma-ray blazars. This sample was derived from the catalog Fermi -LAT pulsars and blazars (Abdo Table 1. X-ray counterparts to 4FGL unassociated sources with extreme X-ray photon indices, ΓX < −1 or ΓX > 5, excluded from the main ML classification effort. Possible stellar counterparts include spectral type and apparent magnitude from SIMBAD if available.
Overall, our training sample is selected for those pulsars and blazars for which gamma-ray, X-ray, and optical brightness data is available, which implicitly limits the training of our classifier to that subsample of pulsars and blazars that have notable gamma-ray, X-ray, and optical flux. This is especially impactful for pulsars, which have a wide range of intrinsic and extrinsic properties, including physical distance to our position in the galaxy. While we found that our training sample of pulsars is representative of catalogued gamma-ray pulsars in the Fermi -LAT pulsar lists in terms of period and spin-down rate, canonical pulsars tend to be younger and more energetic than the general pulsar population. For pulsars, this leads to a training sample that preferentially includes nearby or energetic pulsars. Fortunately, our research sample of unassociated sources is similarly limited by X-ray and UV/optical flux, and has a simi-lar fraction of subtypes of pulsars (65% canonical, 35% millisecond) compared to the overall Fermi -LAT pulsar sample (55% canonical, 45% millisecond). Any sources deemed likely pulsars by our classification efforts are astrophysically similar to our training sample.
Preferring parameters that are distance-independent, we expressed the X-ray and V-band brightness in terms of a ratio with the gamma-ray flux. While log F X /F γ is a simple ratio of two fluxes 2 , the conversion of V magnitude to flux was slightly more complicated. Using the magnitude equation we converted the V magnitudes to fluxes using a reference flux and magnitude m 2 and F 2 . However, because we are taking the logarithm of the ratio between V-band flux to gamma-ray flux, and because ML procedures first rescale and recenter each parameter, the exact reference flux and magnitude used herein are irrelevant. For our purposes, we adopt the conversion from Bessell et al. (1998).
The parameters for each training or research source include: • X-ray photon index, Γ X • Gamma-ray photon index, Γ γ (PL_Index in the 4FGL catalog) • The logarithm of gamma-ray flux, log F γ , in erg/s/cm 2 (Energy_Flux in the 4FGL catalog) • The logarithm of X-ray to gamma-ray flux ratio, • The significance of the curvature in the gammaray spectrum, henceforth simply curvature (PLEC_SigCurv in the 4FGL catalog) • The year-over-year gamma-ray variability index (Variability_Index in the 4FGL catalog) Normalized histograms of all parameters for both the known pulsars and blazars and the unassociated spectra are shown in Figure 2. For the 174 X-ray sources that are the unique XRT/UVOT possible counterpart to the respective unassociated targets in our research sample, the gamma-ray, X-ray, and UV/optical parameters are given in Table 3. For the unassociated targets with multiple possible counterparts in the uncertainty ellipse, the parameters are given in Table 4.
It is worth noting that while the histograms in Figure 2 show notable overlaps between the unassociated and known samples, there are still certain qualitative differences between the samples. The unassociated sources have systematically lower gamma-ray flux than the known blazars and pulsars, meaning that the variability and curvature estimates from the 4FGL catalog have weaker photon statistics. A bright Fermi -LAT blazar with significant variability on many timescales might not have observable variability if moved to a greater distance for no other reason than the Poissondistributed arrival statistics of the photons. Still, the current flux-limited samples of pulsars and blazars (for example, the flux-and redshift-limited sample of Fermi blazars in Ghisellini et al. 2017) suggests that there are some objects left out of the Fermi association lists simply due to being slightly dimmer.
Because there are many more known blazars than known pulsars in the training sample, we used Synthetic Minority Over-sampling Technique (SMOTE) (Chawla et al. 2002) to generate additional pulsars that mirror the distribution of real pulsars with a k-nearest neighbors approach. Previous classification efforts have shown that unbalanced training datasets can lead to classifiers that are biased against the underrepresented class (for example, Last et al. 2017). The result of the SMOTE expansion is a training set with an equal number of blazars and pulsars, with the artificial pulsars being generated via SMOTE using the real pulsar distribution. The final database is eventually split into training and validation subsamples which contain real known blazars, real known pulsars, and SMOTEgenerated 'known pulsars' that mirror the distributions of spectral properties of the real known pulsars.

Neural Network Design and Training
Building the NNC method, we constructed an approach where the unassociated sources will be assigned a blazar probability P bzr depending on their similarity with known pulsars and known blazars. P bzr = 0 denotes a 0% probability of the object being a blazar, while P bzr = 1 corresponds to 100% probability the object is a blazar. The NNC had seven input nodes (one for each parameter of the training dataset), one hidden layer with four neurons, and one output node returning the predicted blazar probability. In this research we use the MLPClassifier function of the scikit-learn Python package for NNC training and validation, training the NNC with the 'adam' optimization approach (Kingma & Ba 2014). To allow for validation and accuracy checks for the trained NNC, we took a random selection (20% of the training dataset) as a validation subsample, leaving the remainder of known pulsars and blazars as the training subsample. The random selection method StratifiedShuffleSplit was used so that the training and validation subsamples would have exactly the same proportion of pulsars and blazars.
To determine when to stop iterating and training the NNC, we relied on the Log-Loss parameter, an error measure in binary classification approaches that is similar to a likelihood estimator. For a sample with true classifications y i = (0, 1) and predicted classification p i ∈ [0, 1], the Log-Loss is given by At each iteration step in training the NNC, we calculate L log of both the training and validation subsamples. As training continues, the training subsample L log should continuously decrease as the NNC approaches a more exacting fit. However, at some point the validation L log begins to increase as the NNC starts to overfit the training sample and lose predictive capabilities on unknown data points. At this point, training is stopped. Figure 3 shows how the NNC training is stopped at approximately 2000 iterations when the validation L log levels out.  After using the training subsample to construct the NNC, we passed the validation subsample through the NNC and recorded the predicted blazar probabilities. Ideally, the NNC would predict blazar probabilities of P bzr = 0 or P bzr = 1 for the validation pulsars or blazars. A validation score for the NNC method depends on the cutoff used to determine what constitutes a "likely" pulsar or plazar. For example, a 90% cutoff would designate any source with P bzr < 0.1 a pulsar, with P bzr > 0.9 a blazar, while a 99% cutoff would only capture sources with 99% confident classifications.
An optimally trained NNC should result in P bzr as close to 0 or 1 as possible for the validation subsample. Typically, NNCs are judged based on the proportion of validation scores above some cutoff value; in a two-class example, the default cutoff for a 'correct' classification is normally 0.5. However, this single cutoff score does not investigate the degree of confidence with which the NNC is classifying the validation sources; a NNC which grades all validation pulsars at P bzr = 0.3 and all validation blazars at P bzr = 0.7 would have a 100% validation score with a cutoff of 0.5 but a 0% score if one strives for greater confidence in classification. Figure 4 shows a generalized approach for determining both the reliability and confidence of classifications from a NNC. Scanning through a logarithmic range of possible cutoff scores approaching unity (with P psr = 1−P bzr for pulsars), the figure shows how the fraction of known pulsars and blazars correctly classified as such during the validation step changes depending on the P bzr cutoff used. This general approach illustrates progress classifying pulsars and blazars at higher degrees of confidence upon adding additional spectral parameters to the samples.
Applying the default 50% cutoff used in previous classification efforts (i.e. a validation blazar must just have a validation P bzr > 0.5 to be considered "correctly" classified) returns an overall accuracy of 99.2% for pulsars and blazars. However, a major improvement of this NNC approach to previous classification efforts (Kaur et al. 2019;Kerby et al. 2021) is its capability to classify validation pulsars and blazars with over 99% confidence. As shown in Figure 4, over 90% of validation pulsars and blazars are classified with greater than 99% confidence in the correct categories, buoying our hope that this NNC routine can pick likely blazars and pulsars from the unassociated sample.
Calculating importances of the different parameters in a neural network ML approach is slightly more difficult than in other related classification schemes. Neurons in an NNC use an activation function which means that certain parameters may not be considered at all for cer-  . The fraction of validation pulsars (blue) and blazars (red) classified with P bzr above different cutoff values after training the NNC. The NNC classifies validation blazars with greater P bzr than validation pulsars, with no validation pulsars having P bzr < 10 −3 but many validation blazars having P bzr > 0.999

XrayInd GammaInd
LogFG LogFXFG LogVar LogCurve LogFVFG Hidden Figure 5. A qualitative visualization of the importances of the spectral features in our NNC routine. The columns represent the four neurons in our single hidden layer, the rows of the upper plot representing the seven spectral features used herein. Squares that are shaded darker have linearly heavier weights in the NNC, and therefore play a more important role in discriminating pulsars from blazars.
tain entries. Still, the NNC is described by a weight matrix, and for a network with only a single hidden layer these weights are appropriate measures of the importance of different features. Figure 5 shows a visualization of the weights of the different parameters in our NNC approach, with darker shaded squares representing more impactful features in classifying blazars and pulsars.

Unassociated Classification Results
After training and validating the NNC, we recorded output P bzr values for each of spectra in the unassociated sample. Of the 205 total fully examined sources, we found 157 with P bzr > 0.99 and 18 with P bzr < 0.01. As some of the unassociated targets have more than one notable X-ray source in the gamma-ray uncertainty ellipse, additional work is needed to decipher the results around those targets including which if any X-ray possible counterpart is the actual partner to the gammaray emission. However, of the 174 X-ray sources that are the unique XRT/UVOT possible counterpart to the unassociated gamma-ray source, 14 have P bzr < 0.01 and 132 have P bzr > 0.99. These portions of our research dataset are highly likely pulsar and blazar candidates. The results of the NNC classification on the 4FGL sources with just a single possible X-ray counterpart are given in the last column of Table 3, organized into likely pulsars (P bzr < 0.01), likely blazars (P bzr > 0.99) and ambiguous sources.
Several of the 4FGL targets examined herein have since been the subject of subsequent discoveries that shed additional light on their nature.  discusses the detection of a redback pulsar near 4FGL J0955.3-3949, which matches both for X-ray spectral parameters and in location on the sky of the X-ray counterpart in our analysis. 4FGL J0212.1+5321 has also been connected to a millisecond redback pulsar (Linares et al. 2017;Li et al. 2016). 4FGL J2039.5-5617 has recently been identified as a redback pulsar (Clark et al. 2021), as has 4FGL J1306-4035 after observations by the Parkes radio telescope (Keane et al. 2018), and 4FGL J1304.4+1203 is unassociated in 4FGL but classified as a pulsar in 4FGL-DR2 (Ballet et al. 2020). These sources are classified as likely pulsars in this work, showing that our NNC can independently predict classifications for sources that in other works are verified. This increases our confidence that our classifier can point towards other promising sources to examine, and that the link between Swift X-ray/UV/optical and Fermi -LAT gamma-ray sources is valuable and well-founded. Additionally, 4FGL J0838.7-2827 and J0523.3-2527, sources that remained ambiguous in this work, have also since been identified as redback pulsars (Halpern et al. 2017;Strader et al. 2014).
To gain additional insight into 4FGL unassociated sources with multiple notable X-ray sources, we leveraged the much higher spatial resolution of the Swift telescopes to do a coordinate search at the positions of X-ray counterparts. For many lower-energy counterparts, this search returned a coincident catalogued object and provided useful information to decide if the XRT/UVOT source is a likely counterpart to the unassociated target. Table 4 includes the NNC classification results of these possible counterparts grouped by 4FGL target, and Table 5 lists the results of our spatial cross-reference. Several possible counterparts are coincident with catalogued stars, X-ray binaries, or galaxies, differences that might be useful for finding the most likely counterpart among multiple near a single 4FGL target. Interestingly, many of the pairs of possible counterparts have very similar P bzr values, suggesting that the 4FGL source in question may be a pulsar or blazar regardless of which X-ray source is truly linked with the gamma-ray emission.
The classification of the unassociated sample is more bifurcated in terms of P bzr compared to previous works (Kaur et al. 2019;Kerby et al. 2021); most of the unassociated sources have P bzr very close to 0 or 1, showing that the NNC makes confident predictions of blazar or pulsar class membership. The NNC's high validation accuracy and confidence suggest that the unassociated sources can be classified properly if they are pulsars and blazars similar to the known pulsars and blazars in our training sample. Even outside the context of the NNC as a whole, the histograms in Figure 2 and the feature weights in Figure 5 show that log F X /F γ and log F V /F γ , flux ratios introduced in this work, are important parameters to distinguish pulsars from blazars, with pulsars having systematically lower values of both.
While previous approaches have classified unassociated sources with gamma-ray spectral parameters alone, Figure 2 shows that X-ray/UV/optical properties of counterparts to 4FGL targets can dramatically improve the discrimination of pulsars from blazars. While the unassociated sources have lower gamma-ray flux than either of the known pulsar/blazar distributions, and the gamma-ray photon index is not a particularly useful discriminatory variable in any capacity, the variability and spectral curvature measures in the Fermi catalog send mixed signals. The unassociated sources have the low variability expected of pulsars, but the low spectral curvature expected of blazars, probably due to their inherently lower photon statistics. Our addition of several new distant-independent properties via Swift-XRT and -UVOT analysis adds three new features; Figure 5 shows that log F V /F G and log F X /F G are features of moderate to major importance in the NNC weight matrix, as is the X-ray photon index. This reinforces the importance of systematic and methodical follow-up observations of unassociated gamma-ray sources at lower energies to uncover their astrophysical justification.
Overall, the spectral analysis and classification using Swift-XRT and -UVOT data herein is valuable not only for characterizing the Fermi unassociated sources, but also for guiding follow-up observations with likely pulsars and blazars for targeted searches while supporting numerous additions to catalogs of known gamma-ray pulsars and blazars.

Comparison with Random Forest Classifier
A random forest (RF) classifier uses an array of decision trees to classify an object into one of several categories. Each individual decision tree consists of several inequalities using the different parameters of a dataset in a "choose your own adventure"-style series of judgements; compounded many times, a RF classifies an unknown object based on the fraction of the constituent trees giving each classification.
Previous works by this group on earlier catalogs of Fermi unassociated sources (Kaur et al. 2019;Kerby et al. 2021) have used RF classifiers to discern likely pulsars from likely blazars among the Fermi -LAT unassociated sources. While we have used an NNC approach in this work, it is worthwhile to compare the NNC P bzr output to the values produced by a RF classifier. Though the NNC and RF approaches make independent predictions for P bzr on the unassociated sources, the two approaches giving similar outputs on the same dataset would increase confidence in the reproducibility of our results even with different classification approaches. Additionally, comparing the two methods illuminates whether one is more descriptive or predictive in its P bzr output on our multiwavelength data.
We applied the RF classifier developed in Kerby et al. (2021) to the 4FGL sample of this paper. Figure 6 shows that while the NNC and RF both classify many sources with P bzr close to 0 or 1, the NNC classifies many sources as high or low P bzr that the RF leaves ambiguous. Indeed, the NNC results only have a few sources with P bzr not close to 0 or 1, and these sources almost uniformly have RF P bzr values around 0.8. This trend suggests that the NNC is more sensitive to pulsars than the RF approach, classifying sources as likely pulsars that in the RF method are ambiguous. The more bifurcated nature of the NNC P bzr values compared to the RF results suggests that our NNC approach tends towards more confident classification and leaves fewer sources ambiguous, producing more immediately verifiable results, rather then hedging with ambiguity.

Summary and Next Steps
In this work we have classified 174 unique gammaray/X-ray/UV/optical sources into 14 likely pulsars, 132 . Comparing the P bzr scores for the 4FGL unassociated sample, using the exact same training and unknown datasets. The histograms of the two datasets have limited Y-axes to show the sparsely-populated bins between the two extremes. In both cases, the number of sources with P bzr close to unity is over 200.
likely blazars, and 28 ambiguous. Using Swift-XRT and -UVOT observations, we built a collection of results, presented in Tables 3 and 4, representing a significant collection of observations within the uncertainty ellipse of the 4FGL unassociated gamma-ray sources. It is likely that the X-ray/UV/optical sources described herein have the same astrophysical origin as the gammarays described in the Fermi -LAT unassociated catalog, so the unique correspondences in Table 3 describe the lower-energy spectra of the astronomical objects behind the gamma-ray emission of Fermi unassociated sources. Next, we built a neural network classification approach to use spectral information to divide the unassociated sample into likely blazars and pulsars using samples of known gamma-ray blazars and pulsars, reaching higher accuracy on validation subsamples and greater confidence in classification of unassociated sources than previous approaches. Of the 174 unique gamma-ray/Xray/UV/optical spectra constructed and described in Table 3, 132 are P bzr > 0.99 likely blazars and 14 are P bzr < 0.01 likely pulsars. Leveraging the advantages of multiwavelength analysis, our new subsamples of likely pulsars and blazars can expand known gamma-ray pulsar and blazar catalogs to include sources with lower gamma-ray luminosity that were previously unassociated.
As Swift continues its observation campaign of the Fermi -LAT unassociated sources, additional X-ray sources will be detected around previously unobserved targets. Planned follow-up observations across the electromagnetic spectrum can continue to investigate interesting sources discovered herein, including likely pulsars and ambiguously classified objects. These additional observations could validate likely pulsar classification with radio observations near the X-ray source to detect radio pulsations, or investigate ambiguous sources in greater detail to discern their true origin.
The subsamples of likely pulsars and blazars classified with our NNC approach are prime candidates for inclusion in population studies of gamma-ray pulsars and blazars. For example, using archival infrared observations of the likely blazars should allow for classification into likely BL Lac or likely FSRQ subsets, illuminating the biases and drawbacks of the Fermi blazar catalogs used to investigate blazars as a class of AGN. The likely pulsar subset, if validated with radio pulsation searches, would expand the list of known gamma-ray pulsars by almost 10%.
Software: Astropy (The Astropy Collaboration et al. (Harris et al. 2020), Matplotlib (Hunter 2007), scikitlearn (Pedregosa et al. 2011), FTools (Blackburn 1995 ACKNOWLEDGMENTS This research has made use of data and/or software provided by the High Energy Astrophysics Science Archive Research Center (HEASARC), which is a service of the Astrophysics Science Division at NASA/GSFC. We gratefully acknowledge the support of NASA grants 80NSSC17K0752 and 80NSSC18K1730. E. Ferrara is supported by NASA under award number 80GSFC17M0002.

2013), numpy
Fermi research at NRL is supported by NASA. Table 3. Fermi-LAT features for the unassociated sample investigated in this work, along with Swift-XRT and -UVOT parameters for likely X-ray/UV/optical counterpart. Only unassociated sources with a single possible counterpart are included in this table, and sources are organized based on P bzr , listing likely pulsars and likely blazars separately. The extraction and derivation of the parameters here are described in sections 3 and 2. The mV estimates were produced using the noted UVOT filter, the closest available to the V band central wavelength.  Table 3 continued  Table 3 continued  Table 3 continued  Table 3 continued  Table 3 continued  Table 3 continued  Table 4. Fermi-LAT features for the unassociated sample investigated in this work, along with Swift-XRT and -UVOT parameters for all possible Xray/UV/optical counterparts. Only unassociated sources with multiple notable X-ray sources within the Fermi-LAT uncertainly ellipse are included here.
The extraction and derivation of the parameters here are described in sections 3 and 2. The mV estimates were produced using the noted UVOT filter, the closest available to the V band central wavelength. Two sources, noted with *asterisks*, are coincident with dim catalogued stars. UV/optical emission from those stars is probably not related to a pulsar or blazar, and thus the P bzr values should not be trusted.  Table 5. Results from a SIMBAD position cross-reference for the possible counterparts to unassociated sources with multiple notable excesses in the gamma-ray uncertainty ellipse, helping discriminate excesses that might be eliminated from consideration as possible counterparts with additional astrophysical information or description.
Target name XRT excess name Cross-reference results