1 Yale University, New Haven, CT; 2Rice University, Houston, TX;
3University of Texas Health Science Center, Houston, TX

# Local Contrastive Learning for Medical Image Recognition

Syed A. Rizvi  BS1 Ruixiang Tang  BSC2 Xiaoqian Jiang  PhD3
Xiaotian Ma  MS3 Xia Hu  PhD2

## Abstract

The proliferation of Deep Learning (DL)-based methods for radiographic image
analysis has created a great demand for expert-labeled radiology data. Recent
self-supervised frameworks have alleviated the need for expert labeling by
obtaining supervision from associated radiology reports. These frameworks,
however, struggle to distinguish the subtle differences between different
pathologies in medical images. Additionally, many of them do not provide
interpretation between image regions and text, making it difficult for
radiologists to assess model predictions. In this work, we propose Local
Region Contrastive Learning (LRCLR), a flexible fine-tuning framework that
adds layers for significant image region selection as well as cross-modality
interaction. Our results on an external validation set of chest x-rays suggest
that LRCLR identifies significant local image regions and provides meaningful
interpretation against radiology text while improving zero-shot performance on
several chest x-ray medical findings.

## 1 Introduction

Advancements in medical imaging technologies have accelerated the creation of
large radiographic image datasets, opening new possibilities for image
analysis in patient care. This development, however, has been in parallel with
an increasing burden on radiologists to interpret larger numbers of medical
images 1. Deep Learning (DL) provides a promising solution for automating the
analysis of these large datasets. Deep Neural Networks (DNNs) have
demonstrated strong representational learning power at multiple layers of
abstraction, without relying on handcrafted features or domain expertise2. In
the medical domain, DL-based systems have shown strong performance on
downstream tasks such as disease classification 3; 4 and pathology
localization 5; 6; 7. Annotating large radiology datasets for training DL
models, however, is cost-prohibitive at scale due to the cost and expertise
required to annotate medical images.

Self-supervised vision-language frameworks address this challenge by obtaining
supervision on images from associated text 8, which comes in the form of
associated radiology reports. In Computer Vision, self-supervised Vision-
language pretraining frameworks 9; 10; 11; 8 train on large corpora of paired
image and text data, learning alignment between image and text embeddings.
Fine-tuning can then be done for downstream tasks such as visual question
answering and image-text retrieval 11. In the medical domain, self-supervised
vision-language architectures have shown promising results in classification
12 and medical image-report retrieval 13. Contrastive frameworks 14; 15; 16 in
particular have shown strong performance on downstream tasks, including
medical image-text retrieval and zero-shot classification, where a model is
able to predict new pathologies that were not explicitly annotated in the
training set 17; 16. Contrastive frameworks operate by encoding images and
text into a shared latent space and define an objective that maximizes
agreement between paired images and text 18; 19.

Although powerful, these architectures may struggle to differentiate radiology
images, which have less variation between different images than natural
images. Medical observations often manifest in small regions of an image, and
often many regions of the image contain uninformative background regions which
do not contribute meaningful information to the global embedding. Several
works have focused on learning both global and local features in input images
to emphasize discriminative information contained in local regions. In natural
images, TransFG 20 analyzed the self-attention maps created by Vision
Transformers (ViTs) 21 to identify significant local regions for fine-grained
visual recognition tasks. Li et al. 22 proposed a visual-semantic reasoning
model to enhance the relationship between local image regions before
contrasting image representations against text data. In radiology image
analysis, GLoRIA 14 proposed a local contrastive loss in addition to the
normal contrastive objective between the global image and text embeddings.
Attention weights are learned between image region embeddings and individual
word embeddings, generating context-aware local image embeddings which are
then used in a contrastive loss against word embeddings. In the attention-
weighted image representations, however, unimportant background image regions
are still contributing information to the contextualized local image
embedding. Additionally, the frameworks above are standalone architectures,
not compatible with other pretrained model architectures.

In addition to downstream performance, the interpretability of computational
systems for medical image analysis tasks is crucial for providing context
behind model predictions. Contrastive architectures 8; 19; 16 gain zero-shot
capabilities when pretrained on large amounts of data, however, cannot give
reasoning or significant features associated with the prediction. In
radiology, computational systems often aid radiologists in analyzing chest
radiographs, necessitating a further degree of cooperation between the
radiologist and the computational model than a single prediction. Highlighting
significant image regions associated with a given text prompt can provide a
visual indicator for the radiologist to check the model’s prediction and
increase radiologist-computer cooperation. In a similar vein, GLoRIA 14
provides attention scores between image regions and text prompts. GLoRIA
cannot, however, highlight significant image regions without contrasting all
image regions against a text prompt. In addition, the architecture is not
immediately compatible with other pretrained contrastive architectures.
Additional modularity and flexibility are needed in order to open up
possibilities for local region selection and local contrastive learning to
other pretrained contrastive models.

In this work, we propose LRCLR, a fine-tuning module aimed at introducing
interpretable local region selection while making minimal changes to existing
pretrained architectures. LRCLR consists of a region selection module, which
selects significant regions in input images solely based on self-attention
matrices from the transformer-based image encoder. Selected regions are passed
along with text tokens to a cross-modal transformer, which contextualizes
significant image regions with texts through self-attention. A local
contrastive objective is added to the global objective, contrasting the image
and text class tokens output from the cross-modal transformer. LRCLR is
compatible with many vision-language architectures as an add-on module during
fine-tuning. We evaluate our approach in a zero-shot classification setting on
an external dataset of chest x-rays and study the interpretability of rare
chest x-ray findings. We demonstrate improved zero-shot classification
performance on rare chest x-ray findings and show interpretability results
between significant local regions and radiology report text.

## 2 Preliminaries

#### Medical Image Recognition

Medical Image Recognition is a diagnostic tool for clinicians in patient care,
involving inspecting multiple modalities of medical image data for patterns of
pathology and medical observations 23. With advancements in imaging
technologies as well as multi-modal image data being generated in patient
examinations, the burden on clinicians has increased to analyze patient
information. Initial studies aimed at using computational power to assist in
medical image analysis emerged as early as the 1960s, based on image analysis
and early Machine Learning (ML) techniques 24. The emergence of Deep Learning,
paired with increased computational power and availability of medical image
data, has accelerated the field of Machine Learning-based medical image
analysis, based on the representation learning ability of Deep Neural Networks
(DNNs) at multiple levels of abstractions 2. Recent DL works based on Deep
Convolutional Neural Networks (CNNs) 25 and Transformers 26 have shown strong
performance on medical image classification 27 and segmentation 28 tasks.

#### Radiology Data

Radiology is a branch of medicine that utilizes diagnostic imaging techniques
to help treat patients, including x-rays, computed tomography (CT) scans,
magnetic resonance images (MRIs), and ultrasounds (US) 29. Common tasks for
ML-based systems in radiology include segmentation, disease detection, brain
function and neurological activity analysis from fMRIs, image-retrieval
systems, and radiology report text analysis 30.

#### Self-Supervised Vision-Language Architectures

Pretraining multimodal vision-language architectures has been shown to improve
performance on downstream vision-language tasks. Vision-language frameworks
can be divided into several categories, depending on their architectures and
encoding strategies. Dual-branch fusion frameworks 9 use two parallel encoders
for different modalities, with cross-modal attention mechanisms for capturing
inter-modal interactions. Single-branch fusion architectures 10; 11 define a
single sequence-based architecture, usually a transformer, which accepts a
sequence of concatenated image and text tokens. The self-attention within the
architecture captures cross-modal interactions, with different objectives
defined on the model output, such as masked language modeling, masked region
modeling, or word-region alignment 10. Dual encoder-based architectures 8
utilize a separate encoder to map images and text into a shared embedding
space, where a contrastive objective is used to align their representations.
This enables zero-shot capabilities at inference time by introducing new text
prompts and images to the trained model.

#### Zero-shot Learning

Zero-shot learning is a paradigm which aims to classify images belonging to
classes that were not explicitly labeled in the training dataset. Without
training labels, supervision for images is obtained from other data
modalities, such as an associated text corpus 31; 32; 8. By learning alignment
between images and text through a contrastive learning objective, new text
prompts can be introduced during inference to find the most similar zero-shot
label to a given image. Recent works have shown promising results on zero-shot
learning in natural images 8 and medical image classification tasks 14; 16.

Figure 1: Overview of Local Contrastive Learning (LRCLR) framework.

## 3 Method

Datasets. Self-supervised fine-tuning was done on the MIMIC-CXR dataset, a
public dataset of 377,100 chest radiographs drawn from 227,835 radiographic
studies 33. Each study comes with a corresponding radiology text report
containing multiple sections written by radiologists: examination, indication,
impression, findings, techniques, and comparison. We fine-tune the model for
20 epochs on the MIMIC-CXR dataset 33. Image-text pairs are constructed by
taking the chest radiograph and the impressions section of the corresponding
radiology report. The zero-shot evaluation was done on the CheXpert dataset, a
public chest X-ray dataset commonly used to train and evaluate the performance
of radiology image classifiers 34. The CheXpert dataset consists of 224,316
chest radiographs taken from 65,240 patients, with a free-text radiology
report associated with each radiograph. Each radiograph is also labeled for
the presence of 14 different medical observations. Following 16, we evaluate
the zero-shot performance of our models on the CheXpert test dataset, which
consists of 500 chest radiographs.

Image and Text Encoding. An overview of the LRCLR framework is presented in
Figure 1. The architecture consists of a Vision Transformer (ViT) 21 based
image encoder and a transformer text encoder. We initialize our image encoder
with the weights of a ViT pretrained on ImageNet 35, which takes as input a
224x224 resolution chest radiograph. The image encoder maps the input image
into a global embedding which summarizes the information in the image. Local
features are encoded in the patch embeddings, which capture information about
subregions within the image. For the text encoder, we use a 12-layer
transformer with a hidden dimension of 512 and 8 attention heads, following
the parameters set in 16. The initialization scheme for the text encoder
follows the original CLIP architecture, and a maximum context length of 77
tokens for the text extracted from the findings section of the radiology
report. 8.

Figure 2: Overview of local region selection module.

Local Region Selection. Accurately selecting local image regions is critical
for differentiating chest radiographs. Medical observations manifest
themselves in small regions of chest x-rays, resulting in subtle variations
between different chest radiographs. We draw inspiration from fine-grained
visual recognition tasks in Computer Vision, which aim to classify subclasses
with small inter-class variation within an object category, such as birds 36
or cars 37. Methods based on region-wise annotations are unfeasible given the
scarcity and cost of annotating chest radiographs. Therefore, we turn to
unsupervised methods for local region detection. Following 20, we utilize the
self-attention matrices learned as a part of the ViT image encoder to select
important local regions. Attention matrices are calculated within each layer
of the image encoder as a product of query and key vectors, which are obtained
from the input sequence of vectors. Given $N$ input vectors
$h_{1}^{l},...,h_{N}^{l}$ as a matrix $H^{l}$ for a given layer:

$Q=H^{l}\cdot W_{query}^{l}$ (1) $K=H^{l}\cdot W_{key}^{l}$ (2)
$a_{l}=softmax(\frac{QK^{T}}{\sqrt{d_{k}}})V$ (3)

Where $W_{query}^{l}$ and $W_{key}^{l}$ are learned projection matrices for
query and key vectors, and $\sqrt{d_{k}}$ is the root of the hidden dimension
of the query and key vectors. We extract self-attention matrices from all L
layers and K heads of the image encoder, and perform matrix multiplication to
obtain a representative matrix for each attention head:

$a_{k}^{final}={\displaystyle\prod_{l=0}^{L}a_{k}^{l}}$ (4)

where k is one of the K attention heads. $a_{final}$ captures how information
about local regions propagates throughout the layers of the ViT, giving a
better attention-based selection compared to a single layer of attention 20.
We take the maximum attended-to index for each attention head and select the
corresponding K image regions. By keeping only selected image regions as well
as the class token, we allow the model to focus on discriminative regions of
images while also keeping the global image context.

Cross-Modality Transformer. Selecting local regions allows the model to focus
on important local regions within an image. However, it cannot by itself
provide any interpretation of the relationship between significant image
regions and the accompanying radiology text. Additionally, the local regions
may be underrepresented in the global image encoding, which contributes to the
alignment between chest radiographs and radiology text. We therefore propose a
cross-modality transformer module as a second add-on for fine-tuning,
accepting as input a sequence of selected image region tokens and text tokens.
A local contrastive loss 8; 15 is calculated between the output image class
token $v_{l}^{i}$ and text class token $t_{l}^{i}$ from the cross-modality
layer:

$L_{local}=\sum_{i=1}^{N}-log\left(\frac{e^{v_{l}^{i}\cdot
t_{l}^{i}}}{\sum_{k=1}^{N}e^{v_{l}^{i}\cdot t_{l}^{k}}}\right)$ (5)

By jointly encoding image regions and text and aligning their representations
through local contrastive learning, we emphasize the interaction between
selected image regions and text. Our overall loss function can then be
calculated as:

$L=L_{global}+\lambda L_{local}$ (6)

Where $L_{local}$ is the global contrastive loss between the global image
class token $v_{g}^{i}$ and global text class token $t_{g}^{i}$ output from
the region selection module. Parameter $\lambda$ is left as a hyperparameter
for weighing local loss relative to global loss. We use $\lambda=0.5$ for our
experiments, chosen through experimentation on our zero-shot task.

Model |  Specificity |  Precision |  Recall |  F1-Score |  AUC
---|---|---|---|---|---
GLoRIA |  61.81 |  28.13 |  65.48 |  34.18 | 55.39
ConVIRT |  59.84 |  25.32 |  67.27 |  32.56 | 58.64
CheXzero |  69.46 |  40.79 |  65.2 |  45.72 | 76.24
CheXzero + LRCLR |  66.86 |  41.04 |  64.25 |  46.79 | 78.66
Table 1: Summary of zero-shot classification performance on the CheXpert test
set. Metrics are calculated over all 14 CheXpert findings. Figure 3: AUC
performance across all 14 CheXpert medical observations.

## 4 Experiments & Results

Zero-Shot Classification. We evaluate the effectiveness of our proposed fine-
tuning module in a zero-shot classification setting on the CheXpert test
dataset following the methodology outlined in Section 3. We first take image-
text pairs of the chest x-rays and zero-shot classification labels associated
with the CheXzero dataset, and encode them using the fine-tuned model. The
encoded global image and text representations are then contrasted against each
other to determine which text labels the image most corresponds with.
Additionally, we extract attention weights from the cross-modality transformer
module to inspect the relationship between selected image regions and zero-
shot labels. To evaluate the performance of baseline models and our proposed
method, predictions are obtained from each model on the CheXpert test dataset,
and metrics are calculated for each of the 14 findings. We report the average
specificity, precision, recall, F1-score, and AUROC for each of the models
across all 14 medical observations. For CheXzero 16, we use the pretrained
CLIP checkpoint from 8, and fine-tune on the MIMIC-CXR dataset 33. For the
GLoRIA 14 baselines, we load the publicly available model and perform zero-
shot evaluation on the CheXpert test dataset. For ConVIRT, since there is no
publicly available code, we implement the method according to 15 and evaluate
it on the CheXpert dataset.

Figure 4: Examples of chest radiographs with top 2 selected regions
highlighted. A single highlighted region indicates that two attention heads
highlighted the same image region with a high score.

We present the comparison of our proposed fine-tuning module against baseline
methods in Table 1. Based on our experimental results, we observe that the
addition of our fine-tuning module led to equal or improved performance in
three of the five evaluated metrics, including a 2.42 improvement in AUROC
across all 14 medical observations in the CheXpert test dataset. A comparison
of AUROC across each of the 14 findings (Figure 3) shows that the region
selection and local contrastive learning resulted in improved performance on
seven of the medical observations, most notably a 27.2 increase in AUC on
Fracture observations, a 7.6 increase on No Finding observations, and a 13.0
increase on Pneumothorax observations.

From these results, we infer that our local region selection module is
effectively selecting significant image regions in the chest radiographs,
allowing the model to focus on more discriminative regions between different
images. We note that the number of selected local regions corresponds to the
number of attention heads in the image encoder, and that two attention heads
may attend highly to the same image region.

In Figure 4, we visualize the attention weights between selected local regions
and zero-shot labels extracted from the cross-modality transformer module. The
highlighted regions indicate significant image regions selected by the region
selection module, colored by attention score from lowest (dark blue) to
highest (dark red). We observe that the attention scores and selected regions
correctly identify significant image regions corresponding to the text prompt,
and provide an ordering when comparing image regions to a given text prompt.
We emphasize that, unlike previous works, our fine-tuning modules are
compatible with many contrastive attention-based architectures, and can be
added in for additional interpretability across image regions and text.
Previous works 38 have studied the interpretability of attention coefficients
in deep networks, and while the scores may not provide complete
interpretation, a well-trained contrastive model can still highlight
significant image regions corresponding to a text prompt.

Figure 5: Comparison of selected regions from the LRCLR framework and baseline
CheXzero model. Baseline selected regions were obtained by extracting
attention weights from the ViT image encoder and following the same procedure
as outlined in Section 3.

Case Study: Fracture. We conducted an analysis of our fine-tuning module on
the fracture medical observation in the CheXpert test set. We identified the
largest improvement in AUROC over baseline methods on this finding. Attention
visualizations, presented in Figure 4, revealed that the trained model often
focused on the chest rib regions in chest radiographs - a common site for
fractures. Such rib fractures can result in additional complications for
patients 39, and chest radiography is a reliable method for their diagnosis.
Our computational aid, the cross-modality module, can provide valuable
interpretability between significant regions and text, and aid radiologists in
identifying critical regions that require further analysis.

We compare our selected image regions to baseline selected regions in Figure
5. Baseline selected regions were obtained by extracting attention weights
from the Vision Transformer image encoder in the baseline CheXzero model, and
following the procedure outlined in Section 3 to obtain highlighted regions
per attention head. We rank the baseline selected regions by the associated
attention score from the corresponding attention head, and choose the top two
scored image regions for comparison to our method. We note that baseline
selected regions do not effectively highlight the rib region in the chest
radiographs, and instead focus on other outer portions of the image.

Case Study: Lung Opacity. We additionally evaluate our fine-tuning module on
the lung opacity medical observation in the CheXpert test set, which also saw
improved performance compared to baseline zero-shot performance. Attention
visualizations from our method 4 highlight regions in the lower lung, where
uncertain areas of white material are present in the radiograph. These regions
may be indications of lung complications. In comparison to baseline selected
regions shown in Figure 5, we observe that our method focuses more on unclear
lung regions, whereas the baseline highlights regions on the edges of the
radiograph.

While insightful, we acknowledge that the selected regions highlighted by our
method do not always correspond well to the intended text prompt. Highlighted
regions may focus on different regions within the radiograph, particularly in
radiographs with many conditions present or in noisy images. We do, however,
note that by ranking image regions with attention corresponding to the given
text prompt, we can filter out some of these uninsightful selected regions.

### 6 Related Works

In this section, we review related literature on contrastive learning methods
and local representational learning.

Contrastive vision-language frameworks. CLIP 8 extracts vision and language
embeddings using image and text encoders which are trained by contrastive
learning by a symmetric cross-entropy loss. Similarity scores are calculated
between image and text embeddings, which are used after pretraining for
downstream tasks such as zero-shot classification and image-to-text retrieval.
ConVIRT 15 utilized the idea of contrastive learning to pretrain on medical
image-text pairs. MedCLIP 40 decouples paired images and texts and uses soft
targets of semantic similarities to learn from unpaired medical images and
text. CheXzero 16 is a direct application of the CLIP model on large-scale
chest X-ray datasets to enable zero-shot classification of unseen findings in
images. Seibold et al. 41 improved the CLIP model by adding a local
contrastive loss function that contrasts local image features to sentence
features in the radiology reports. LoVT 42 extracts localized representations
by the local alignment between output feature maps (image sub-region features)
and max pooling of token embeddings (sentence embeddings). BioViL 43 develops
a new vision-and-language pre-training model and proposes a phrase grounding
task to evaluate the local alignments between images and texts. MGCA 44
further uses three levels of embedding alignments, i.e., pathological region,
instance, and disease, to combine local, global, and cluster representations
for downstream tasks.

Global-Local Representational Learning. Chaitanya et al. proposed a framework
45 that uses self-supervised learning (SSL) and contrastive learning for semi-
supervised segmentation of volumetric medical images with limited annotations.
They leverage structural similarity across images and a local contrastive loss
to learn distinctive local representations for per-pixel segmentation. GLoRIA
14 learns a global image and text embedding, and additionally learns attention
weights between individual region and word embeddings through a local
contrastive objective. Although similar to our work, our proposed method
differs from GLoRIA in that we perform self attention-based local region
selection to discard unimportant local regions before the cross-modal
transformer module and local contrastive loss. Additionally, we formulate our
framework as a flexible fine-tuning module, which can be added to any
contrastive-based architecture for added interpretability.

### 7 Conclusion

We introduce a flexible fine-tuning module for contrastive transformer-based
frameworks, designed to enhance cross-modal interpretability in chest
radiograph analysis tasks. Our two-part module, named LRCLR, comprises a
region selection layer and a cross-modal interaction layer, offering
interpretable scores between image regions and associated text. Visualizations
of attention scores and selected regions exhibit the model’s ability to
identify and emphasize crucial image regions, thereby enabling cross-modal
interpretability for radiologists analyzing chest radiographs. Our work
highlights the versatility of fine-tuning to incorporate interpretability in
contrastive architectures. We suggest extending local region selection to non-
attention-based image encoders and exploring various unsupervised region
selection procedures for attention-based methods to improve this approach.

## Acknowledgement

XJ is CPRIT Scholar in Cancer Research (RR180012), and he was supported in
part by Christopher Sarofim Family Professorship, UT Stars award, UTHealth
startup, the National Institute of Health (NIH) under award number
R01AG066749, R01LM013712, and U01TR002062, and the National Science Foundation
(NSF) #2124789

## References

  * 1 Bhargavan M, Kaye AH, Forman HP, Sunshine JH.  Workload of radiologists in United States in 2006–2007 and trends since 1991–1992.  Radiology. 2009;252(2):458-67.
  * 2 LeCun Y, Bengio Y, Hinton G.  Deep learning.  nature. 2015;521(7553):436-44.
  * 3 Esteva A, Kuprel B, Novoa RA, Ko J, Swetter SM, Blau HM, et al.  Dermatologist-level classification of skin cancer with deep neural networks.  nature. 2017;542(7639):115-8.
  * 4 McKinney SM, Sieniek M, Godbole V, Godwin J, Antropova N, Ashrafian H, et al.  International evaluation of an AI system for breast cancer screening.  Nature. 2020;577(7788):89-94.
  * 5 Chen H, Ni D, Qin J, Li S, Yang X, Wang T, et al.  Standard plane localization in fetal ultrasound via domain transferred deep neural networks.  IEEE journal of biomedical and health informatics. 2015;19(5):1627-36.
  * 6 Ghesu FC, Georgescu B, Mansi T, Neumann D, Hornegger J, Comaniciu D.  An artificial agent for anatomical landmark detection in medical images.  In: Medical Image Computing and Computer-Assisted Intervention-MICCAI 2016: 19th International Conference, Athens, Greece, October 17-21, 2016, Proceedings, Part III 19. Springer; 2016. p. 229-37.
  * 7 Rajpurkar P, Irvin J, Zhu K, Yang B, Mehta H, Duan T, et al.  Chexnet: Radiologist-level pneumonia detection on chest x-rays with deep learning.  arXiv preprint arXiv:171105225. 2017.
  * 8 Radford A, Kim JW, Hallacy C, Ramesh A, Goh G, Agarwal S, et al.  Learning transferable visual models from natural language supervision.  In: International conference on machine learning. PMLR; 2021. p. 8748-63.
  * 9 Tsai YHH, Bai S, Liang PP, Kolter JZ, Morency LP, Salakhutdinov R.  Multimodal transformer for unaligned multimodal language sequences.  In: Proceedings of the conference. Association for Computational Linguistics. Meeting. vol. 2019. NIH Public Access; 2019. p. 6558.
  * 10 Chen YC, Li L, Yu L, El Kholy A, Ahmed F, Gan Z, et al.  Uniter: Universal image-text representation learning.  In: Computer Vision–ECCV 2020: 16th European Conference, Glasgow, UK, August 23–28, 2020, Proceedings, Part XXX. Springer; 2020. p. 104-20.
  * 11 Kim W, Son B, Kim I.  Vilt: Vision-and-language transformer without convolution or region supervision.  In: International Conference on Machine Learning. PMLR; 2021. p. 5583-94.
  * 12 Yan B, Pei M.  Clinical-BERT: Vision-Language Pre-training for Radiograph Diagnosis and Reports Generation.  In: Proceedings of the AAAI Conference on Artificial Intelligence. vol. 36; 2022. p. 2982-90.
  * 13 Moon JH, Lee H, Shin W, Kim YH, Choi E.  Multi-modal understanding and generation for medical images and text via vision-language pre-training.  IEEE Journal of Biomedical and Health Informatics. 2022;26(12):6070-80.
  * 14 Huang SC, Shen L, Lungren MP, Yeung S.  Gloria: A multimodal global-local representation learning framework for label-efficient medical image recognition.  In: Proceedings of the IEEE/CVF International Conference on Computer Vision; 2021. p. 3942-51.
  * 15 Zhang Y, Jiang H, Miura Y, Manning CD, Langlotz CP.  Contrastive learning of medical visual representations from paired images and text.  In: Machine Learning for Healthcare Conference. PMLR; 2022. p. 2-25.
  * 16 Tiu E, Talius E, Patel P, Langlotz CP, Ng AY, Rajpurkar P.  Expert-level detection of pathologies from unannotated chest X-ray images via self-supervised learning.  Nature Biomedical Engineering. 2022:1-8.
  * 17 Xian Y, Lampert CH, Schiele B, Akata Z.  Zero-shot learning—a comprehensive evaluation of the good, the bad and the ugly.  IEEE transactions on pattern analysis and machine intelligence. 2018;41(9):2251-65.
  * 18 Chauhan G, Liao R, Wells W, Andreas J, Wang X, Berkowitz S, et al.  Joint modeling of chest radiographs and radiology reports for pulmonary edema assessment.  In: Medical Image Computing and Computer Assisted Intervention–MICCAI 2020: 23rd International Conference, Lima, Peru, October 4–8, 2020, Proceedings, Part II 23. Springer; 2020. p. 529-39.
  * 19 Chen T, Kornblith S, Norouzi M, Hinton G.  A simple framework for contrastive learning of visual representations.  In: International conference on machine learning. PMLR; 2020. p. 1597-607.
  * 20 He J, Chen JN, Liu S, Kortylewski A, Yang C, Bai Y, et al.  Transfg: A transformer architecture for fine-grained recognition.  In: Proceedings of the AAAI Conference on Artificial Intelligence. vol. 36; 2022. p. 852-60.
  * 21 Dosovitskiy A, Beyer L, Kolesnikov A, Weissenborn D, Zhai X, Unterthiner T, et al.  An image is worth 16x16 words: Transformers for image recognition at scale.  arXiv preprint arXiv:201011929. 2020.
  * 22 Li K, Zhang Y, Li K, Li Y, Fu Y.  Visual semantic reasoning for image-text matching.  In: Proceedings of the IEEE/CVF International conference on computer vision; 2019. p. 4654-62.
  * 23 Chan HP, Samala RK, Hadjiiski LM, Zhou C.  Deep learning in medical image analysis.  Deep Learning in Medical Image Analysis: Challenges and Applications. 2020:3-21.
  * 24 Li Q, Nishikawa RM.  Computer-aided detection and diagnosis in medical imaging.  Taylor & Francis; 2015.
  * 25 Krizhevsky A, Sutskever I, Hinton GE.  Imagenet classification with deep convolutional neural networks.  Communications of the ACM. 2017;60(6):84-90.
  * 26 Vaswani A, Shazeer N, Parmar N, Uszkoreit J, Jones L, Gomez AN, et al.  Attention is all you need.  Advances in neural information processing systems. 2017;30.
  * 27 Sowrirajan H, Yang J, Ng AY, Rajpurkar P.  Moco pretraining improves representation and transferability of chest x-ray models.  In: Medical Imaging with Deep Learning. PMLR; 2021. p. 728-44.
  * 28 Cao H, Wang Y, Chen J, Jiang D, Zhang X, Tian Q, et al.  Swin-unet: Unet-like pure transformer for medical image segmentation.  In: Computer Vision–ECCV 2022 Workshops: Tel Aviv, Israel, October 23–27, 2022, Proceedings, Part III. Springer; 2023. p. 205-18.
  * 29 Ker J, Wang L, Rao J, Lim T.  Deep learning applications in medical image analysis.  Ieee Access. 2017;6:9375-89.
  * 30 Wang S, Summers RM.  Machine learning and radiology.  Medical image analysis. 2012;16(5):933-51.
  * 31 Karpathy A, Fei-Fei L.  Deep visual-semantic alignments for generating image descriptions.  In: Proceedings of the IEEE conference on computer vision and pattern recognition; 2015. p. 3128-37.
  * 32 Socher R, Ganjoo M, Manning CD, Ng A.  Zero-shot learning through cross-modal transfer.  Advances in neural information processing systems. 2013;26.
  * 33 Johnson AE, Pollard TJ, Berkowitz SJ, Greenbaum NR, Lungren MP, Deng Cy, et al.  MIMIC-CXR, a de-identified publicly available database of chest radiographs with free-text reports.  Scientific data. 2019;6(1):317.
  * 34 Irvin J, Rajpurkar P, Ko M, Yu Y, Ciurea-Ilcus S, Chute C, et al.  Chexpert: A large chest radiograph dataset with uncertainty labels and expert comparison.  In: Proceedings of the AAAI conference on artificial intelligence. vol. 33; 2019. p. 590-7.
  * 35 Russakovsky O, Deng J, Su H, Krause J, Satheesh S, Ma S, et al.  Imagenet large scale visual recognition challenge.  International journal of computer vision. 2015;115:211-52.
  * 36 Van Horn G, Mac Aodha O, Song Y, Cui Y, Sun C, Shepard A, et al.  The inaturalist species classification and detection dataset.  In: Proceedings of the IEEE conference on computer vision and pattern recognition; 2018. p. 8769-78.
  * 37 Krause J, Stark M, Deng J, Fei-Fei L.  3d object representations for fine-grained categorization.  In: Proceedings of the IEEE international conference on computer vision workshops; 2013. p. 554-61.
  * 38 Jain S, Wallace BC.  Attention is not explanation.  arXiv preprint arXiv:190210186. 2019.
  * 39 Nazal Y, et al.  Rib fracture: Different radiographic projections.  Medical Science Monitor. 2012;77(4):13-6.
  * 40 Wang Z, Wu Z, Agarwal D, Sun J.  MedCLIP: Contrastive Learning from Unpaired Medical Images and Text.  In: Proceedings of the 2022 Conference on Empirical Methods in Natural Language Processing. Abu Dhabi, United Arab Emirates: Association for Computational Linguistics; 2022. p. 3876-87.  Available from: https://aclanthology.org/2022.emnlp-main.256.
  * 41 Seibold C, Reiß S, Sarfraz MS, Stiefelhagen R, Kleesiek J.  Breaking with fixed set pathology recognition through report-guided contrastive training.  In: Medical Image Computing and Computer Assisted Intervention–MICCAI 2022: 25th International Conference, Singapore, September 18–22, 2022, Proceedings, Part V. Springer; 2022. p. 690-700.
  * 42 Müller P, Kaissis G, Zou C, Rueckert D.  Joint learning of localized representations from medical images and reports.  In: Computer Vision–ECCV 2022: 17th European Conference, Tel Aviv, Israel, October 23–27, 2022, Proceedings, Part XXVI. Springer; 2022. p. 685-701.
  * 43 Boecking B, Usuyama N, Bannur S, Castro DC, Schwaighofer A, Hyland S, et al.  Making the most of text semantics to improve biomedical vision–language processing.  In: Computer Vision–ECCV 2022: 17th European Conference, Tel Aviv, Israel, October 23–27, 2022, Proceedings, Part XXXVI. Springer; 2022. p. 1-21.
  * 44 Wang F, Zhou Y, Wang S, Vardhanabhuti V, Yu L.  Multi-Granularity Cross-modal Alignment for Generalized Medical Visual Representation Learning.  In: Oh AH, Agarwal A, Belgrave D, Cho K, editors. Advances in Neural Information Processing Systems; 2022. Available from: https://openreview.net/forum?id=Yul402KcD5d.
  * 45 Chaitanya K, Erdil E, Karani N, Konukoglu E.  Contrastive learning of global and local features for medical image segmentation with limited annotations.  Advances in Neural Information Processing Systems. 2020;33:12546-58.


# Evidence that Ultra-High-Energy Gamma Rays are a Universal Feature Near
Powerful Pulsars

A. Albert  Physics Division, Los Alamos National Laboratory, Los Alamos, NM,
USA  R. Alfaro  Instituto de Física, Universidad Nacional Autónoma de México,
Ciudad de Mexico, Mexico  C. Alvarez  Universidad Autónoma de Chiapas, Tuxtla
Gutiérrez, Chiapas, México  J.D. Álvarez  Universidad Michoacana de San
Nicolás de Hidalgo, Morelia, Mexico  J.R. Angeles Camacho  Instituto de
Física, Universidad Nacional Autónoma de México, Ciudad de Mexico, Mexico
J.C. Arteaga-Velázquez  Universidad Michoacana de San Nicolás de Hidalgo,
Morelia, Mexico  K.P. Arunbabu  Instituto de Geofísica, Universidad Nacional
Autónoma de México, Ciudad de Mexico, Mexico  D. Avila Rojas  Instituto de
Física, Universidad Nacional Autónoma de México, Ciudad de Mexico, Mexico
H.A. Ayala Solares  Department of Physics, Pennsylvania State University,
University Park, PA, USA  V. Baghmanyan  Institute of Nuclear Physics Polish
Academy of Sciences, PL-31342 IFJ-PAN, Krakow, Poland  E. Belmont-Moreno
Instituto de Física, Universidad Nacional Autónoma de México, Ciudad de
Mexico, Mexico  S.Y. BenZvi  Department of Physics & Astronomy, University of
Rochester, Rochester, NY , USA  C. Brisbois  Department of Physics, University
of Maryland, College Park, MD, USA  K.S. Caballero-Mora  Universidad Autónoma
de Chiapas, Tuxtla Gutiérrez, Chiapas, México  T. Capistrán  Instituto de
Astronomía, Universidad Nacional Autónoma de México, Ciudad de Mexico, Mexico
A. Carramiñana  Instituto Nacional de Astrofísica, Óptica y Electrónica,
Puebla, Mexico  S. Casanova  Institute of Nuclear Physics Polish Academy of
Sciences, PL-31342 IFJ-PAN, Krakow, Poland  U. Cotti  Universidad Michoacana
de San Nicolás de Hidalgo, Morelia, Mexico  J. Cotzomi  Facultad de Ciencias
Físico Matemáticas, Benemérita Universidad Autónoma de Puebla, Puebla, Mexico
S. Coutiño de León  Instituto Nacional de Astrofísica, Óptica y Electrónica,
Puebla, Mexico  E. De la Fuente  Departamento de Física, Centro Universitario
de Ciencias Exactase Ingenierias, Universidad de Guadalajara, Guadalajara,
Mexico  C. de León  Universidad Michoacana de San Nicolás de Hidalgo, Morelia,
Mexico  R. Diaz Hernandez  Instituto Nacional de Astrofísica, Óptica y
Electrónica, Puebla, Mexico  B.L. Dingus  Physics Division, Los Alamos
National Laboratory, Los Alamos, NM, USA  M.A. DuVernois  Department of
Physics, University of Wisconsin-Madison, Madison, WI, USA  M. Durocher
Physics Division, Los Alamos National Laboratory, Los Alamos, NM, USA  J.C.
Díaz-Vélez  Departamento de Física, Centro Universitario de Ciencias Exactase
Ingenierias, Universidad de Guadalajara, Guadalajara, Mexico  R.W. Ellsworth
Department of Physics, University of Maryland, College Park, MD, USA  K. Engel
Department of Physics, University of Maryland, College Park, MD, USA  C.
Espinoza  Instituto de Física, Universidad Nacional Autónoma de México, Ciudad
de Mexico, Mexico  K.L. Fan  Department of Physics, University of Maryland,
College Park, MD, USA  M. Fernández Alonso  Department of Physics,
Pennsylvania State University, University Park, PA, USA  N. Fraija  Instituto
de Astronomía, Universidad Nacional Autónoma de México, Ciudad de Mexico,
Mexico  A. Galván-Gámez  Instituto de Astronomía, Universidad Nacional
Autónoma de México, Ciudad de Mexico, Mexico  J.A. García-González
Tecnologico de Monterrey, Escuela de Ingeniería y Ciencias, Ave. Eugenio Garza
Sada 2501, Monterrey, N.L., Mexico, 64849  F. Garfias  Instituto de
Astronomía, Universidad Nacional Autónoma de México, Ciudad de Mexico, Mexico
G. Giacinti  Max-Planck Institute for Nuclear Physics, 69117 Heidelberg,
Germany  M.M. González  Instituto de Astronomía, Universidad Nacional Autónoma
de México, Ciudad de Mexico, Mexico  J.A. Goodman  Department of Physics,
University of Maryland, College Park, MD, USA  J.P. Harding  Physics Division,
Los Alamos National Laboratory, Los Alamos, NM, USA  S. Hernandez  Instituto
de Física, Universidad Nacional Autónoma de México, Ciudad de Mexico, Mexico
B. Hona  Department of Physics and Astronomy, University of Utah, Salt Lake
City, UT, USA  D. Huang  Department of Physics, Michigan Technological
University, Houghton, MI, USA  F. Hueyotl-Zahuantitla  Universidad Autónoma de
Chiapas, Tuxtla Gutiérrez, Chiapas, México  P. Hüntemeyer  Department of
Physics, Michigan Technological University, Houghton, MI, USA  A. Iriarte
Instituto de Astronomía, Universidad Nacional Autónoma de México, Ciudad de
Mexico, Mexico  A. Jardin-Blicq  Max-Planck Institute for Nuclear Physics,
69117 Heidelberg, Germany  Department of Physics, Faculty of Science,
Chulalongkorn University, 254 Phayathai Road, Pathumwan, Bangkok 10330,
Thailand  National Astronomical Research Institute of Thailand (Public
Organization), Don Kaeo, MaeRim, Chiang Mai 50180, Thailand  V. Joshi
Erlangen Centre for Astroparticle Physics, Friedrich-Alexander-Universität
Erlangen-Nürnberg, Erlangen, Germany  D. Kieda  Department of Physics and
Astronomy, University of Utah, Salt Lake City, UT, USA  A. Lara  Instituto de
Geofísica, Universidad Nacional Autónoma de México, Ciudad de Mexico, Mexico
W.H. Lee  Instituto de Astronomía, Universidad Nacional Autónoma de México,
Ciudad de Mexico, Mexico  J. Lee  University of Seoul, Seoul, Rep. of Korea
H. León Vargas  Instituto de Física, Universidad Nacional Autónoma de México,
Ciudad de Mexico, Mexico  J.T. Linnemann  Department of Physics and Astronomy,
Michigan State University, East Lansing, MI, USA  A.L. Longinotti  Instituto
Nacional de Astrofísica, Óptica y Electrónica, Puebla, Mexico  Instituto de
Astronomía, Universidad Nacional Autónoma de México, Ciudad de Mexico, Mexico
G. Luis-Raya  Universidad Politecnica de Pachuca, Pachuca, Hgo, Mexico  J.
Lundeen  Department of Physics and Astronomy, Michigan State University, East
Lansing, MI, USA  K. Malone  Physics Division, Los Alamos National Laboratory,
Los Alamos, NM, USA  V. Marandon  Max-Planck Institute for Nuclear Physics,
69117 Heidelberg, Germany  O. Martinez  Facultad de Ciencias Físico
Matemáticas, Benemérita Universidad Autónoma de Puebla, Puebla, Mexico  J.
Martínez-Castro  Centro de Investigación en Computación, Instituto Politécnico
Nacional, México City, México.  J.A. Matthews  Dept of Physics and Astronomy,
University of New Mexico, Albuquerque, NM, USA  P. Miranda-Romagnoli
Universidad Autónoma del Estado de Hidalgo, Pachuca, Mexico  J.A. Morales-Soto
Universidad Michoacana de San Nicolás de Hidalgo, Morelia, Mexico  E. Moreno
Facultad de Ciencias Físico Matemáticas, Benemérita Universidad Autónoma de
Puebla, Puebla, Mexico  M. Mostafá  Department of Physics, Pennsylvania State
University, University Park, PA, USA  A. Nayerhoda  Institute of Nuclear
Physics Polish Academy of Sciences, PL-31342 IFJ-PAN, Krakow, Poland  L.
Nellen  Instituto de Ciencias Nucleares, Universidad Nacional Autónoma de
Mexico, Ciudad de Mexico, Mexico  M. Newbold  Department of Physics and
Astronomy, University of Utah, Salt Lake City, UT, USA  M.U. Nisa  Department
of Physics and Astronomy, Michigan State University, East Lansing, MI, USA  R.
Noriega-Papaqui  Universidad Autónoma del Estado de Hidalgo, Pachuca, Mexico
L. Olivera-Nieto  Max-Planck Institute for Nuclear Physics, 69117 Heidelberg,
Germany  N. Omodei  Department of Physics, Stanford University: Stanford, CA
94305–4060, USA  A. Peisker  Department of Physics and Astronomy, Michigan
State University, East Lansing, MI, USA  Y. Pérez Araujo  Instituto de
Astronomía, Universidad Nacional Autónoma de México, Ciudad de Mexico, Mexico
E.G. Pérez-Pérez  Universidad Politecnica de Pachuca, Pachuca, Hgo, Mexico
C.D. Rho  University of Seoul, Seoul, Rep. of Korea  Y.J. Roh  University of
Seoul, Seoul, Rep. of Korea  D. Rosa-González  Instituto Nacional de
Astrofísica, Óptica y Electrónica, Puebla, Mexico  E. Ruiz-Velasco  Max-Planck
Institute for Nuclear Physics, 69117 Heidelberg, Germany  H. Salazar  Facultad
de Ciencias Físico Matemáticas, Benemérita Universidad Autónoma de Puebla,
Puebla, Mexico  F. Salesa Greus  Institute of Nuclear Physics Polish Academy
of Sciences, PL-31342 IFJ-PAN, Krakow, Poland  Instituto de Física
Corpuscular, CSIC, Universitat de València, E-46980, Paterna, Valencia, Spain
A. Sandoval  Instituto de Física, Universidad Nacional Autónoma de México,
Ciudad de Mexico, Mexico  M. Schneider  Department of Physics, University of
Maryland, College Park, MD, USA  H. Schoorlemmer  Max-Planck Institute for
Nuclear Physics, 69117 Heidelberg, Germany  J. Serna-Franco  Instituto de
Física, Universidad Nacional Autónoma de México, Ciudad de Mexico, Mexico
A.J. Smith  Department of Physics, University of Maryland, College Park, MD,
USA  R.W. Springer  Department of Physics and Astronomy, University of Utah,
Salt Lake City, UT, USA  P. Surajbali  Max-Planck Institute for Nuclear
Physics, 69117 Heidelberg, Germany  M. Tanner  Department of Physics,
Pennsylvania State University, University Park, PA, USA  K. Tollefson
Department of Physics and Astronomy, Michigan State University, East Lansing,
MI, USA  I. Torres  Instituto Nacional de Astrofísica, Óptica y Electrónica,
Puebla, Mexico  R. Torres-Escobedo  Departamento de Física, Centro
Universitario de Ciencias Exactase Ingenierias, Universidad de Guadalajara,
Guadalajara, Mexico  Tsung-Dao Lee Institute $\&$ School of Physics and
Astronomy, Shanghai Jiao Tong University, Shanghai, China  R. Turner
Department of Physics, Michigan Technological University, Houghton, MI, USA
F. Ureña-Mena  Instituto Nacional de Astrofísica, Óptica y Electrónica,
Puebla, Mexico  L. Villaseñor  Facultad de Ciencias Físico Matemáticas,
Benemérita Universidad Autónoma de Puebla, Puebla, Mexico  T. Weisgarber
Department of Physics, University of Wisconsin-Madison, Madison, WI, USA  E.
Willox  Department of Physics, University of Maryland, College Park, MD, USA
H. Zhou  Tsung-Dao Lee Institute $\&$ School of Physics and Astronomy,
Shanghai Jiao Tong University, Shanghai, China missing<EMAIL_ADDRESS>
###### Abstract

The highest-energy known gamma-ray sources are all located within 0.5 degrees
of extremely powerful pulsars. This raises the question of whether ultra-high-
energy (UHE; $>$ 56 TeV) gamma-ray emission is a universal feature expected
near pulsars with a high spin-down power. Using four years of data from the
High Altitude Water Cherenkov (HAWC) Gamma-Ray Observatory, we present a
joint-likelihood analysis of ten extremely powerful pulsars to search for
subthreshold UHE gamma-ray emission correlated with these locations. We report
a significant detection ($>$ 3$\sigma$), indicating that UHE gamma-ray
emission is a generic feature of powerful pulsars. We discuss the emission
mechanisms of the gamma rays and the implications of this result. The
individual environment, such as the magnetic field and particle density in the
surrounding area, appears to play a role in the amount of emission.

††journal: The Astrophysical Journal Letters

=500

Kelly Malone

##  1 Introduction

Ultra-high-energy (UHE; $>$ 56 TeV) gamma-ray emission can be created via
hadronic or leptonic processes. In the hadronic mechanism, a neutral pion
decays into two gamma rays. In the leptonic mechanism, a lower-energy photon
scatters off a relativistic electron via inverse Compton scattering. The
electron transfers part of its energy to the gamma ray, resulting in a higher-
energy photon.

Traditionally, it was thought that UHE gamma-ray sources would be hadronic in
nature, as leptonic emission is suppressed in this energy regime due to the
Klein-Nishina (KN) effect. However, present-day gamma-ray observatories have
the sensitivity required to detect leptonic sources above 56 TeV. See, for
example, the detections of the Crab Nebula as well as several known pulsar
wind nebulae (PWN) (Abeysekara et al., 2019; Amenomori et al., 2019;
Abeysekara et al., 2020; Abdalla et al., 2019). The astrophysical spectrum of
a leptonic source typically exhibits significant curvature due to the KN
effect (Moderski et al., 2005). Conversely, hadronic sources follow the
spectrum of their parent cosmic-ray population, which may or may not include a
cutoff or curvature.

The High Altitude Water Cherenkov (HAWC) Observatory is a gamma-ray
observatory with an wide instantaneous field-of-view ($\sim$2 steradians) and
sensitivity to energies between a few hundred GeV and a few hundred TeV. It is
sensitive to sources with declinations between -26 and +64 degrees (Smith,
2015; Abeysekara et al., 2017, 2019).

The first HAWC catalog of UHE sources (Abeysekara et al. (2020); hereafter
referred to as the “eHWC” catalog) contains nine sources emitting above 56
TeV, three of which continue above 100 TeV. The highest-energy sources all
exhibit curvature in the spectrum. Additionally, all nine sources are located
within 0.5 degrees of pulsars. For eight of the nine sources, the pulsar has
an extremely high spin-down power ($\dot{E}>10^{36}$ erg/s). This is much
higher than the number of high-$\dot{E}$ pulsars that would be expected to be
found near UHE gamma-ray sources (Abeysekara et al., 2020). Emission near
pulsars is expected to be powered by a PWN or TeV halo and is therefore
dominantly leptonic, even at the high energies studied here (Breuhaus et al.,
2020; Sudoh et al., 2021). While young pulsars are expected to have an
associated supernova remnant (SNR) with accompanying hadronic emission, there
has only been one detection of gamma-ray emission from an SNR to UHE, and
leptonic emission from this source has not been conclusively ruled out (Albert
et al., 2020a). SNR theory starts to run into technical problems accelerating
particles to these energies (Gabici, 2017). SNR acceleration to UHE energies
may not be possible (Zeng et al., 2021).

The proximity of these gamma-ray sources to the most powerful pulsars, along
with the curvature in their spectra, invites the question of whether UHE
gamma-ray emission is a generic feature expected near these sources. This is
investigated in this paper through a joint-likelihood analysis of pulsars with
$\dot{E}>10^{36}$ erg/s to search for subthreshold sources in the HAWC data.
While each source is too weak to be individually detected in HAWC’s standard
catalog search, analyzing the regions jointly may lead to a general detection
for this source class.

The paper is organized as follows: Section 2 describes the details of the
joint-likelihood method. Section 3 contains the results of the analysis.
Section 4 discusses implications of the results.

##  2 Analysis method

###  2.1 Joint-Likelihood Method

In this paper, we perform a joint analysis of several high-$\dot{E}$ pulsars.
The analysis uses a binned maximum-likelihood method that searches for a
gamma-ray excess above the background, using a simple model that describes the
UHE emission from pulsars based on various observables such as the distance
and pulsar age. This is performed using the HAWC Accelerated Likelihood
(HAL)111https://github.com/threeML/hawc_hal plugin to the Multi-Mission
Maximum Likelihood Framework (3ML) (Vianello et al., 2015). We describe the
background rejection and likelihood method in Abeysekara et al. (2019). HAWC’s
background comes from two main sources: the Galactic diffuse background
stemming from unresolved gamma-ray sources, and cosmic rays that are detected
by HAWC and survive gamma/hadron separation cuts. Note that neither background
is very large above 56 TeV. The Galactic diffuse background is largely due to
Inverse Compton scattering, which is not very prominent at these energies. The
fraction of cosmic rays surviving gamma/hadron separation cuts is also very
low: $\sim$0.001 (Abeysekara et al., 2017).

Data were collected over 1343 days between June 2015 and June 2019. The data
are binned in quarter-decade bins of reconstructed gamma-ray energy using the
“ground parameter” energy estimator, one of two energy estimators currently
used by HAWC (Abeysekara et al., 2019).

We use the last three quarter-decade energy bins from Abeysekara et al.
(2019), restricting the analysis to reconstructed energies between 56 and 316
TeV. These are the highest energies probed by HAWC. At these energies, there
is very little diffuse emission and less source confusion than at lower
energies.

The spectral model is a power law:

$\frac{dN}{dE}=A_{i}K\left(\frac{E}{E_{0}}\right)^{-\alpha},$ (1)

where $A_{i}$ accounts for the model-dependent relative flux of each source
(see Section 2.3), $K$ is the normalization, and $E_{0}$ is the pivot energy,
which is fixed at the center of each energy bin. The three values used for
$E_{0}$ in this analysis are 74.99 TeV, 133.35 TeV, and 237.14 TeV. Setting
the pivot energy at the center of each bin makes the analysis relatively
insensitive to the choice of $\alpha$, which is fixed at 2.5.

We fit each source (see Table 1) individually in each energy bin, placing a
step function at the boundaries of the bin to ensure that events that are mis-
reconstructed in energy are excluded. We assume the sources are spatially
extended with a Gaussian morphology. The width is fixed at $\sigma$=0.34∘. The
values of $\alpha$ and $\sigma$ are approximately the average values for the
highest-energy gamma-ray sources from the eHWC catalog (Abeysekara et al.,
2020).

To obtain a joint-likelihood result, the log-likelihood profiles for each
individual source, in each energy bin, are added and the value of $K$
(hereafter called $\hat{K}$) that optimizes this log-likelihood profile is
found. The total flux normalization, $\kappa$, for the ten sources combined is
then just:

$\kappa=\sum_{i=1}^{n_{sources}}\hat{K}A_{i}.$ (2)

Using this flux normalization, the total differential flux from all ten
sources is:

$\frac{dN}{dE}=\kappa\left(\frac{E}{E_{0}}\right)^{-\alpha}.$ (3)

We calculate a test statistic (TS) to show how significant the value of
$\kappa$ is:

$\mathrm{TS}=2\mathrm{ln}\frac{L_{S+B}(\kappa)}{L_{B}},$ (4)

where $L_{S+B}$ is the maximum likelihood for the signal-plus-background
hypothesis and $L_{B}$ is the likelihood for the background-only hypothesis.

After an overall best-fit value of $\kappa$ is determined, it is used as an
input to a Markov chain Monte Carlo and a distribution of the value is
determined. If the TS in a given energy bin is significant (TS $>4$), a
Bayesian credible interval (68$\%$ containment) is shown, obtained from the
estimated distribution of the parameter $\kappa$. The model only has one
degree of freedom, so using Wilks Theorem (Wilks, 1938), a TS value of 4
corresponds to 2$\sigma$. Otherwise, 95$\%$ credible interval quasi-
differential upper limit is plotted. A uniform prior is used in the Bayesian
analysis.

In the bins where an upper limit is determined, the range of expected upper
limits are obtained from simulations of Poisson-fluctuated background.

###  2.2 Source Selection and Dataset

Figure 1: $\sqrt{TS}$ for $\hat{E}>$ 56 TeV emission from HAWC, for the field-
of-view used in this analysis. We assume a point source morphology. The
locations of the ten sub-threshold pulsars used in this analysis are labeled.

For this analysis, we define the source list by selecting all pulsars from the
ATNF pulsar database, version
1.62222https://www.atnf.csiro.au/research/pulsar/psrcat/ (Manchester et al.,
2005) in the inner Galactic plane that are within HAWC’s field-of-view ($|b|<$
1∘, 5∘ $<l<$ 90∘) and have $\dot{E}>$ 1036 erg/s. There are 24 pulsars that
meet this criteria (see Table 1). Pulsars with a high $\dot{E}$ are centered
around the Galactic plane, so the choice to concentrate on $|b|<1^{\circ}$
only removes 3 additional pulsars from the analysis. We use this list of 24
pulsars to determine which models of emission are reasonable (see Section
2.3).

After making this determination, this list is then downselected to search for
sub-threshold gamma-ray sources. First, we remove all pulsars that are located
within a degree of sources from the eHWC catalog (Abeysekara et al., 2020).
This removes pulsars that already have significant UHE emission detected in
their vicinity. Since the gamma-ray emission is modeled as extended in nature,
this also removes pulsars whose associated emission may overlap with those
known sources, which would require more detailed modeling.

We remove three additional pulsars from the source selection. PSR J1813-1749
is removed because HAWC has added more data since the publication of
Abeysekara et al. (2020) and has detected a new UHE source (eHWC-J1813-176)
$\sim$0.2∘ away from the pulsar. PSR J1913+1011 and PSR J1930+1852 have been
removed because the known TeV emission in those regions is likely from a SNR,
not a PWN (Abdalla et al., 2018a, b). Since the majority of the emission at
these energies is expected to be from a PWN or TeV halo, this is done to
prevent the introduction of a different, predominantly hadronic source class
into the analysis.

The final list is composed of ten pulsars that are candidates for sub-
threshold gamma-ray emission. Figure 1 shows HAWC’s $>$ 56 TeV map with these
sources labeled.

PSR name | RA (∘) | Dec (∘) | Age ($\frac{P}{2\dot{P}}$) (kyr) |  $\dot{E}$ ($\times$1036 erg/s) | Distance (kpc) | P (s) |  $\dot{P}$ (ss-1) | Subthreshold?
---|---|---|---|---|---|---|---|---
B1800-21 | 270.96 | -21.62 | 15.8 | 2.2 | 4.40 | 0.134 | 1.35 $\times 10^{-13}$ | $\checkmark$
J1809-1917 | 272.43 | -19.29 | 51.3 | 1.8 | 3.27 | 0.083 | 2.55 $\times 10^{-14}$ |
J1811-1925 | 272.87 | -19.42 | 23.3 | 6.4 | 5.00 | 0.065 | 4.4 $\times 10^{-14}$ |
J1813-1749 | 273.40 | -17.83 | 5.6 | 56 | 4.70 | 0.045 | 1.27 $\times 10^{-13}$ |
J1826-1256 | 276.54 | -12.94 | 14.4 | 3.6 | 1.55 | 0.110 | 1.21 $\times 10^{-13}$ |
B1823-13 | 276.55 | -13.58 | 21.4 | 2.8 | 3.61 | 0.101 | 7.53 $\times 10^{-14}$ |
J1828-1101 | 277.08 | -11.03 | 77.1 | 1.6 | 4.77 | 0.072 | 1.48 $\times 10^{-14}$ | $\checkmark$
J1831-0952 | 277.89 | -9.87 | 128 | 1.1 | 3.68 | 0.067 | 8.32 $\times 10^{-15}$ | $\checkmark$
J1833-1034 | 278.39 | -10.57 | 4.85 | 34 | 4.10 | 0.062 | 2.02 $\times 10^{-13}$ | $\checkmark$
J1837-0604 | 279.43 | -6.08 | 33.8 | 2.0 | 4.77 | 0.096 | 45.1 $\times 10^{-15}$ |
J1838-0537 | 279.73 | -5.62 | 4.89 | 6.0 | 2.0333Pseudo-distance from Pletsch et al. (2012) | 0.146 | 4.72 $\times 10^{-13}$ |
J1838-0655 | 279.51 | -6.93 | 22.7 | 5.5 | 6.60 | 0.070 | 4.93 $\times 10^{-14}$ | $\checkmark$
J1844-0346 | 281.14 | -3.78 | 11.6 | 4.2 | 2.40444Pseudo-distance derived from Eq. 3 of Wu et al. (2018) | 0.113 | 1.55 $\times 10^{-13}$ |
J1846-0258 | 281.60 | -2.98 | 0.728 | 8.1 | 5.8 | 0.327 | 7.11 $\times 10^{-12}$ |
J1849-0001 | 282.23 | -0.02 | 42.9 | 9.8 | 7.0555Distance estimate from Gotthelf et al. (2011) | 0.039 | 1.42 $\times 10^{-14}$ |
J1856+0245 | 284.21 | 2.76 | 20.6 | 4.6 | 6.32 | 0.081 | 6.21 $\times 10^{-14}$ | $\checkmark$
J1907+0602 | 286.98 | 6.04 | 19.5 | 2.8 | 2.37 | 0.107 | 8.68 $\times 10^{-14}$ |
J1913+1011 | 288.33 | 10.19 | 169 | 2.9 | 4.61 | 0.036 | 3.37 $\times 10^{-15}$ |
J1928+1746 | 292.18 | 17.77 | 82.6 | 1.6 | 4.34 | 0.069 | 1.32 $\times 10^{-14}$ | $\checkmark$
J1930+1852 | 292.63 | 18.87 | 2.89 | 12 | 7.00 | 0.137 | 7.51 $\times 10^{-13}$ |
J1935+2025 | 293.92 | 20.43 | 20.9 | 4.7 | 4.60 | 0.080 | 6.08 $\times 10^{-14}$ | $\checkmark$
B1937+21 | 294.91 | 21.58 | 2.35e5 | 1.1 | 3.50 | 0.002 | 1.05$\times 10^{-19}$ | $\checkmark$
J2021+3651 | 305.27 | 36.85 | 17.2 | 3.4 | 1.8 | 0.104 | 9.57 $\times 10^{-14}$ |
J2022+3842 | 305.59 | 38.70 | 8.94 | 30 | 10.00 | 0.049 | 8.61 $\times 10^{-14}$ | $\checkmark$
Table 1: Information on the pulsars. All information comes from the ATNF
database, version 1.62 (Manchester et al., 2005) except for some distance
estimates which are not included in the pulsar database (see footnotes). Here,
$\dot{E}$ is the spin-down energy loss rate, $P$ is the barycentric period,
and $\dot{P}$ is the time derivative of the period. The checkmark in the last
column denotes the ten pulsars included in the sub-threshold analysis.

###  2.3 Models

The parameter $A_{i}$ in Equation 1 describes the relative contribution each
pulsar receives in the analysis. This parameter can be used to test different
models of gamma-ray emission near pulsars. For the models that rely on pulsar
parameters, the relevant quantities are taken from the ATNF pulsar catalog,
version 1.62 (Manchester et al., 2005). We consider a variety of different
models. In all descriptions, “emission” refers solely to gamma-ray emission
above 56 TeV.

  * •

No model: Here, $A_{i}$ for all sources is set equal to 1. All sources are
treated equally and the emission is expected to be uncorrelated with pulsar
parameters such as distance.

  * •

$\boldsymbol{1/d^{2}}$: In this model, $A_{i}$ is set to 1/$d^{2}$, where $d$
is the distance to the pulsar. This model assumes that closer sources will
produce observable emission.

  * •

$\boldsymbol{\dot{E}/d^{2}}$: Here, the 1/$d^{2}$ model discussed above is
multiplied by the spin-down power of the pulsar. Therefore, closer, more-
energetic pulsars have more gamma-ray emission.

  * •

Inverse age: In this model, $A_{i}$ is the inverse of the spin-down age. This
is defined as $P/(2\dot{P})$, where $P$ and $\dot{P}$ are the period and the
time derivative of the period, respectively. In this model, younger sources
are more likely to have detectable emission.

  * •

Flux at 7 TeV: Here, $A_{i}$ is the HAWC flux at 7 TeV. This model assumes
that sources that are bright at multi-TeV energies should also give off
detectable emission above 56 TeV. The 0.5∘ extended source map from the third
catalog of HAWC sources (3HWC) (Albert et al., 2020b) is used to extract these
values. $A_{i}$ is computed by averaging the flux from all pixels within a 0.5
degree radius of the pulsar.

##  3 Results

We first run the joint-likelihood analysis (as described in Section 2) with
the full list of 24 pulsars to determine which models of emission are
reasonable. Table 2 gives the total TS for each scenario. For each model, we
can also calculate the expected gamma-ray flux above 56 TeV from an arbitrary
pulsar using Equation 1. The values of $K$, derived from Equation 2 are also
given in Table 2. In all cases, the TS is much higher than 5$\sigma$. Two
scenarios, the flux at 7 TeV and 1/$d^{2}$, perform better than the “no model
case”. The model based on the inverse age of the pulsar performs significantly
worse than the others, but the significance, $\sqrt{TS}$ is still well above
the 5$\sigma$ level.

Model |  $K$ (56 $<$ E $<$ 100 TeV) |  $K$ (100 $<$ E $<$ 177 TeV) |  $K$ (177 $<$ E $<$ 316 TeV) | Units for $A_{i}$ | Total TS
---|---|---|---|---|---
No model | 2.55 $\times$ 10-16 | 4.65 $\times$ 10-17 | 9.46 $\times$ 10-18 | dimensionless | 633
1/$d^{2}$ | 2.44 $\times$ 10-16 | 4.50 $\times$ 10-16 | 8.12 $\times$ 10-17 | kpc-2 | 734
$\dot{E}/d^{2}$ | 4.06 $\times$ 10-13 | 6.88 $\times$ 10-14 | 1.04 $\times$ 10-14 | erg m-2 s-1 | 568
Inverse age | 6.36 $\times$ 10-13 | 1.10 $\times$ 10-13 | 1.65 $\times$ 10-14 | yr-1 | 152
Flux at 7 TeV | 5.42 $\times$ 10-3 | 9.11 $\times$ 10-4 | 1.90 $\times$ 10-4 | TeV-1 cm-2 s-1 | 886
Table 2: The proportionality constants needed to calculate the expected gamma-
ray flux near a given pulsar. The dimensions for $A_{i}K$ are energy-1
distance-1 time-1. With the value of $K$ known and $A_{i}$ in the units given
by the last column, the expected gamma-ray flux can easily be calculated using
Equation 1. $K$ is reported at the pivot energy in each bin. Proportionality
constants are desrived from the full sample of 24 pulsars.

Satisfied that the chosen models are adequate, we then run the joint-
likelihood analysis using the downselected list of ten pulsars (as described
in Section 2.2) to search for sub-threshold leptonic emission from pulsars.
Results are given in Table 3. Each of the models give a total TS value above
9, corresponding to a detection of more than 3$\sigma$. For some models, the
TS is much higher. The 1/$d^{2}$ model gives the highest TS: 41.3
(6.4$\sigma$). The same two models as before, 1/$d^{2}$ and the gamma-ray flux
at 7 TeV, perform better than the “no model” case, but in this case 1/$d^{2}$
performs the best, with a TS slightly higher than the “gamma-ray flux at 7
TeV” scenario.

Model | TS (56 $<$ E $<$ 100 TeV) | TS (100 $<$ E $<$ 177 TeV) | TS (177 $<$ E $<$ 316 TeV) | Total TS
---|---|---|---|---
No model | 27.9 | 8.33 | 1.59 | 37.9
1/$d^{2}$ | 31.9 | 9.08 | 1.29 | 41.3
$\dot{E}$/$d^{2}$ | 9.58 | 5.24 | 0.00 | 14.8
Inverse age | 9.19 | 3.79 | 0.03 | 13.0
Flux at 7 TeV | 26.3 | 9.61 | 3.62 | 39.6
Table 3: The test statistic for the joint-likelihood analysis for each model,
using the ten sub-threshold sources. Note that $\kappa$ is fit individually in
each energy bin so the last column is not a joint TS.  Figure 2: The total
combined flux from the ten sub-threshold candidates for the scenario where the
pulsars have a relative contribution defined by 1/$d^{2}$, where $d$ is the
distance between the pulsar and the Earth. For the first two energy bins, the
TS $>$ 4 so Bayesian credible intervals (68$\%$ containment) are plotted. In
the last energy bin, there is no significant detection so a 95$\%$ upper limit
is plotted. The dark and light grey bands are 68$\%$ and 90$\%$ containment
for expected upper limits from Poisson fluctuated background. The dotted lines
are systematic uncertainties on the central value (i.e., the solid red line
may be as low or as high as the dotted line once systematic uncertainties are
included, see section 3.1)

Figure 2 shows the total flux from all ten sub-threshold candidates for this
best-case scenario. The figures for the other models can be seen in the
Appendix. Table 4, also located in the Appendix, gives the total flux
normalizations for each model.

The Appendix also includes a discussion of how often TS this high would be
expected from stacking random points on the sky.

###  3.1 Systematic Uncertainties

Systematic uncertainties are broken down into two categories: detector
systematics and modeling systematics. Detector systematics, described in
Abeysekara et al. (2019), stem from mis-modeling of detector quantities such
as the photomultiplier tube threshold and charge in simulated Monte Carlo
events. Each systematic is treated independently; the results are added in
quadrature to get a total uncertainty. This process is repeated in each energy
bin. Depending on the energy bin and model assumed, the detector systematic
ranges from 10$\%$ to 25$\%$.

Modeling systematic uncertainties investigate how the analysis would change if
the emission near the sub-threshold pulsars is different from what is assumed
in the main analysis. Several modeling systematics are considered:

  * •

Spectral indices in the power law ($\alpha$ in Equation 1) of 2.0 and 3.0

  * •

Replacing the power-law spectral model (Equation 1) with a power-law with an
exponential cutoff:

$\frac{dN}{dE}=A_{i}K\left(\frac{E}{E_{0}}\right)^{-\alpha}e^{-E/E_{cut}}.$
(5)

Here, $\alpha$ is fixed at 2.5 and $E_{cut}$ is fixed at 60 TeV, which are
average values from Abeysekara et al. (2020)

  * •

Changing the Gaussian width to 0.23∘ and 0.45∘. These are $\pm$ 1 standard
deviation from the average extension of the sources from HAWC’s eHWC catalog
(Abeysekara et al., 2020).

The modeling systematics are larger than the detector systematics, driven
predominantly by the assumed source size. In the last energy bin, assuming a
power law with an exponential cutoff instead of a power law is also a dominant
effect.

Depending on energy and model, the modeling systematic ranges from 13$\%$ to
34$\%$. The sum of the detector and modeling systematics, added in quadrature,
are denoted with a dotted black line in Figure 2 and in Figures 3 through 5 in
the Appendix.

As described in Albert et al. (2020b), the absolute pointing uncertainty of
HAWC is declination-dependent and may be larger than previously thought;
perhaps as large as 0.3∘ at the edges of the field-of-view. This means that
the TS values presented in this paper may be underestimated; the size of this
effect is $\sim$10$\%$. The flux may also be underestimated. The flux
underestimate is a function of energy and ranges from 9$\%$ in the 56-100 TeV
bin, decreasing to a negligible amount (0.6$\%$) in the 177-316 TeV bin.

##  4 Discussion

Regardless of which model is used, it is clear that the areas around
high-$\dot{E}$ pulsars show hints of emitting at ultra-high energies ($>$ 56
TeV). Interestingly, some models based on the pulsar parameters, such as
$\dot{E}$ or inverse age, have much lower TS values than the “no model” case,
implying that they do not describe the emission well. Two models that perform
better than the “no model” case are the gamma-ray flux at 7 TeV and 1/$d^{2}$.

It is unclear at this time why some pulsar parameters do not seem to be good
predictors of emission. Some of these parameters have fairly large
uncertainties, which could be a contributing factor. For example, the
characteristic age of the pulsar can differ from the true age.

However, the uncertainties on the pulsar distances are relatively small. For
the ATNF pulsar database, 95$\%$ of pulsars will have a relative error of less
than a factor of 0.9 in their distance estimate (10$\%$ uncertainty) (Yao et
al., 2017).

Alternatively, the individual environment that each pulsar is in could play a
large factor in the amount of UHE emission. While the pulsar itself is the
particle accelerator, diffusion of electrons and positrons away from the
pulsar is strongly dependent on quantities such as the density and magnetic
field of the environment. Note that this may be only true for the
high-$\dot{E}$ pulsars studied here, which are relatively young; this means
the magnetic fields are likely to be affected by the pulsar age and SNR
interactions. For weaker, older pulsars, the magnetic field density and
environment are instead likely dominated by the ISM.

Also, the emission mechanisms are still uncertain. While it is commonly
assumed that the bulk of emission from the vicinity of a pulsar is from a PWN
and therefore predominantly leptonic, a hadronic contribution cannot be a
priori discarded. While emission from associated SNRs are unlikely, some have
raised the possibility of more exotic hadronic emission mechanisms in or near
PWN (Amato et al., 2003; Di Palma et al., 2017). Several of the HAWC sources
known to emit above 56 TeV, most notably eHWC J1908+063 and eHWC J1825-134,
have molecular clouds nearby (Duvidovich et al., 2020; Voisin et al., 2016).
These molecular clouds may be serving as a target for a portion of the gamma-
ray production.

Multi-wavelength and multi-messenger campaigns can help disentangle emission
mechanisms. A neutrino detection coincident with one of these pulsars would be
a smoking gun for hadronic emission mechanisms in or near PWN. However, a
recent stacked analysis looking for neutrino emission from PWN by IceCube did
not yield a detection (Aartsen et al., 2020).

Electromagnetic multi-wavelength studies could also be helpful. A leptonic
source emitting above 56 TeV will have a different signature at lower energies
than a hadronic one. For example, 100 TeV gamma rays imply a synchrotron peak
in the keV regime, assuming a 3 $\mu$G field (Hinton & Hofmann, 2009), with
the emission extending up to the MeV energy range. If the emission is instead
predominantly hadronic, there will be no such peak at these energies. Proposed
experiments such as AMEGO (McEnery et al., 2019) will be important in
distinguishing emission mechanisms.

##  5 Conclusions

In this study, we have searched for UHE gamma-ray emission in the vicinity of
pulsars with an $\dot{E}>10^{36}$ erg/s. We find, with high significance ($>$
3$\sigma$) , that UHE gamma-ray emission is a generic feature in the vicinity
of this class of pulsars. 1/$d^{2}$ is the model that gives the highest TS; a
source that is closer to us is more likely to have observed UHE emission.
Other pulsar parameters do not seem to be good predictors of emission. This
implies that the environments the pulsars are located in may play a role in
the amount of emission.

The TS values obtained are higher than would be expected from combining random
points in the Galactic plane. There is the possibility that all known gamma-
ray sources emit above this energy threshold, albeit at an extremely low flux.

Multimessenger and multiwavelength studies are needed to disentagle the origin
of the UHE emission from these pulsars.

We acknowledge the support from: the US National Science Foundation (NSF); the
US Department of Energy Office of High-Energy Physics; the Laboratory Directed
Research and Development (LDRD) program of Los Alamos National Laboratory;
Consejo Nacional de Ciencia y Tecnología (CONACyT), México, grants 271051,
232656, 260378, 179588, 254964, 258865, 243290, 132197, A1-S-46288,
A1-S-22784, cátedras 873, 1563, 341, 323, Red HAWC, México; DGAPA-UNAM grants
IG101320, IN111315, IN111716-3, IN111419, IA102019, IN110621; VIEP-BUAP; PIFI
2012, 2013, PROFOCIE 2014, 2015; the University of Wisconsin Alumni Research
Foundation; the Institute of Geophysics, Planetary Physics, and Signatures at
Los Alamos National Laboratory; Polish Science Centre grant,
DEC-2017/27/B/ST9/02272; Coordinación de la Investigación Científica de la
Universidad Michoacana; Royal Society - Newton Advanced Fellowship 180385;
Generalitat Valenciana, grant CIDEGENT/2018/034; Chulalongkorn University’s
CUniverse (CUAASC) grant. Thanks to Scott Delay, Luciano Díaz and Eduardo
Murrieta for technical support.

## References

  * Aartsen et al. (2020) Aartsen, M. G., Ackermann, M., Adams, J., et al. 2020, The Astrophysical Journal, 898, 117, doi: 10.3847/1538-4357/ab9fa0
  * Abdalla et al. (2018a) Abdalla, H., Abramowski, A., Aharonian, F., et al. 2018a, Astronomy and Astrophysics, 612, A8, doi: 10.1051/0004-6361/201730737
  * Abdalla et al. (2018b) —. 2018b, Astronomy and Astrophysics, 612, A1
  * Abdalla et al. (2019) Abdalla, H., Aharonian, F., Benkhali, F. A., Angüner, E. O., et al. 2019, Astrononomy and Astrophysics, 621
  * Abeysekara et al. (2017) Abeysekara, A., Albert, A., Alfaro, R., et al. 2017, The Astrophysical Journal, 843, 39
  * Abeysekara et al. (2019) —. 2019, The Astrophysical Journal, 881, 134
  * Abeysekara et al. (2020) —. 2020, Phys. Rev. Lett., 124, 021102, doi: 10.1103/PhysRevLett.124.021102
  * Albert et al. (2020a) Albert, A., Alfaro, R., Alvarez, C., et al. 2020a, ApJ, 896, L29, doi: 10.3847/2041-8213/ab96cc
  * Albert et al. (2020b) —. 2020b, The Astrophysical Journal, 905, 76, doi: 10.3847/1538-4357/abc2d8
  * Amato et al. (2003) Amato, E., Guetta, D., & Blasi, P. 2003, Astronomy and Astrophysics, 402, 827, doi: 10.1051/0004-6361:20030279
  * Amenomori et al. (2019) Amenomori, M., Bao, Y. W., Bi, X. J., et al. 2019, Physical Review Letters, 051101, doi: 10.1103/PhysRevLett.123.051101
  * Breuhaus et al. (2020) Breuhaus, M., Hahn, J., Romoli, C., et al. 2020, arXiv e-prints, arXiv:2010.13960.  https://arxiv.org/abs/2010.13960
  * Di Palma et al. (2017) Di Palma, I., Guetta, D., & Amato, E. 2017, The Astrophysical Journal, 836, 159, doi: 10.3847/1538-4357/836/2/159
  * Duvidovich et al. (2020) Duvidovich, L., Petriella, A., & Giacani, E. 2020, MNRAS, 491, 5732, doi: 10.1093/mnras/stz3414
  * Gabici (2017) Gabici, S. 2017, in American Institute of Physics Conference Series, Vol. 1792, 6th International Symposium on High Energy Gamma-Ray Astronomy, 020002
  * Gotthelf et al. (2011) Gotthelf, E. V., Halpern, J. P., Terrier, R., & Mattana, F. 2011, Astrophys. J. Lett., 729
  * Hinton & Hofmann (2009) Hinton, J., & Hofmann, W. 2009, The Annual Review of Astronomy and Astrophysics, 47, 523
  * Huentemeyer et al. (2019) Huentemeyer, P., Abreu, P., Albert, A., et al. 2019, in Bulletin of the American Astronomical Society, Vol. 51, 109
  * Linden et al. (2017) Linden, T., Auchettl, K., Bramante, J., et al. 2017, Phys. Rev. D, 96, 103016, doi: 10.1103/PhysRevD.96.103016
  * Manchester et al. (2005) Manchester, R. N., Hobbs, G. B., Teoh, A., & Hobbs, M. 2005, The Astronomical Journal, 129, 1993, doi: 10.1086/428488
  * McEnery et al. (2019) McEnery, J., van der Horst, A., Dominguez, A., et al. 2019, in Bulletin of the American Astronomical Society, Vol. 51, 245
  * Moderski et al. (2005) Moderski, R., Sikora, M., Coppi, P. S., & Aharonian, F. 2005, Monthly Notices of the Royal Astronomical Society, 363, 954
  * Pletsch et al. (2012) Pletsch, H. J., Guillemot, L., Allen, B., et al. 2012, The Astrophysical Journal Letters, 755
  * Smith (2015) Smith, A. J. 2015, Proceedings of Science (34th International Cosmic Ray Conference), 966
  * Sudoh et al. (2021) Sudoh, T., Linden, T., & Hooper, D. 2021, arXiv e-prints, arXiv:2101.11026.  https://arxiv.org/abs/2101.11026
  * Vianello et al. (2015) Vianello, G., Lauer, R. J., Younk, P., et al. 2015, Proceedings of Science (34th International Cosmic Ray Conference), 1042
  * Voisin et al. (2016) Voisin, F., Rowell, G., Burton, M. G., et al. 2016, MNRAS, 458, 2813, doi: 10.1093/mnras/stw473
  * Wilks (1938) Wilks, S. S. 1938, The Annals of Mathematical Statistics, 9, 60, doi: 10.1214/aoms/1177732360
  * Wu et al. (2018) Wu, J., Clark, C. J., Pletsch, H. J., et al. 2018, The Astrophysical Journal, 854, 1, doi: 10.3847/1538-4357/aaa411
  * Yao et al. (2017) Yao, J. M., Manchester, R. N., & Wang, N. 2017, The Astrophysical Journal, 835, 29, doi: 10.3847/1538-4357/835/1/29
  * Zeng et al. (2021) Zeng, H., Xin, Y., Zhang, S., & Liu, S. 2021, arXiv e-prints, arXiv:2102.03465.  https://arxiv.org/abs/2102.03465

###  .1 Additional joint-likelihood results

Figures 3 through 5 are analogous to Figure 2 in the main text, but for the
other four models that have been investigated (described in Section 2.3). All
figures contain the combined flux for the joint-likelihood analysis of the ten
sub-threshold pulsars. Table 4 contains the total combined flux normalizations
for each of the models that have been considered.

Figure 3: Identical to Figure 2 from the main text, but for the “no model” case.  Figure 4: Identical to Figure 2 from the main text, but for the $\dot{E}/d^{2}$ model.  Figure 5: Identical to Figure 2 from the main text, but for the inverse age model.  Figure 6: Identical to Figure 2 from the main text, but for the model defined by the gamma-ray flux at 7 TeV.  Model |  $\kappa$ (56 $<$ E $<$ 100 TeV) |  $\kappa$ (100 $<$ E $<$ 177 TeV) |  $\kappa$ (177 $<$ E $<$ 316 TeV)
---|---|---|---
No model | 8.47${{}_{-1.78}^{+1.88}}$ $\times$ 10-16 | 1.57${}_{-0.56}^{+0.63}$ $\times$ 10-16 | 7.56 $\times$ 10-17
1/$d^{2}$ | 8.38${}_{-1.72}^{+1.76}$ $\times$ 10-16 | 1.54${}_{-0.55}^{+0.60}$ $\times$ 10-16 | 7.07 $\times$ 10-17
$\dot{E}$/$d^{2}$ | 3.45${}_{-1.18}^{+1.27}$ $\times$ 10-16 | 7.41${}_{-3.24}^{+4.10}$ $\times$ 10-17 | 3.39 $\times$ 10-17
Inverse age | 3.99${}_{-1.39}^{+1.48}$ $\times$ 10-16 | 1.59 $\times$ 10-16 | 4.48 $\times$ 10-17
Flux at 7 TeV | 7.80${}_{-1.69}^{+1.78}$ $\times$ 10-16 | 1.56${}_{-0.55}^{+0.57}$ $\times$ 10-16 | 7.98 $\times$ 10-17
Table 4: The total combined flux normalization for the sub-threshold joint-
likelihood analysis for each of the models ($\kappa$ from Equation 3). The
units are TeV-1 cm-2 s-1. The flux normalization is reported at the pivot
energy, which is the center of each bin. Values without uncertainties are
upper limits, otherwise the values correspond to the 68$\%$ containment for
the Bayesian credible interval. Uncertainties are statistical only.

###  .2 Testing of Random Backgrounds

We investigate how often randomly chosen non-source locations give TS values
as high as in the study presented in the main text. Sets of ten randomly
chosen points from the analysis region ($4^{\circ}$ $<l<$ 90∘;
$|b|<1^{\circ}$) are run through the joint-likelihood analysis. Points that
are within a degree of either a known high-energy source or any of the
selected ATNF pulsars are excluded. Because the gamma-ray emission is assumed
to be spatially extended, this is necessary to avoid contributions from the
known sources and the pulsars of interest.

This analysis is performed for 2,877 sets of ten random source locations. Only
five of these randomly chosen sets have a TS greater than 37.9. This is the TS
from the sub-threshold joint-likelihood “no model” case and is used as the
comparison because the other models require known pulsar information. This
means that the results that we have obtained are significant at the 3$\sigma$
level.

We also investigate sets of ten randomly chosen sources from HAWC’s third
catalog (3HWC) (Albert et al., 2020b) to see if any of HAWC’s previously
detected TeV sources emit at UHE. Once again, sources that are known to emit
above 56 TeV and sources within a degree of the ATNF pulsars used in the
nominal analysis are removed. Note that the majority of HAWC sources are
leptonic in origin (Linden et al., 2017).

Due to the relatively small number of remaining 3HWC sources after this
downselection (48 of the original 65 3HWC sources), it is not possible to run
thousands of sets of ten sources without repeating a large subset of the
sources. We instead run 100 sets of ten randomly chosen 3HWC sources.

None of these 100 trials have a TS above the value from the sub-threshold
joint-likelihood unmodeled case. The highest TS from this study is 26.9 and
the mean is 12.2 (standard deviation 5.5). This implies that known gamma-ray
sources that are not located near high-$\dot{E}$ pulsars are unlikely to emit
at high energies, or if they do, their fluxes are very low and combining ten
sources is not enough for a detection.

We also explore combining 20 3HWC sources at a time, instead of the ten that
were combined in the preceding paragraph. The entire TS distribution is
shifted to higher values. The highest TS is 45.9 (higher than in the “no
model” case); the mean is 24.6 and the standard deviation is 7.2. The higher
TS values when more 3HWC sources are combined implies that UHE emission may be
a generic feature of known gamma-ray sources, but with a very low flux. This
should be investigated, but may be hard to do with current-generation
experiments. Proposed experiments such as SWGO (Huentemeyer et al., 2019),
will be important here.