:
126 2021 Machine Learning for Healthcare

# MedAug: Contrastive learning leveraging patient metadata improves
representations for chest X-ray interpretation

Yen Nhi Truong Vu1*<EMAIL_ADDRESS>Richard Wang1*<EMAIL_ADDRESS>Niranjan Balachandar2*<EMAIL_ADDRESS>Can Liu1<EMAIL_ADDRESS>Andrew Y. Ng1<EMAIL_ADDRESS>Pranav Rajpurkar1<EMAIL_ADDRESS>

*Equal Contribution
1Department of Computer Science, Stanford University
2School of Medicine, Stanford University

###### Abstract

Self-supervised contrastive learning between pairs of multiple views of the
same image has been shown to successfully leverage unlabeled data to produce
meaningful visual representations for both natural and medical images.
However, there has been limited work on determining how to select pairs for
medical images, where availability of patient metadata can be leveraged to
improve representations. In this work, we develop a method to select positive
pairs coming from views of possibly different images through the use of
patient metadata. We compare strategies for selecting positive pairs for chest
X-ray interpretation including requiring them to be from the same patient,
imaging study or laterality. We evaluate downstream task performance by fine-
tuning the linear layer on 1% of the labeled dataset for pleural effusion
classification. Our best performing positive pair selection strategy, which
involves using images from the same patient from the same study across all
lateralities, achieves a performance increase of 14.4% in mean AUC from the
ImageNet pretrained baseline. Our controlled experiments show that the keys to
improving downstream performance on disease classification are (1) using
patient metadata to appropriately create positive pairs from different images
with the same underlying pathologies, and (2) maximizing the number of
different images used in query pairing. In addition, we explore leveraging
patient metadata to select hard negative pairs for contrastive learning, but
do not find improvement over baselines that do not use metadata. Our method is
broadly applicable to medical image interpretation and allows flexibility for
incorporating medical insights in choosing pairs for contrastive learning.

##  1 Introduction

Self-supervised contrastive learning has recently made significant strides in
enabling the learning of meaningful visual representations through unlabeled
data (Wu et al., 2018b; Hjelm et al., 2018; Chen et al., 2020b, a). In medical
imaging, previous work has found performance improvement when applying
contrastive learning to chest X-ray interpretation (Sowrirajan et al., 2020;
Sriram et al., 2021; Azizi et al., 2021), dermatology classification (Azizi et
al. (2021)) and MRI segmentation (Chaitanya et al., 2020). Despite the early
success of these applications, there is limited work on determining how to
improve upon standard contrastive algorithms using medical information
(Sowrirajan et al., 2020; Chaitanya et al., 2020; Zhang et al., 2020; Kiyasseh
et al., 2020).

In contrastive learning, the selection of pairs controls the information
contained in learned representations, as the loss function dictates that
representations of positive pairs are pulled together while those of negative
pairs are pushed apart (Oord et al., 2018). For natural images where there are
no other types of annotations, positive pairs are created using different
augmented views of the same image while negative pairs are views of different
images (Chen et al., 2020a). Tian et al. (2020) argue that good positive pairs
are those that contain minimal mutual information apart from common downstream
task information. In the natural image setting, Tamkin et al. (2020) train a
generative model which learns to produce multiple positive views from a single
input. However, previous contrastive learning studies on medical imaging have
not systematically investigated how to leverage patient metadata available in
medical imaging datasets to select positive pairs that go beyond augmentations
of the same image.

In this work, we propose a method to treat different images that share common
properties found in patient metadata as positive pairs in the context of
contrastive learning. We demonstrate the application of this method to a chest
X-ray interpretation task. Similar to the concurrent work by Azizi et al.
(2021), we experiment with requiring positive pairs to come from the same
patient as these images likely share highly similar pathological features.
However, our method incorporates these positive pairs with possibly different
images directly as part of the view generation scheme in a single contrastive
pretraining stage, as opposed to Azizi et al. (2021), which adds a second
pretraining stage where a positive pair must be formed by two distinct images.
Further, we go beyond the simple strategy of forming a positive pair using any
two data points coming from the same patient as in Azizi et al. (2021) and
Kiyasseh et al. (2020) and experiment with other metadata such as study number
and laterality to identify a pair of images that are likely to have the same
pathologies. Although study number has also been leveraged successfully in
Sriram et al. (2021) to create a sequence of pretrained embeddings
representating patient disease progression, our work differs in that we use
this information specifically to choose positive pairs during the contrastive
pretraining stage.

We conduct MoCo-pretraining (Chen et al., 2020b) using these different
criteria and evaluate the quality of pretrained representations by freezing
the base model and fine-tuning a linear layer using 1% of the labeled dataset
for the task of pleural effusion classification. Our contributions are:

  1. 1.

We develop a method, MedAug, to use patient metadata to select positive pairs
in contrastive learning, and apply this method to chest X-rays for the
downstream task of pleural effusion classification.

  2. 2.

Our best pretrained representation achieves a performance increase of 14.4% in
mean AUC compared to the ImageNet pretrained baseline, showing that using
patient metadata to select positive pairs from multiple images can notably
improve representations.

  3. 3.

We perform comparative empirical analysis with label information and show that
(1) using positive pairs of different images from the same patients that share
underlying pathologies improves pretrained representations, and (2) increasing
the number of distinct images selected to form positive pairs per image query
improves the quality of pretrained representations.

  4. 4.

We perform an exploratory analysis of strategies to select negative pairs
using patient metadata, and do not find improvement over the default strategy
that does not use metadata.

fig:method

Figure 1: Selecting positive pairs for contrastive learning with patient
metadata

### Generalizable Insights about Machine Learning in the Context of Healthcare

Our work presents methods to improve contrastive pretraining for medical image
interpretation. We evaluate our methods on the task of chest X-ray
classification, but our approach is generalizable to other medical imaging
tasks. First, we demonstrate that leveraging patient metadata – specifically
information about the patient, study, and laterality of the images – to
identify pathological similarity and define positive pairs for contrastive
learning improves pretrained representations. This can be applicable to tasks
in medical image interpretation where such metadata is available. Second, we
show that using two distinct images as positive pairs for contrastive learning
provides better pretrained representations than only using two augmentations
of the same image. Therefore, other contrastive learning methods for medical
imaging may benefit from expanding the construction of positive pairs by using
patient metadata to inform this expansion. Third, our analysis indicates that
increasing the number of images available for selecting positive pairs
improves pretrained representations. In general, in contrastive learning for
medical imaging, using domain-specific information to increase the pool of
images to select positive pairs may improve pretrained representations and
downstream task performance.

##  2 Cohort

We use CheXpert, a large collection of de-identified chest X-ray images (Irvin
et al., 2019). The dataset consists of 224,316 images from 65,240 patients
labeled for the presence or absence of 14 radiological observations. We use
these images for pretraining and random samples of 1% of these images for
fine-tuning. The test set consists of 500 additional labeled images from 500
studies not included in the training set. We perform fine-tuning experiments
on the downstream task of pleural effusion classification, which was selected
based on quality of ground truth labels, clinical importance and prevalence.
We evaluate model performance using test set AUC after fine-tuning.

##  3 Methods

###  3.1 Selecting positive pairs for contrastive learning with patient
metadata

Given an input image $\bm{x}$, encoder $g$, and a set of augmentations
$\mathcal{T}$, most contrastive learning algorithms involve minimizing the
InfoNCE loss

$\displaystyle\mathcal{L}(\bm{x})=-\log\frac{\exp[g(\bm{\tilde{x}}_{1})\cdot
g(\bm{\tilde{x}}_{2})]}{\exp[g(\bm{\tilde{x}}_{1})\cdot
g(\bm{\tilde{x}}_{2})]+\sum_{i=1}^{K}\exp[g(\bm{\tilde{x}}_{1})\cdot
g(\bm{z}_{i})]}.$ (1)

Here, the positive pair
$(\bm{\tilde{x}}_{1}=t_{1}(\bm{x}),\bm{\tilde{x}}_{2}=t_{2}(\bm{x}))$ with
$t_{1},t_{2}\in\mathcal{T}$ are augmentations of the input image $\bm{x}$,
while the negative pairs $(\bm{\tilde{x}}_{1},\bm{z}_{i})$, $1\leq i\leq K$
are pairs of augmentations of different images, with $\bm{z}_{i}$ coming from
either a queue in the case of MoCo or the minibatch in the case of SimCLR.
Recognizing that many augmentation strategies available for natural images are
not applicable to medical images, Sowrirajan et al. (2020) restrict
$\mathcal{T}$ to be the set of simple augmentations such as horizontal
flipping and random rotation between -10 to 10 degrees. As a result, their
method can be thought of as instance discrimination, as $\bm{\tilde{x}}_{1}$
and $\bm{\tilde{x}}_{2}$ must come from the same image input.

In this work, we propose MedAug, a method to use multiple images as a way to
increase the number of positive pair choices. Beyond the disease labels, we
can use patient metadata such as patient number, study number, laterality,
patient historical record etc. to create appropriate positive pairs. Formally,
we can use patient metadata to obtain an enhanced augmentation set dependent
on $\bm{x}$ as follows

$\displaystyle\mathcal{T}_{\text{enhanced}}(\bm{x})$
$\displaystyle=\begin{cases}\\{t_{i}(\bm{x}^{\prime})|t_{i}\in\mathcal{T},\bm{x}^{\prime}\in\mathcal{S}_{c}(\bm{x})\\}&\text{if
}\mathcal{S}_{c}(\bm{x})\neq\emptyset\\\
\mathcal{T}(\bm{x})&\text{otherwise}\end{cases}$ (2)

where $\mathcal{S}_{c}(\bm{x})$ is the set of all images satisfying some
predefined criteria $c$ in relation to the properties of $\bm{x}$. The
criteria for using the metadata could be informed by clinical insights about
the downstream task of interest.

We apply this method on chest X-ray interpretation and pretrain ResNet-18
models using MoCo v2 (Chen et al., 2020b) with learning rate of $10^{-4}$,
batch size of 16, temperature of 0.2, MLP projection, and 20 epochs. Since the
downstream task is disease classification, we experiment with using
$\mathcal{S}_{\text{same patient}}(\bm{x})$ since images from the same patient
are likely to share high amount of pathological features. We also experiment
with further applying criteria on study numbers as well as laterality. An
example application of the method is illustrated in Figure LABEL:fig:method.
Given an image query, the method collects all images from the same patient as
the query, chooses a subset of the collection using the given criteria, then
sample an image from the subset. Finally, an augmentation of the image
together with the augmentation of an initial query image form a positive pair.

###  3.2 Fine-tuning and evaluation

We evaluate the pretrained representations by (1) training a linear classifier
on outputs of the frozen encoder using labeled data and (2) end-to-end fine-
tuning. Pretrained checkpoints are selected with k-nearest neighbors algorithm
(Wu et al., 2018a) based on Faiss similarity search and clustering library
(Johnson et al., 2017). To simulate label scarcity encountered in medical
contexts, we fine-tune using only 1% of the labeled dataset. The fine-tuning
experiments are repeated on 5 randomly drawn 1% splits from the labeled
dataset to provide an understanding of the model’s performance variance. We
report the mean AUC and standard deviation over these five 1% fine-tuning
splits. Following Irvin et al. (2019), we use a learning rate of $3\times
10^{-5}$, batch size of 16 and 95 epochs for training.

##  4 Experiments

###  4.1 Positive pair selection

Our formulation of using any set of images $\mathcal{S}_{c}(\bm{x})$ from the
same patient to enhance the set of augmentations for contrastive learning
provides the flexibility of experimenting with different criteria $c$ for
constraining $\mathcal{S}_{c}(\bm{x})$. We experiment with limiting
$\mathcal{S}_{c}(\bm{x})$ using properties found in the metadata of the query
$\bm{x}$. In particular, we focus on two properties:

##### Study number.

The study number of an image associated with a particular patient reflects the
session in which the image was taken. We experiment with three different
criteria on study number:

  1. 1.

All studies: no restriction on $\mathcal{S}_{\text{all studies}}(\bm{x})$ is
dependent on the study number of $\bm{x}$

  2. 2.

Same study: only images from the same study with $\bm{x}$ belong to
$\mathcal{S}_{\text{same study}}(\bm{x})$

  3. 3.

Distinct studies: only images with different study number from $\bm{x}$ belong
to $\mathcal{S}_{\text{distinct studies}}(\bm{x})$

##### Laterality.

Chest X-rays can be of either frontal (AP/PA) view or lateral view.

  1. 1.

All lateralities: no restriction on $\mathcal{S}_{\text{all
lateralities}}(\bm{x})$ is dependent on the laterality of $\bm{x}$

  2. 2.

Same laterality: only images from the same laterality with $\bm{x}$ belong to
$\mathcal{S}_{\text{same laterality}}(\bm{x})$

  3. 3.

Distinct lateralities: only images with a different laterality from that of
$\bm{x}$ belongs to $\mathcal{S}_{\text{distinct lateralities}}(\bm{x})$

tab:positive  Baseline models Linear End-to-end ImageNet baseline $0.766\pm
0.009$ $0.858\pm 0.011$ MoCo v2 baseline (Sowrirajan et al., 2020) $0.847\pm
0.007$ $0.881\pm 0.017$ MoCo v2 baseline with random crop scale $0.864\pm
0.005$ $0.890\pm 0.026$ Criteria $c$ for creating $\mathcal{S}_{c}(\bm{x})$
Linear End-to-end Same patient, same study, same laterality $0.862\pm 0.004$
$0.894\pm 0.013$ Same patient, same study, distinct lateralities $0.865\pm
0.008$ $0.897\pm 0.008$ Same patient, same study $0.876\pm 0.013$ $0.902\pm
0.007$ Same patient, all studies $0.859\pm 0.006$ $0.877\pm 0.012$ Same
patient, distinct studies $0.848\pm 0.007$ $0.874\pm 0.013$ Same patient, same
study with random crop scale $\mathbf{0.883\pm 0.005}$ $\mathbf{0.906\pm
0.015}$

Table 1: Except for criteria $c$ that involve images from different studies,
using images from the same patient to select positive pairs result in improved
AUC in downstream pleural effusion classification.

##### Results.

We report the results of experiments using these criteria in Table
LABEL:tab:positive. Except from when $\mathcal{S}_{c}(\bm{x})$ includes images
with different study numbers from $\bm{x}$, where there is a drop in
performance, we see consistent large improvement from the baseline in
Sowrirajan et al. (2020). The best result is obtained when using
$\mathcal{S}_{\text{same study, all lateralities}}(\bm{x})$, the set of images
from the same patient and same study as that of $\bm{x}$, regardless of
laterality. Incorporating this augmentation strategy while holding other
settings from Sowrirajan et al. (2020) constant results in respective gains of
0.029 (3.4%) and 0.021 (2.4%) in AUC for the linear model and end-to-end
model. We also experiment with including random crop augmentation from MoCo v2
(Chen et al., 2020b), where the scaling is modified to be $[0.95,1.0]$ in
order to avoid cropping out areas of interest in the lungs. Adding this
augmentation to the same patient, same study strategy, we obtain our best
pretrained model, which achieves a linear fine-tuning AUC of 0.883 and an end-
to-end fine-tuning AUC of 0.906 on the test set, significantly outperforming
previous baselines. Results for repeated experiments across other CheXpert
competition tasks are included in Appendix A.

###  4.2 Comparative Empirical Analysis

We perform comparative analysis with labels to understand how different
criteria on patient metadata affect downstream performance results seen in
Table LABEL:tab:positive.

####  4.2.1 All studies v.s. same study

We hypothesize that the drop in transfer performance when moving from using
images with the same study number to using images regardless of study number
is because $\mathcal{S}_{\text{all studies}}(\bm{x})$ may contain images of a
different disease pathology than that seen in $\bm{x}$. As a result, the model
is asked to push the representation of a diseased image close to the
representation of a non-diseased image, causing poor downstream performance.
To test this hypothesis, we carry out an oracle experiment with
$\mathcal{S}_{\text{all studies, same label}}(\bm{x})$, the set of images from
the same patient and with the same downstream label as that of $\bm{x}$,
regardless of study number.

##### Results.

Table LABEL:tab:cheat_studies shows that the model pretrained with
$\mathcal{S}_{\text{all studies, same label}}(\bm{x})$ achieves a respective
improvement of 0.034 and 0.022 in AUC over $\mathcal{S}_{\text{all
studies}}(\bm{x})$ strategy for the linear model and end-to-end model. This
experiment supports our hypothesis that positive pairs from images with
different downstream labels hurt performance.

tab:cheat_studies  Criteria $c$ for creating $\mathcal{S}_{c}(\bm{x})$ Linear
End-to-end Same patient, all studies $0.859\pm 0.006$ $0.877\pm 0.012$ Same
patient, all studies, same disease label as $\bm{x}$ $\mathbf{0.893\pm 0.009}$
$\mathbf{0.899\pm 0.010}$

Table 2: Experiment with and without using downstream labels shows that
positive pairs with different labels hurt downstream classification
performance.

####  4.2.2 All studies v.s. distinct studies

There is a further performance drop when moving from using images across all
studies of the same patient to images with a different study number from the
current query image (Table LABEL:tab:positive). This finding may also support
our hypothesis because there is a larger proportion of positive pairs of
different disease pathologies in pairs of images from strictly different
studies. To make sure this result holds independently of the different number
of available images to form pair per query, we repeated these experiments
while forcing $\left|\mathcal{S}_{\text{same study, all
lateralities}}(\bm{x})\right|=\left|\mathcal{S}_{\text{same study, same
laterality}}(\bm{x})\right|$ via random subset pre-selection. Further, we only
use distinct images as a pair, i.e. skipping any $\bm{x}$ with
$\mathcal{S}_{c}(\bm{x})=\emptyset$ in (2) in order to remove any possible
contribution from positive pairs formed from the same image.

##### Results.

Table LABEL:tab:control_distinct_studies shows the same patient, all studies
strategy (AUC = 0.848) outperforms the same patient, distinct studies strategy
(AUC = 0.792) even when the size of $\mathcal{S}_{c}(\bm{x})$ is controlled.
This supports the hypothesis that a higher proportion of positive pairs with
different disease pathologies hurts downstream task performance. Downstream
performance from the $\mathcal{S}_{\text{distinct studies}}(\bm{x})$ is likely
lower than that of $\mathcal{S}_{\text{all studies}}(\bm{x})$ because there is
a higher proportion of positive pairs with different disease labels in
$\mathcal{S}_{\text{distinct studies}}(\bm{x})$. Figure 2 shows that there is
almost 9% of $\bm{x}$ where $\mathcal{S}_{\text{distinct studies}}(\bm{x})$
contains only images with a different disease label from $\bm{x}$, whereas
this scenario does not appear for $\mathcal{S}_{\text{all studies}}(\bm{x})$.

tab:control_distinct_studies  Criteria $c$ for creating $\mathcal{S}(\bm{x})$
Linear End-to-end Same patient, distinct studies $0.792\pm 0.007$ $0.841\pm
0.013$ Same patient, all studies (size controlled) $\mathbf{0.848\pm 0.009}$
$\mathbf{0.863\pm 0.010}$

Table 3: Experiments where we force positive pairs to come from different
images and control the size of $\mathcal{S}_{c}(\bm{x})$ shows that higher
proportion of pairs with different downstream labels contribute to lower
downstream performance. Figure 2: Histogram showing the distribution of the
proportions of positive pairs with different disease labels in
$\mathcal{S}_{\text{distinct studies}}(\bm{x})$ versus $\mathcal{S}_{\text{all
studies}}(\bm{x})$.

####  4.2.3 All lateralities v.s. distinct lateralities v.s. same laterality

First, we hypothesize that the drop in performance from the all lateralities
to the same laterality strategy could be due to $\mathcal{S}_{\text{same
study, same laterality}}(\bm{x})$ having smaller size. To test this, we carry
out an experiment in which $\mathcal{S}_{\text{same study, all
lateralities}}(\bm{x})$ is constrained by $\left|\mathcal{S}_{\text{same
study, same laterality}}(\bm{x})\right|$, the number of images from the same
study and has the same laterality as $\bm{x}$.

tab:control_laterality_1  Criteria $c$ for creating $\mathcal{S}_{c}(\bm{x})$
Linear End-to-end Same patient, same study, same laterality $0.862\pm 0.004$
$0.894\pm 0.013$ Same patient, same study, all lateralities (size controlled)
$0.860\pm 0.004$ $0.899\pm 0.011$ Same patient, same study, all lateralities
(no control) $\mathbf{0.876\pm 0.013}$ $\mathbf{0.902\pm 0.007}$

Table 4: Experiments with all lateralities where we control the size of
$S_{\text{same study, all lateralities}}$ show that the size of
$\mathcal{S}_{c}(\bm{x})$ affects downstream performance.

tab:control_laterality_2  Criteria $c$ for creating $\mathcal{S}_{c}(\bm{x})$
Linear End-to-end Same patient, same study, same laterality $0.856\pm 0.016$
$0.878\pm 0.015$ Same patient, same study, distinct lateralities
$\mathbf{0.866\pm 0.015}$ $\mathbf{0.882\pm 0.017}$

Table 5: Experiments to compare same v.s. distinct lateralities with size
restriction on $\mathcal{S}_{c}(\bm{x})$ shows no significant difference.

Our second hypothesis is that mutual information in images with different
lateralities is lower, which benefits retaining only information important to
the downstream task, as shown in Tian et al. (2020). We test this by training
two models on images that include at least one counterpart from the other
laterality. We pretrain one model with $\mathcal{S}_{\text{same study, same
laterality}}(\bm{x})$ containing only images with the same laterality as
$\bm{x}$, and the other model with $\mathcal{S}_{\text{same study, distinct
lateralities}}(\bm{x})$ containing only images with different laterality from
$\bm{x}$. To prevent the effect of different sizes of
$\mathcal{S}_{c}(\bm{x})$, we force that $\left|\mathcal{S}_{\text{same study,
same laterality}}(\bm{x})\right|=\left|\mathcal{S}_{\text{same study, distinct
lateralities}}(\bm{x})\right|$ via random subset pre-selection.

##### Results.

Table LABEL:tab:control_laterality_1 shows that once we control for the size
of $\mathcal{S}_{c}(\bm{x})$, there is no significant difference between using
images from the same laterality (AUC = 0.862) or from all lateralities (AUC =
0.860). However, the model pretraining with all images from all lateralities
achieves much larger downstream AUC of 0.876. Thus, it supports our first
hypothesis that the size of $\mathcal{S}_{c}(\bm{x})$ influences pretrained
representation quality. Table LABEL:tab:control_laterality_2 shows that once
we control for the size of $\mathcal{S}_{c}(\bm{x})$, the model pretrained
with images from different lateralities only gain $0.010$ AUC in linear fine-
tuning performance and a non-significant $0.004$ in end-to-end performance.
This experiment shows that the effect of mutual information from different
lateralities on pretrained representation quality is less pronounced.

###  4.3 Negative pair selection

We explore strategies using metadata in the CheXpert dataset to define
negative pairs. Similar to our method of defining positive pairs, we take
advantage of metadata available in the dataset to select the negative pairs.
However, unlike positive pair selection, where only a single pair is required
for each image, an image has to pair with the entire queue to select negative
pairs. This property makes selecting negative pairs from the same patient as
done in selecting positive pairs not suitable because only a small number of
images are available for a patient. We instead use a more general property –
laterality – across the patients to define negative pairs to retain sufficient
negative pairs in the loss function (1). Similarly, other metadata such as age
and sex may be exploited for the same purpose.

The default negative pair selection strategy is to select all keys from the
queue that are not views of the query image. However, we hypothesize that
negative pairs with the same laterality are “hard” negative pairs that are
more difficult to distinguish and provide more accurate pretrained
representations for the downstream task. We describe our four strategies
briefly as follows and in more detail in Appendix B. Our first strategy is to
only select images from the queue with the same laterality as the query to
create negative pairs. Our second strategy is to reweight the negative logits
based on laterality so in effect queries with each laterality (frontal and
lateral) equally contribute to the loss and the queue size remains fixed as in
the original MoCo approach. Following a similar idea in Kalantidis et al.
(2020), our third strategy is to sample a portion of negative pairs with the
same laterality for each query and append them to the queue for loss
computation. Our fourth strategy is to create synthetic negatives for
additional hard negative pairs. Unlike Kalantidis et al. (2020), we do not
determine hardness of negative pairs based on similarities of representations.
Instead, we use existing metadata (image laterality) to approximate hardness
of an negative pair. We evaluate the performance of each of these negative
pair strategies combined with the positive pair strategy of “same patient,
same study, all lateralities”.

##### Results.

Results are given in Table LABEL:tab:negative_results. The default negative
pair selection strategy (AUC = 0.876) is not outperformed by any of the
metadata-exploiting negative pair selection strategies including same
laterality only (AUC = 0.872), same laterality reweighted (AUC = 0.864 ), same
laterality appended (AUC = 0.875) and same laterality synthetic (AUC = 0.870).
Thus, our exploratory analysis does not indicate sufficient evidence for
performance improvement using strategies that incorporate metadata, but
further experiments with other metadata sources may be required to further
understand this relationship.

tab:negative_results  Negative Pairs Strategy Linear Default $\mathbf{0.876\pm
0.013}$ Same Laterality only $0.872\pm 0.011$ Same Laterality (reweighted)
$0.864\pm 0.006$ Same Laterality (appended) $0.875\pm 0.006$ Same Laterality
(synthetic) $0.870\pm 0.004$

Table 6: Experiments with the default negative pair definition (different
images) and various negative pair selection strategies.

##  5 Discussion

We introduce MedAug, a method to use patient metadata to select positive pairs
for contrastive learning, and demonstrate the utility of this method on a
chest X-ray interpretation task.

Can we improve performance by leveraging metadata to choose positive pairs?
Yes. Our best pretrained strategy with multiple images from the same patient
and same study obtains an increase of 3.4% in linear fine-tuning AUC in
comparison to the instance discrimination approach implemented in Sowrirajan
et al. (2020). A similar result has been shown by Kiyasseh et al. (2020) for
ECG signal interpretation. Azizi et al. (2021) also found improvement in
dermatology classification when applying a second contrastive pretraining
stage where strictly distinct images from the same patient are selected as
positive pairs.

Unlike previous work, our empirical analysis on using images from all studies
and distinct studies shows that simply choosing images from the same patient
may hurt downstream performance. We show that using appropriate metadata such
as study number to select positive pairs that share underlying disease
information is needed to obtain the best representation for the downstream
task of disease classification. For future studies, it is of interest to
experiment with other metadata such as age group, medical history, etc. and
how they can inform on tasks other than disease classification.

Our analysis using different criteria on laterality shows that the number of
images selected to form positive pairs plays an important role, while the
effect of mutual information is less clear. Given time and resources, it would
be informative to experiment with how the maximum number of distinct images
chosen per query affect downstream performance.

Can we improve performance by leveraging metadata to choose hard negative
pairs? Not necessarily. We perform an exploratory analysis of strategies to
leverage patient metadata to select negative pairs, and do not find them to
outperform the baseline.

In closing, our work demonstrates the potential benefits of incorporating
patient metadata into self-supervised contrastive learning for medical images,
and can be extended to a broader set of tasks (Rajpurkar et al., 2020;
Uyumazturk et al., 2019).

##### Limitations

In the CheXpert dataset, any two images from the same patient and study will
always have the same set of ground truth pathology labels. For medical image
datasets with different constraints regarding patient metadata, it remains
future work to determine the positive pair selection strategies that are
clinically relevant and produce good pretrained representations. Our approach
is not applicable to datasets lacking patient metadata altogether. For
datasets with limited data per patient, future work could be to cluster data
using the images and available metadata into larger groups, and define
positive pairs based on cluster assignments. While the results show that our
metadata-based contrastive learning methods are generalizable across all
CheXpert competition tasks, it remains future work to assess performance in
other datasets.

Furthermore, our strategies for negative pair selection did not improve
pretrained representations. Our strategies leveraged information regarding
image laterality. However, future work is required to whether negative pair
selection strategies using other metadata such as image view (anteroposterior
or posteroanterior), patient age or patient sex, or strategies using
similarity metrics can improve negative pair selection.

## References

  * Azizi et al. (2021) Shekoofeh Azizi, Basil Mustafa, Fiona Ryan, Zachary Beaver, Jan Freyberg, Jonathan Deaton, Aaron Loh, Alan Karthikesalingam, Simon Kornblith, Ting Chen, et al.  Big self-supervised models advance medical image classification.  _arXiv preprint arXiv:2101.05224_ , 2021.
  * Chaitanya et al. (2020) Krishna Chaitanya, Ertunc Erdil, Neerav Karani, and Ender Konukoglu.  Contrastive learning of global and local features for medical image segmentation with limited annotations.  _arXiv preprint arXiv:2006.10511_ , 2020.
  * Chen et al. (2020a) Ting Chen, Simon Kornblith, Kevin Swersky, Mohammad Norouzi, and Geoffrey Hinton.  Big self-supervised models are strong semi-supervised learners.  _arXiv preprint arXiv:2006.10029_ , 2020a.
  * Chen et al. (2020b) Xinlei Chen, Haoqi Fan, Ross Girshick, and Kaiming He.  Improved baselines with momentum contrastive learning.  _arXiv preprint arXiv:2003.04297_ , 2020b.
  * Hjelm et al. (2018) R Devon Hjelm, Alex Fedorov, Samuel Lavoie-Marchildon, Karan Grewal, Phil Bachman, Adam Trischler, and Yoshua Bengio.  Learning deep representations by mutual information estimation and maximization.  _arXiv preprint arXiv:1808.06670_ , 2018.
  * Irvin et al. (2019) Jeremy Irvin, Pranav Rajpurkar, Michael Ko, Yifan Yu, Silviana Ciurea-Ilcus, Chris Chute, Henrik Marklund, Behzad Haghgoo, Robyn Ball, Katie Shpanskaya, et al.  Chexpert: A large chest radiograph dataset with uncertainty labels and expert comparison.  In _Proceedings of the AAAI Conference on Artificial Intelligence_ , volume 33, pages 590–597, 2019.
  * Johnson et al. (2017) Jeff Johnson, Matthijs Douze, and Hervé Jégou.  Billion-scale similarity search with gpus.  _arXiv preprint arXiv:1702.08734_ , 2017.
  * Kalantidis et al. (2020) Yannis Kalantidis, Mert Bulent Sariyildiz, Noe Pion, Philippe Weinzaepfel, and Diane Larlus.  Hard negative mixing for contrastive learning, 2020.
  * Kiyasseh et al. (2020) Dani Kiyasseh, Tingting Zhu, and David A Clifton.  Clocs: Contrastive learning of cardiac signals.  _arXiv preprint arXiv:2005.13249_ , 2020.
  * Oord et al. (2018) Aaron van den Oord, Yazhe Li, and Oriol Vinyals.  Representation learning with contrastive predictive coding.  _arXiv preprint arXiv:1807.03748_ , 2018.
  * Rajpurkar et al. (2020) Pranav Rajpurkar, Allison Park, Jeremy Irvin, Chris Chute, Michael Bereket, Domenico Mastrodicasa, Curtis P Langlotz, Matthew P Lungren, Andrew Y Ng, and Bhavik N Patel.  Appendixnet: Deep learning for diagnosis of appendicitis from a small dataset of ct exams using video pretraining.  _Scientific reports_ , 10(1):1–7, 2020.
  * Sowrirajan et al. (2020) Hari Sowrirajan, Jingbo Yang, Andrew Y Ng, and Pranav Rajpurkar.  Moco pretraining improves representation and transferability of chest x-ray models.  _arXiv preprint arXiv:2010.05352_ , 2020.
  * Sriram et al. (2021) A Sriram, M Muckley, K Sinha, F Shamout, J Pineau, KJ Geras, L Azour, Y Aphinyanaphongs, N Yakubova, and W Moore.  Covid-19 prognosis via self-supervised representation learning and multi-image prediction.  2021\.
  * Tamkin et al. (2020) Alex Tamkin, Mike Wu, and Noah Goodman.  Viewmaker networks: Learning views for unsupervised representation learning, 2020.
  * Tian et al. (2020) Yonglong Tian, Chen Sun, Ben Poole, Dilip Krishnan, Cordelia Schmid, and Phillip Isola.  What makes for good views for contrastive learning.  _arXiv preprint arXiv:2005.10243_ , 2020.
  * Uyumazturk et al. (2019) Bora Uyumazturk, Amirhossein Kiani, Pranav Rajpurkar, Alex Wang, Robyn L. Ball, Rebecca Gao, Yifan Yu, Erik Jones, Curtis P. Langlotz, Brock Martin, Gerald J. Berry, Michael G. Ozawa, Florette K. Hazard, Ryanne A. Brown, Simon B. Chen, Mona Wood, Libby S. Allard, Lourdes Ylagan, Andrew Y. Ng, and Jeanne Shen.  Deep learning for the digital pathologic diagnosis of cholangiocarcinoma and hepatocellular carcinoma: Evaluating the impact of a web-based diagnostic assistant, 2019.
  * Wu et al. (2018a) Zhirong Wu, Yuanjun Xiong, Stella Yu, and Dahua Lin.  Unsupervised feature learning via non-parametric instance-level discrimination, 2018a.
  * Wu et al. (2018b) Zhirong Wu, Yuanjun Xiong, Stella X Yu, and Dahua Lin.  Unsupervised feature learning via non-parametric instance discrimination.  In _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition_ , pages 3733–3742, 2018b.
  * Zhang et al. (2020) Yuhao Zhang, Hang Jiang, Yasuhide Miura, Christopher D Manning, and Curtis P Langlotz.  Contrastive learning of medical visual representations from paired images and text.  _arXiv preprint arXiv:2010.00747_ , 2020.

##  Appendix A Repetition Across Tasks

We repeated each linear fine-tuning experiment for the other 4 CheXpert
competition classification tasks: Atelectasis, Cardiomegaly, Consolidation,
and Edema. The results are given in Table LABEL:tab:new. For all 5 CheXpert
competition classification tasks (see Table LABEL:tab:positive for Pleural
Effusion), we found that our best metadata-based approach to contrastive
learning, which was the “same patient, same study” positive pair selection
criterion, resulted in higher downstream classification performance than did
the ImageNet baseline and MoCo v2 baselines.

tab:new

Table 7: Linear fine-tuning AUCs for baseline methods and positive pair
selection criteria repeated across the tasks of Atelectasis, Cardiomegaly,
Consolidation, and Edema classification.

Baseline models | Atelectasis | Cardiomegaly | Consolidation | Edema
---|---|---|---|---
ImageNet baseline | 0.612 | 0.634 | 0.615 | 0.810
MoCo v2 baseline (Sowrirajan et al., 2020) | 0.671 | 0.735 | 0.699 | 0.847
MoCo v2 baseline with random crop scale | 0.708 | 0.730 | 0.693 | 0.858
Criteria $c$ for creating $\mathcal{S}_{c}(\bm{x})$ | Atelectasis | Cardiomegaly | Consolidation | Edema
| Same patient, same study, same laterality
---
0.738 | 0.718 | 0.672 | 0.866
Same patient, same study, distinct lateralities | 0.712 | 0.800 | 0.703 | 0.883
| Same patient, same study
---
0.762 | 0.785 | 0.721 | 0.889
| Same patient, all studies
---
0.727 | 0.786 | 0.805 | 0.879
| Same patient, distinct studies
---
0.676 | 0.757 | 0.760 | 0.846
| Same patient, same study with random crop scale
---
0.721 | 0.779 | 0.801 | 0.866

##  Appendix B Negative Pairs

Following the loss function in equation (1), we denote the exponential sum of
the negative pairs by $G$

$\displaystyle\mathcal{L}(\bm{x})=-\log\frac{\exp[g(\bm{\tilde{x}}_{1})\cdot
g(\bm{\tilde{x}}_{2})]}{\exp[g(\bm{\tilde{x}}_{1})\cdot
g(\bm{\tilde{x}}_{2})]+G\left(\bm{\tilde{x}}_{1},\bm{z}_{i}\right)}.$ (3)

where

$\displaystyle G\left(\bm{\tilde{x}}_{1},\bm{z}_{i}\right)=\sum_{\bm{z}_{i}\in
Q}\exp[g(\bm{\tilde{x}}_{1}))\cdot g(\bm{z}_{i}))]$ (4)

We follow the MoCo setup and denote $Q$ as the queue. Let
$\mathcal{S}(\bm{x})$ be the set of image representations in $Q$ that have the
same laterality as $\bm{x}$. We use the symbol $\mathbin{\|}$ to denote list
concatenation. We describe each of our negative pair selection strategies as
follows:

  1. 1.

(Same laterality only) For each query, we select keys in the queue that have
the same laterality as the query. Specifically, we replace $G$ in equation (4)
by $G^{l}$

$\displaystyle
G^{l}\left(\bm{\tilde{x}}_{1},\bm{z}_{i}\right)=\sum_{\bm{z}_{i}\in\mathcal{S}(\bm{x})}\exp[g(\bm{\tilde{x}}_{1}))\cdot
g(\bm{z}_{i}))]$ (5)

  2. 2.

(Same laterality reweighted) The first strategy excluded the keys in the queue
that have different laterality from the query. Here we set a target hard
negative weight and reweight each $exp$ term to achieve the target weight. Let

$\displaystyle
G^{w}\left(\bm{\tilde{x}}_{1},\bm{z}_{i}\right)=\sum_{\bm{z}_{i}\in\mathcal{S}(\bm{x})}w_{i}^{s}\exp[g(\bm{\tilde{x}}_{1}))\cdot
g(\bm{z}_{i}))]+\sum_{\bm{z}_{i}\in\mathcal{S}\left(x\right)^{c}}w_{i}^{d}\exp[g(\bm{\tilde{x}}_{1}))\cdot
g(\bm{z}_{i}))]$ (6)

where $t$ is the target hard negative weight and
$r=\frac{|\mathcal{S}(\bm{x})|}{|Q|}$ is the proportion of the negative keys
in the queue that have the same laterality as $\bm{x}$. Then
$w_{i}^{d}=\frac{1-t}{1-r}$ and $w_{i}^{s}=\frac{t}{r}$ for all $i$. In our
experiments, we set $t=0.1$ to allocate 90% of the weight to hard negatives.
This allows us to include all negative pairs in the contrastive loss but place
emphasis on hard negative pairs with the same laterality.

  3. 3.

(Same laterality appended) For each query, we select a random sample of the
keys that have the same laterality and append them to the existing queue

$Q=[z_{1},z_{2},...,z_{K}]$

where $K$ is the queue size. Let

$\displaystyle A=\left\\{z_{i_{1}},z_{i_{2}},...,z_{i_{m}}\right\\}\subset
S\left(x\right)$ (7)

be the random sample of keys with the same laterality as the query. The new
queue is

$Q^{a}=Q\mathbin{\|}A$

and

$G^{a}\left(\bm{\tilde{x}}_{1},\bm{z}_{i}\right)=\sum_{\bm{z}_{i}\in
Q^{a}}\exp[g(\bm{\tilde{x}}_{1}))\cdot g(\bm{z}_{i}))]$

replaces $G$ in equation (4).

  4. 4.

(Same laterality synthetic) For each query, in addition to appending samples
of the keys from $S(\bm{x})$, we use the samples to generate synthetic keys
and append them to the queue. We randomly sample $m$ pairs
$\left(\bm{s}_{i},\bm{s}_{j}\right)\in
A=\left\\{\bm{z}_{i_{1}},\bm{z}_{i_{2}},...,\bm{z}_{i_{m}}\right\\}$ and call
this set of pairs $B$.

For each pair $\left(\bm{s}_{i},\bm{s}_{j}\right)\in B$, we uniformly sample a
number $u$ between 0 and 1 and let

$h=u\cdot\bm{s}_{i}+(1-u)\cdot\bm{s}_{j}$

A synthetic image representation is defined as the normalized vector
$\frac{h}{\|h\|}$. Let $H$ be the set of these $m$ synthetic image
representations and

$Q^{h}=Q\mathbin{\|}B\mathbin{\|}H$

is the new queue. $G$ in equation (1) is replaced by

$G^{h}\left(\bm{\tilde{x}}_{1},\bm{z}_{i}\right)=\sum_{z_{i}\in
Q^{h}}\exp[g(\bm{\tilde{x}}_{1}))\cdot g(\bm{z}_{i}))]$

Note that unlike Kalantidis et al. (2020), we only construct synthetic images
once.


# Path integral molecular dynamics for anyons, bosons and fermions

Yunuo Xiong  College of Science, Zhejiang University of Technology, Hangzhou
31023, China  Hongwei Xiong<EMAIL_ADDRESS>College of Science, Zhejiang
University of Technology, Hangzhou 31023, China

###### Abstract

In this article we develop a general method to numerically calculate physical
properties for a system of anyons with path integral molecular dynamics. We
provide a unified method to calculate the thermodynamics of identical bosons,
fermions and anyons. Our method is tested and applied to systems of anyons,
bosons and fermions in a two-dimensional harmonic trap. We also consider a
method to calculate the energy for fermions as an application of the path
integral molecular dynamics to simulate the anyon model.

Keywords: anyons; bosons; fermions; path integral; molecular dynamics

##  I introduction

Quantum many-body systems are fascinating and reveal much of the interesting
phenomenon of nature. In particular, systems consisting of identical particles
are of great interests since quantum mechanics dictates that fundamental
exchange symmetries be present in such systems. Insightful effects such as
Berezinskii-Kosterlitz-Thouless transition Kosterlitz ; Hadzibabic  are found
resulting from the exchange symmetry between identical particles.
Traditionally, people have considered only two kinds of fundamental exchange
symmetry in a system of identical particles, with the wave function either be
completely symmetric or anti-symmetric, corresponding to bosons and fermions
respectively. However, when studying quantum many-body systems in two
dimensions, researchers have found that other kinds of exchange symmetry are
possible Leinaas . This leads to the discovery of anyons Leinaas ; FluxW ;
a-wilczek ; Wilczek ; khare . Anyons manifest novel phenomenon such as
fractional quantum Hall effects Tsui  and their existences in two dimensional
systems are shown by recent experiments Bar ; Nak ; Nak2 .

Aside from experimental study of many-body quantum systems, numerical
simulations of such systems are also desirable. There are two main ways to
perform the simulations based on path integral formalism feynman ; kleinert :
path integral Monte Carlo CeperRMP ; boninsegni1 ; boninsegni2 ; Dornheim  and
path integral molecular dynamics (PIMD) chandler ; Parrinello ; Tuckerman .
Path integral Monte Carlo has been applied successfully to large bosonic
systems boninsegni1 ; boninsegni2 , by sampling directly the partition
function corresponding to the system under study. However, PIMD has found few
applications in the field of identical particles until recently, when an
efficient polynomial algorithm Hirshberg  to perform PIMD has been developed,
enabling applications of PIMD to large bosonic system Hirshberg ; Xiong ;
XiongSpinor  and fermionic system HirshbergFermi ; Xiong1 . Both of these
approaches suffer from severe efficiency issues when applied to identical
fermions, due to the so-called fermion sign problem ceperley ; loh ; troyer ;
lyubartsev ; vozn ; Yao1 ; Yao2 ; Yao3 .

In this paper we extend PIMD to perform numercial simulations for a system of
identical anyons. In particular, we calculate the energies and densities for
anyons and study their dependences on the statistical properties for the
anyons. Our work provided a unified method to calculate various properties of
anyons, bosons and fermions. The paper is organized as follows: In Sec. II, we
start by reviewing path integral formulations for anyons. We then present PIMD
and show how to extend the technique to perform simulations for a system of
identical anyons. In Sec. III, we present our simulation results for the
energies and densities for anyons and study their dependences on the
statistical properties of the anyons. In Sec. IV, by taking advantage of the
extra freedom in choosing the statistical properties for anyons, we suggest a
way to calculate the energy of a few fermions by an extrapolation with a
simple function based on reasonable physical considerations. In Sec. V, we
give a concise summary and discussion of our work.

##  II theory

###  II.1 Path integral for anyons, bosons and fermions

In this paper, we will use "boson" gauge to consider the anyons in a general
way Myrheim . The derivation from Bose-Einstein statistics in this bosonic
gauge is due to a "statistics" field, which is a vector potential based on the
flux-tube model for anyons by Wilczek FluxW . The Bose statistics may change
continuously into Fermi statistics when the statistics angle $\theta$ changes
from boson value $0$ to fermion value $\pi$. The fractional statistics is the
situation of $0<\theta<\pi$. In this case, it is anticipated that there would
be a general numerical method which can be applied to both bosons and fermions
in a unified way.

Figure 1: The flux-tube model for anyons.

In the flux-tuble model for anyons FluxW , the anyons turn out to be point
particles having a fictitious magnetic flux. In this case, the flux is located
exactly at the position of the $j$th particle, and the fictitious magnetic
field for this particle is

$B(\mathbf{r})=-2\pi\nu\delta(\mathbf{r}-\mathbf{r}_{j}).$ (1)

Here $\nu$ is a parameter for the statistical properties of the anyon. $\nu=0$
for bosons and $\nu=1$ for fermions. The statistics angle $\theta=\nu\pi$. For
$N$ anyons, the total magnetic field is

$B(\mathbf{r})=-2\pi\nu\sum_{j=1}^{N}\delta(\mathbf{r}-\mathbf{r}_{j}).$ (2)

From the relation between $B$ and the vector potential $\mathbf{A}$

$B=\partial_{x}A_{y}-\partial_{y}A_{x},$ (3)

we have the following vector potential at $\mathbf{r}$ generated by the $j$th
particle,

$A_{x}^{j}(\mathbf{r})=\nu\frac{y-y_{j}}{(x-x_{j})^{2}+(y-y_{j})^{2}},$
$A_{y}^{j}(\mathbf{r})=-\nu\frac{x-x_{j}}{(x-x_{j})^{2}+(y-y_{j})^{2}}.$ (4)

The N-particle Hamiltonian is

$\hat{H}_{N}=\frac{1}{2m}\sum_{j=1}^{N}\left(\hat{\mathbf{p}}_{j}-\hbar{\mathbf{a}}_{j}\right)^{2}+V.$
(5)

Here $\hat{\mathbf{p}}_{j}$ is the canonical momentum of the $j$th particle.
$V$ is the interaction potential including both the external potential and
interaction potential between anyons. The statistical interaction is due to
the vector potential ${\mathbf{a}}_{j}$ for the $j$th particle. Because this
vector potential is due to all other particles, we have

$a_{jx}(\mathbf{r}_{j})=\nu\sum_{k\neq
j}\frac{y_{j}-y_{k}}{(x_{j}-x_{k})^{2}+(y_{j}-y_{k})^{2}},$
$a_{jy}(\mathbf{r}_{j})=-\nu\sum_{k\neq
j}\frac{x_{j}-x_{k}}{(x_{j}-x_{k})^{2}+(y_{j}-y_{k})^{2}}.$ (6)

###  II.2 The path integral for anyons

The partition function for $N$ anyons is

$Z_{A}(\beta)=Tre^{-\beta\hat{H}_{N}}.$ (7)

Here $\beta=1/k_{B}T$ with $k_{B}$ the Boltzmann constant and $T$ the
temperature of the system. Because we use the bosonic gauge to consider the
anyons, we have

$Z_{A}(\beta)=\frac{1}{N!}\sum_{p\in S_{N}}\int d\mathbf{r}_{1}\cdots
d\mathbf{r}_{N}\left<p\\{\mathbf{r}_{1},\cdots,\mathbf{r}_{N}\\}\right|e^{-\beta\hat{H}_{N}}\left|\mathbf{r}_{1},\cdots,\mathbf{r}_{N}\right>.$
(8)

Here $p$ denotes $N!$ permutation of $S_{N}$ permutation group. The integral
is about the usual $\mathbf{R}^{dN}$ space.

We divide $\beta$ into $P$ intervals (let us denote the interval length by
$\Delta\beta=\beta/P$) and expand each of the evolution operators into
products of the form

$e^{-\beta\hat{H}}=e^{-\Delta\beta\hat{H}}\cdots e^{-\Delta\beta\hat{H}}.$ (9)

In the limit of $P\rightarrow\infty$, we have

$Z_{A}(\beta)=\frac{1}{N!}\sum_{p\in
S_{N}}\int_{\mathcal{C}(p)}\mathcal{D}[\mathbf{r}_{1}(\tau),\cdots,\mathbf{r}_{N}(\tau)]e^{-\frac{S}{\hbar}-i\theta
Q}.$ (10)

Here

$S=\int_{0}^{\hbar\beta}d\tau\left(\frac{m}{2}\sum_{j=1}^{N}\left|\frac{d\mathbf{r}_{j}}{d\tau}\right|^{2}+V(\mathbf{r}_{1}\cdots\mathbf{r}_{N})\right),$
(11)

and

$-\theta Q=-\pi\nu Q=\sum_{j=1}^{N}\int
d\mathbf{r}_{j}\cdot\mathbf{a}_{j}(\mathbf{r}_{j}).$ (12)

This equation stands for the sum of a series of line integrals. To evaluate
such line integrals some parametrization of the curve must be chosen. Of
course the variable $\tau$ provides a natural parametrization and this was
used in our numerical calculation, but generally any parametrization works and
leads to the same result for the integral. We chose to write the line integral
in the general form in Eq. (12), whereas the discrete form based on $\tau$
parametrization will be presented in Eq. (15). In the above equation for
$Z_{A}(\beta)$, the domain of integration $\mathcal{C}(p)$ represents
$\\{\mathbf{r}_{1}(\hbar\beta),\cdots\mathbf{r}_{N}(\hbar\beta)\\}=p\\{\mathbf{r}_{1}(0),\cdots,\mathbf{r}_{N}(0)\\}$.

In our numerical calculations, $P$ is always finite. In this case, the
approximate $Z_{A}(\beta)$ can be written as

$Z_{A}(\beta)=\left(\frac{mP}{2\pi\hbar^{2}\beta}\right)^{PdN/2}\int
e^{-\beta(V_{B}^{(N)}+\frac{1}{P}U)}e^{-i\theta\Phi}d\mathbf{R}_{1}...d\mathbf{R}_{N},$
(13)

where $\mathbf{R}_{i}$ represents the collection of ring polymer coordinates
$(\mathbf{r}_{i}^{1},...,\mathbf{r}_{i}^{P})$ corresponding to the $i$th
particle. The system considered has $d$ spatial dimensions. For anyons
considered in this work, $d=2$. $V_{B}^{(N)}$ considers exchange effects by
describing all the possible ring polymer configurations. $U$ is the
interaction between different particles, which is given by

$U=\sum_{l=1}^{P}V(\mathbf{r}_{1}^{l},...,\mathbf{r}_{N}^{l}).$ (14)

In addition,

$-\theta\Phi(\mathbf{R}_{1},...,\mathbf{R}_{N})=\sum_{j=1}^{N}\sum_{l=1}^{P-1}(\mathbf{r}_{j}^{l+1}-\mathbf{r}_{j}^{l})\cdot\mathbf{a}_{j}(\mathbf{r}_{1}^{l},...,\mathbf{r}_{N}^{l}),$
(15)

where $\mathbf{a}_{j}(\mathbf{r}_{1}^{l},...,\mathbf{r}_{N}^{l})$ is given by

$a_{jx}(\mathbf{r}_{1}^{l},...,\mathbf{r}_{N}^{l})=\nu\sum_{k\neq
j}\frac{r_{jy}^{l}-r_{ky}^{l}}{(r_{jx}^{l}-r_{kx}^{l})^{2}+(r_{jy}^{l}-r_{ky}^{l})^{2}},$
$a_{jy}(\mathbf{r}_{1}^{l},...,\mathbf{r}_{N}^{l})=-\nu\sum_{k\neq
j}\frac{r_{jx}^{l}-r_{kx}^{l}}{(r_{jx}^{l}-r_{kx}^{l})^{2}+(r_{jy}^{l}-r_{ky}^{l})^{2}}.$
(16)

The expression of $V_{B}^{(N)}$ is the key to consider the exchange effects of
identical bosons, which will be introduced in due course.

###  II.3 Path integral molecular dynamics

To carry out the numerical calculations, we should give the expression of
$V_{B}^{(N)}$. Fortunately, in Ref. Hirshberg , it is given for path integral
molecular dynamics of bosons. In our bosonic gauge for anyons, because the
factor $\theta Q$ is symmetric about the exchange of $\mathbf{r}_{i}$ and
$\mathbf{r}_{j}$, the recursion formula of $V_{B}^{(N)}$ can be used without
any change.

The recursion formula is

$e^{-\beta
V_{B}^{(\alpha)}}=\frac{1}{\alpha}\sum_{k=1}^{\alpha}e^{-\beta(E_{\alpha}^{(k)}+V_{B}^{(\alpha-k)})},$
(17)

where $V_{B}^{(0)}=0$ and $E_{\alpha}^{(k)}$ is given by

$E_{\alpha}^{(k)}(\mathbf{R}_{\alpha-k+1},...,\mathbf{R}_{\alpha})=\frac{1}{2}m\omega_{P}^{2}\sum_{l=\alpha-k+1}^{\alpha}\sum_{j=1}^{P}(\mathbf{r}_{l}^{j+1}-\mathbf{r}_{l}^{j})^{2},$
(18)

where $\mathbf{r}_{l}^{P+1}=\mathbf{r}_{\alpha-k+1}^{1}$ if $l=\alpha$ and
$\mathbf{r}_{l}^{P+1}=\mathbf{r}_{l+1}^{1}$ otherwise. $E_{\alpha}^{(k)}$ may
be interpreted as partial configuration of particles in which the last $k$
particles are fully connected. In particular, $E_{N}^{(N)}$ represents the
longest ring configuration for $N$ particles.

From the above recursion formula, we may get

$V_{B}^{(1)}\rightarrow V_{B}^{(2)}\rightarrow\cdots
V_{B}^{(\alpha)}\cdots\rightarrow V_{B}^{(N)}.$ (19)

Therefore we see that the recursive formula allows us to compute contributions
from all ring polymer configurations correctly. Moreover, it is easy to see
that evaluating $V_{B}^{(N)}$ only takes $O(PN^{3})$ time, an exponential
speedup over simply summing all configurations.

We will use molecular dynamics and the massive and separate Nosé-Hoover chains
Nose1 ; Nose2 ; Hoover ; Martyna ; Berne ; Jang  to get the distribution
$e^{-\beta(V_{B}^{(N)}+\frac{1}{P}U)}$. To perform molecular dynamics we need
the gradient of $V_{B}^{(N)}$, which may also be computed recursively from Eq.
(17) as

$-\nabla_{\mathbf{r}_{l}^{j}}V_{B}^{(\alpha)}=-\frac{\sum_{k=1}^{\alpha}\left[\nabla_{\mathbf{r}_{l}^{j}}E_{\alpha}^{(k)}+\nabla_{\mathbf{r}_{l}^{j}}V_{B}^{(\alpha-k)}\right]e^{-\beta(E_{\alpha}^{(k)}+V_{B}^{(\alpha-k)})}}{\sum_{k=1}^{\alpha}e^{-\beta(E_{\alpha}^{(k)}+V_{B}^{(\alpha-k)})}},$
(20)

where $\nabla_{\mathbf{r}_{l}^{j}}E_{\alpha}^{(k)}$ is

$\nabla_{\mathbf{r}_{l}^{j}}E_{\alpha}^{(k)}=m\omega_{P}^{2}\left[(\mathbf{r}_{l}^{j}-\mathbf{r}_{l}^{j+1})+(\mathbf{r}_{l}^{j}-\mathbf{r}_{l}^{j-1})\right],$
(21)

with the boundary conditions
$\mathbf{r}_{l}^{P+1}=\mathbf{r}_{\alpha-k+1}^{1}$ if $l=\alpha$ and
$\mathbf{r}_{l}^{P+1}=\mathbf{r}_{l+1}^{1}$ otherwise;
$\mathbf{r}_{l}^{0}=\mathbf{r}_{\alpha}^{P}$ if $l=\alpha-k+1$ and
$\mathbf{r}_{l}^{0}=\mathbf{r}_{l-1}^{P}$ otherwise. If $l$ is outside the
interval $[\alpha-k+1,\alpha]$ then
$\nabla_{\mathbf{r}_{l}^{j}}E_{\alpha}^{(k)}=0$. By recursion calculation from
$\alpha=1$ to $\alpha=N$, we may get $-\nabla_{\mathbf{r}_{l}^{j}}V_{B}^{(N)}$
for the numerical calculations of the molecular dynamics.

The above equations completely define a molecular dynamics algorithm for
particles which scales as $O(PN^{3})$, enabling the application of PIMD to
large particle system. In order to extract physical quantities from our
simulations we use various estimators. For example, from the following
equation

$\left<E\right>=-\frac{1}{Z_{A}{(\beta)}}\frac{\partial
Z_{A}{(\beta)}}{\partial\beta},$ (22)

the average energy is given by

$\left<E\right>=\frac{PdN}{2\beta}+\frac{\left<U\right>}{P}+\left<V_{B}^{(N)}+\beta\frac{\partial
V_{B}^{(N)}}{\partial\beta}\right>.$ (23)

$V_{B}^{(N)}+\beta\frac{\partial V_{B}^{(N)}}{\partial\beta}$ may be evaluated
as

$V_{B}^{(N)}+\beta\frac{\partial
V_{B}^{(N)}}{\partial\beta}=\frac{\sum_{k=1}^{N}[V_{B}^{(N-k)}+\beta\frac{\partial
V_{B}^{(N-k)}}{\partial\beta}-E_{N}^{(k)}]e^{-\beta(E_{N}^{(k)}+V_{B}^{(N-k)})}}{\sum_{k=1}^{N}e^{-\beta(E_{N}^{(k)}+V_{B}^{(N-k)})}},$
(24)

with $V_{B}^{(0)}+\beta\frac{\partial V_{B}^{(0)}}{\partial\beta}=0$.

The average density is simply given by

$\rho(\textbf{r})=\left<\frac{1}{P}\sum_{j=1}^{P}\sum_{k=1}^{N}\delta(\mathbf{r}_{l}^{j}-\textbf{r})\right>.$
(25)

It is worthy to point out that the statistical interaction term $e^{-i\theta
Q}$ will not contribute to the expression of the energy estimator and density
estimator. The key reason why $e^{-i\theta Q}$ does not affect the energy
estimator is that the magnetic field in the region beyond the point where all
particles are located is $0$.

Nevertheless, $e^{-i\theta Q}$ will influence the average energy and average
density. In our molecular dynamics simulation, we only get the distribution of
$e^{-\beta(V_{B}^{(N)}+\frac{1}{P}U)}$. In this case, assume that
$<>^{\prime}$ denotes the average value in terms of this distribution, for an
observable $\hat{A}$, we have

$\left<\hat{A}\right>=\frac{\left<\epsilon_{A}e^{-i\theta
Q}\right>^{\prime}}{\left<e^{-i\theta Q}\right>^{\prime}}.$ (26)

Here $\epsilon_{A}$ represents the estimator for $\hat{A}$.

Because the N-particle Hamiltonian operator is hermitian, for another
hermitian operator $\hat{A}$, the average value $\left<\hat{A}\right>$ should
be real. This provides the general criteria to judge whether the numerical
result is reasonable. Of course, for $\theta\neq 0$, we will have sign
problem. Because $\theta=\pi$ is the situation of fermions, this means the
famous fermion sign problem ceperley ; loh ; troyer ; lyubartsev ; vozn ; Yao1
; Yao2 ; Yao3 . Because $\theta=0$ is the situation of bosons, this means that
bosons will not experience sign problem. The sign problem will become more
severe with the increasing of $\theta$ between $0$ and $\pi$. Compared with
the result of Ref. Hirshberg  for bosons and Ref. HirshbergFermi  for
fermions, here we give a general formula for anyons, and provide the chance to
study the continuous transition from bosons to fermions.

##  III Thermodynamics and sign problem in PIMD for anyon model

In all our importance sampling, we use massive Nosé-Hoover chain Nose1 ; Nose2
; Hoover ; Martyna ; Berne ; Jang  to establish constant temperature for our
simulations, where each separate degree of freedom of our system has been
coupled to a Nosé-Hoover thermostat. Specifically, we use the number of beads
$P=12$ and $2\times 10^{7}$ MD steps to achieve convergence. In all of the
following we checked convergence with respect to the number of beads and MD
steps performed. For details of how to assure the convergence and the general
strategy to choose the number of beads and MD steps, one may refer to the
supplementary material in Hirshberg  and also our previous works Xiong ;
Xiong1 .

Without the loss of generality, we consider the following interaction
potential

$V(\mathbf{r}_{1},\cdots,\mathbf{r}_{N})=\frac{1}{2}m\omega^{2}\sum_{j=1}^{N}|\mathbf{r}_{j}|^{2}+\sum_{l>k}^{N}\frac{\lambda}{|\mathbf{r}_{l}-\mathbf{r}_{k}|}.$
(27)

We will use the length unit $l_{0}=\sqrt{\hbar/m\omega}$ and energy unit
$E_{0}=\hbar\omega$. In this case, the dimensionless $\tilde{\beta}$ is

$\tilde{\beta}=\beta E_{0},$ (28)

while the dimensionless $\tilde{\lambda}$ is

$\tilde{\lambda}=\frac{\lambda}{E_{0}l_{0}}.$ (29)

We have

$\frac{V}{E_{0}}=\frac{1}{2}\sum_{j=1}^{N}|\tilde{\mathbf{r}}_{j}|^{2}+\sum_{l>k}^{N}\frac{\tilde{\lambda}}{|\tilde{\mathbf{r}}_{l}-\tilde{\mathbf{r}}_{k}|},$
(30)

with $\tilde{\mathbf{r}}_{k}=\mathbf{r}_{k}/l_{0}$.

Figure 2: Plot of the energy as a function of $\nu$ for 4 particles.
$\tilde{\lambda}$ is equal to 0.5. The simulation is performed at a
temperature of $\tilde{\beta}=1$. We used $P=12$ beads and $2\times 10^{7}$ MD
steps to assure convergence. The statistical fluctuations are negligible for
independent calculations. Hence, we did not give error bar in all our figures.

As a first test, we apply our method for general anyons by varying the
parameter $\nu$ from 0 (bosons) to 1 (fermions), and plot the energy as a
function of $\nu$, as shown in Fig. 2 for $N=4$ particles. We used a nonzero
interaction between particles to reduce sign problem HirshbergFermi ; Dornheim
. The energies calculated are taken to be the real part of the complex average
energy estimator presented above. The imaginary part of the energy estimator
provides an indication for the accuracy of our simulations, a small imaginary
part implies that our energy estimation is accurate. In all our simulations in
this work, the imaginary part of the energy estimator is negligible. In our
simulations, $|{ImE/ReE}|<0.002$. At $\nu=1$ for fermions, our result agrees
with the numerical result in Ref. HirshbergFermi  based on PIMD. Hence, our
simulations are correct in the whole region from bosons to fermions.

Figure 3: Plot of the magnitude of $\left<e^{-i\theta Q}\right>^{\prime}$ as a
function of $\nu$ for 4 particles. For bosons ($\nu=0$) $\left<e^{-i\theta
Q}\right>^{\prime}$ is always 1.0. As we increase $\nu$, the magnitude of
$\left<e^{-i\theta Q}\right>^{\prime}$ continuously drops to a small number,
fermions ($\nu=1$) correspond to the smallest magnitude and most severe sign
problem. The inset for $\ln\left<e^{-i\theta Q}\right>^{\prime}$ shows that
$\left<e^{-i\theta Q}\right>^{\prime}$ decreases exponentially with increasing
$v$.

In order to provide a further investigation into sign problem, we also plot
the magnitude of $\left<e^{-i\theta Q}\right>^{\prime}$ as a function of
$\nu$, as shown in Fig. 3. A smaller magnitude means that the sign problem is
more severe, and more MD iterations are needed to achieve convergence. We
found that the sign problem is especially severe for larger number of
particles and lower temperatures. To assure the accuracy, we used $2\times
10^{7}$ MD steps for importance sampling, and more MD steps were used to check
the convergence.

Figure 4: Distribution of the angle $\theta Q$, from $-\pi$ to $\pi$, for
interacting fermions with $\tilde{\lambda}=0.5$.

We also plot the distribution of the angle $\theta Q$. For fermions with non-
zero interaction, the distribution takes the form of a symmetric and Gaussian
function, as illustrated in Fig. 4. In fact, the symmetric probability
distribution of $\theta Q$ is the inevitable outcome of the path integral
formalism of the partition function given by Eq. (10).

Figure 5: The density function for anyons, normalized to unity. Blue, orange
and green represent $\nu=0$, $\nu=0.5$ and $\nu=1.0$ respectively. As $\nu$
increases the peak of the density decreases and the distribution of anyons
becomes broader.

Finally, we also investigate the density function for anyons. As $\nu$
increases from 0 to 1, the energy increases, and therefore the density should
become broader because of Pauli exclusion principle. This phenomenon is shown
in Fig. 5.

##  IV Calculation of the energy of fermions by anyon model

For anyon systems, it is clear that the sign problem becomes more severe with
the increasing of $\nu$ between $0\leq\nu\leq 1$. In Fig. 3, it is shown
clearly that $\left<e^{-i\theta Q}\right>^{\prime}$ decreases exponentially
with increasing $v$, which leads to the sign problem ceperley ; loh ; troyer ;
lyubartsev ; vozn ; Yao1 ; Yao2 ; Yao3 . In our numerical simulations, we can
calculate more accurately the energy for $0\leq\nu\leq 0.5$, while less
accurately for $\nu\geq 1/2$ (in particular $\nu=1$ for fermions) because of
the sign problem. Another numerical inaccuracy in our simulations lies in the
approximate calculation of $Q$. Nevertheless, it is expected that $E(\nu)$ is
an analytical function of $\nu$ for finite number of particles. When numerical
simulation is considered, the analyticity of $E(\nu)$ may require finite
temperature so that $P$ is finite. In the thermodynamic limit or zero
temperature with infinite $P$, because there are infinite number of beads, we
can not assure the analyticity of $E(\nu)$. Of course, it is not necessary to
worry about this in practical numerical simulation with finite number of
beads. Hence, a reasonable fitting of $E(\nu)$ for small $\nu$ may give us
accurate result for the energy of fermions $E(\nu=1)$.

From the partition function for anyons, it is quite interesting to notice that
$E(\nu)$ is a periodic function Myrheim  of $\nu$ with period $2$. From the
symmetry analysis, we also have $E(\nu)=E(-\nu)$. In addition, it is a
monotonically increasing function of $\nu$ for $0\leq\nu\leq 1$ because of
different exchange forces between bosons and fermions. The exchange forces for
fermions are repulsive, while the exchange forces for bosons are attractive.
This suggests the following simple function for $E(\nu)$.

$E(\nu)=E(\nu=0)+a\left[\sin\left(\frac{\pi\nu}{2}\right)\right]^{b}.$ (31)

Here $a$ and $b$ are two positive real numbers for fitting. The choice of
$\sin\left(\frac{\pi\nu}{2}\right)$ in the above expression is due to the
general consideration of $E(\nu)$.

Because small value of $\nu$ has less severe sign problem, we may consider to
calculate $E(\nu)$ below a value $\nu_{c}$, and use the above function to fit
the data to get $E(\nu)$ for fermions with $\nu=1$. For the same parameters in
Fig. 2, we use the numerical results of $E(0)$, $E(0.1)$, $E(0.2)$, $E(0.3)$,
$E(0.4)$, $E(0.5)$ for data fitting. In Fig. 6(a), the solid line is the best
fitting with function (31). We do find that the fitting can give accurate
energy for fermions, while suffering from much less severe sign problem. The
error of the fitting for the fermions is $2\%$, compared with the result of
direct numerical simulation. If the whole region is considered, the mean
deviation is smaller than $1\%$.

Figure 6: (a) For 4 interacting anyons with the same parameters of Fig. 2, the
solid line illustrates the fitting with the data shown by black circles. The
blue circle gives the numerical result of $E(\nu=1)$ for fermions, which
agrees with the extrapolation in our fitting. (b) For two ideal anyons in a
harmonic trap, the solid line shows the fitting of the data shown by black
circles in our numerical simulation. As a comparison, the blue circle is the
result of numerical simulation, while the red circle is the exact result. We
see that the fitting to get the energy of fermions agrees better with the red
circle. The slight difference is due to the numerical error in calculating
$e^{-i\theta Q}$.

To test further our idea, we also consider here two ideal anyons which have
analytical result for the partition function. We consider here two anyons in a
two-dimensional harmonic trap $\frac{1}{2}m\omega^{2}(x^{2}+y^{2})$. In unit
of $\hbar\omega$ and with the parameter $\tilde{\beta}=1$, in Fig. 6(b), the
black circles give the average energy $E(\nu)$ as a function of $\nu$. As
expected, we see a monotonic increase of the energy from $\nu=0$ for bosons to
$\nu=1/2$ for anyons. For this anyon model, we have analytical result of the
energy spectrum. The analytical result of the partition function is Myrheim

$Z(\beta,\nu)=\frac{\cosh((1-\nu)\beta)}{8\sinh^{2}(\frac{\beta}{2})\sinh^{2}\beta}.$
(32)

The solid line in Fig. 6(b) gives the average energy based on the fitting with
function (31). The red circle gives the energy of two fermions based on the
analytical result of the partition function, which agrees well with the
extrapolation of our fitting. The error of the fitting for the fermions is
$1\%$, compared the exact analytical result. If the whole region is
considered, the mean deviation is smaller than $0.5\%$. The slight difference
is due to the systematic error of the numerical error of $Q$ because of the
discretization in imaginary time and the function for fitting. In the limit of
zero temperature, the derivative of $E(\nu)$ about $v$ is nonzero Myrheim  at
$v=0$ and $v=1$, which is different from the situation of PIMD simulation with
finite number of beads. Hence, we should consider carefully the validity of
PIMD for zero temperature, when the practical application of our method is
considered.

The above two examples show that we have hope to calculate the energy of
fermions by anyon models with less severe sign problem. Of course, the present
idea can not be applied to large fermion systems because in this case, there
are still insurmountable obstacles due to the sign problem even for small
$\nu$.

##  V summary and discussion

In this work we showed how to extend traditional PIMD simulations chandler ;
Parrinello ; Tuckerman ; Hirshberg ; HirshbergFermi ; Miura ; Cao ; Cao2 ;
Jang2 ; Ram ; Poly ; Craig ; Braa ; Haber ; Thomas  to enable the studies of
general anyons. Through simple modifications of various physical estimators,
we are able to extract physical quantities from simulations for a system of
identical interacting anyons. Our methodology may be easily extended to study
realistic physical systems of anyons with more complex interactions. We look
forward to the applications of PIMD for anyons in the exciting field of many-
body quantum systems.

By developing a numerical algorithm for anyons, we also provide further
insights into the fermion sign problem. Using the extra freedom in choosing
the statistics parameter $\nu$, combined with the choice of interaction
parameters, we may extrapolate regime in which the sign problem is less severe
into regime with most severe sign problem HirshbergFermi ; DornheimMod  (i.e.
weakly interacting fermions), hopefully attenuating the fermion sign problem
further. Although the present method based on anyon models can not be applied
to large fermion system, our work suggests that we may consider to construct
more appropriate fictitious particle XiongSFP  with a parameter to connect
bosons and fermions, similarly to anyon models connecting continuously bosons
and fermions.

Acknowledgements This work is partly supported by the National Natural Science
Foundation of China under grants number 11175246, and 11334001.

DATA AVAILABILITY

The data that support the findings of this study are available from the
corresponding author upon reasonable request. The code of this study is openly
available in GitHub (https://github.com/xiongyunuo/PIMD-Anyon).

## References

## References

  * (1) J. M. Kosterlitz and D. J. Thouless, Metastability and Phase Transitions in Two Dimensional Systems, J. Phys. C 6, 1181 (1973).
  * (2) Z. Hadzibabic, P. Krüger, M. Cheneau, B. Battelier, and J. Dalibard, Berezinskii-Kosterlitz-Thouless crossover in a trapped atomic gas, Nature 41, 1118 (2006).
  * (3) J. M. Leinaas and J. Myrheim, On the theory of identical particles, Il Nuovo Cimento B. 37, 1 (1977).
  * (4) F. Wilczek, Magnetic Flux, Angular Momentum, and Statistics, Phys. Rev. Lett. 48, 1144 (1982).
  * (5) F. Wilczek, Quantum Mechanics of Fractional-Spin Particles, Phys. Rev. Lett. 49, 957 (1982).
  * (6) F. Wilczek, Fractional Statistics and Anyon Superconductivity (World Scientific, Teaneck, 1990).
  * (7) A. Khare, Fractional Statistics and Quantum Theory (World Scientific, Singapore, 2005).
  * (8) D. C. Tsui, H. L. Stormer, and A. C. Gossard, Two-Dimensional Magnetotransport in the Extreme Quantum Limit, Phys. Rev. Lett. 48, 1559 (1982).
  * (9) H. Bartolomei, et al., Fractional statistics in anyon collisions, Science 368, 173 (2020).
  * (10) J. Nakamura, S. Fallahi, H. Sahasrabudhe, R. Rahman, S. Liang, G. C. Gardner, and M. J. Manfra, Aharonov-Bohm interference of fractional quantum Hall edge modes, Nat. Phys. 15, 563 (2019).
  * (11) J. Nakamura, S. Liang, G. C. Gardner, and M. J. Manfra, Direct observation of anyonic braiding statistics, Nat. Phys. 16, 931 (2020).
  * (12) R. P. Feynman and A. R. Hibbs, Quantum Mechanics and Path Integrals (Dover Publications, New York, 2010).
  * (13) H. Kleinert, Path Integrals in Quantum Mechanics, Statistics, Polymer Physics, and Financial Markets (World Scientific, Singapore, 2009).
  * (14) D. M. Ceperley, Path integrals in the theory of condensed helium, Rev. Mod. Phys. 67, 279 (1995).
  * (15) M. Boninsegni, N. V. Prokofev, and B. V. Svistunov, Worm algorithm and diagrammatic Monte Carlo: A new approach to continuous-space path integral Monte Carlo simulations, Phys. Rev. E 74, 036701 (2006).
  * (16) M. Boninsegni, N. V. Prokofev, and B. V. Svistunov, Worm Algorithm for Continuous-Space Path Integral Monte Carlo Simulations, Phys. Rev. Lett. 96, 070601 (2006).
  * (17) T. Dornheim, The Fermion Sign Problem in Path Integral Monte Carlo Simulations: Quantum Dots, Ultracold Atoms, and Warm Dense Matter, Phys. Rev. E 100, 023307 (2019).
  * (18) D. Chandler and P. G. Wolynes, Exploiting the isomorphism between quantum theory and classical statistical mechanics of polyatomic fluids, J. Chem. Phys. 74, 4078 (1981).
  * (19) M. Parrinello and A. Rahman, Study of an F center in molten KCl, J. Chem. Phys. 80, 860 (1984).
  * (20) M. E. Tuckerman, Statistical Mechanics: Theory and Molecular Simulation (Oxford University, New York, 2010).
  * (21) B. Hirshberg, V. Rizzi, and M. Parrinello, Path integral molecular dynamics for bosons, Proc. Natl. Acad. Sci. U. S. A. 116, 21445 (2019).
  * (22) Y. N. Xiong and H. W. Xiong, Path integral molecular dynamics simulations for Green’s function in a system of identical bosons, J. Chem. Phys. 156, 134112 (2022).
  * (23) Y. L. Yu, S. J. Liu, H. W. Xiong, and Y. N. Xiong, Path integral molecular dynamics for thermodynamics and Green’s function of ultracold spinor bosons, arXiv:2207.07653 (2022).
  * (24) B. Hirshberg, M. Invernizzi, and M. Parrinello, Path integral molecular dynamics for fermions: Alleviating the sign problem with the Bogoliubov inequality, J. Chem. Phys. 152, 171102 (2020).
  * (25) Y. N. Xiong and H. W. Xiong, Numerical calculation of Green’s function and momentum distribution for spin-polarized fermions by path integral molecular dynamics, J. Chem. Phys. 156, 204117 (2022).
  * (26) D. M. Ceperley, Path Integral Monte Carlo Methods for Fermions, Monte Carlo and Molecular Dynamics of Condensed Matter Systems, K. Binder and G. Ciccotti (Eds.), Bologna (Italy) (1996).
  * (27) E. Y. Loh, J. E. Gubernatis, R. T. Scalettar, S. R. White, D. J. Scalapino, and R. L. Sugar, Sign problem in the numerical simulation of many-electron systems, Phys. Rev. B 41, 9301 (1990).
  * (28) M. Troyer and U. J. Wiese, Computational Complexity and Fundamental Limitations to Fermionic Quantum Monte Carlo Simulations, Phys. Rev. Lett. 94, 170201 (2005).
  * (29) A. P. Lyubartsev, Simulation of excited states and the sign problem in the path integral Monte Carlo method, J. Phys. A: Math. Gen. 38, 6659 (2005).
  * (30) M. A. Voznesenskiy, P. N. Vorontsov-Velyaminov, and A. P. Lyubartsev, Path-integral-expanded-ensemble Monte Carlo method in treatment of the sign problem for fermions, Phys. Rev. E 80, 066702 (2009).
  * (31) Z. X. Li, Y. F. Jiang, and H. Yao, Solving the fermion sign problem in quantum Monte Carlo simulations by Majorana representation, Phys. Rev. B 91, 241117(R) (2015).
  * (32) Z. X. Li, Y. F. Jiang, and H. Yao, Majorana-Time-Reversal Symmetries: A Fundamental Principle for Sign-Problem-Free Quantum Monte Carlo Simulations, Phys. Rev. Lett. 117, 267002 (2016).
  * (33) Z. X. Li and H. Yao, Sign-Problem-Free Fermionic Quantum Monte Carlo: Developments and Applications, Ann. Rev. Cond. Mat. Phys. 10, 337 (2019).
  * (34) J. Myrheim, in Anyons edited by A. Comtet, T. Jolicoeur, S. Ouvry, and F. David (Springer-Verlag, 1999), 69, pp. 265-413.
  * (35) S. Nosé, A molecular dynamics method for simulations in the canonical ensemble, Mol. Phys. 52, 255 (1984).
  * (36) S. Nosé, A unified formulation of the constant temperature molecular dynamics methods, J. Chem. Phys. 81, 511 (1984).
  * (37) W. G. Hoover, Canonical dynamics: Equilibrium phase-space distributions, Phys. Rev. A 31, 1695 (1985).
  * (38) G. J. Martyna, M. L. Klein, and M. Tuckerman, Nosé-Hoover chains: The canonical ensemble via continuous dynamics, J. Chem. Phys. 97, 2635 (1992).
  * (39) Mark E. Tuckerman and Bruce J. Berne, Efficient molecular dynamics and hybrid Monte Carlo algorithms for path integrals, J. Chem. Phys. 99, 2796 (1993).
  * (40) S. Jang and G. A. Voth, Simple reversible molecular dynamics algorithms for Nosé-Hoover chain dynamics, J. Chem. Phys. 107, 9514 (1997).
  * (41) S. Miura and S. Okazaki, Path integral molecular dynamics for Bose-Einstein and Fermi-Dirac statistics. J. Chem. Phys. 112, 10116 (2000).
  * (42) J. Cao and G. A. Voth, The formulation of quantum statistical mechanics based on the Feynman path centroid density. I. Equilibrium properties, J. Chem. Phys. 100, 5093 (1994).
  * (43) J. Cao and G. A. Voth, The formulation of quantum statistical mechanics based on the Feynman path centroid density. II. Dynamical properties, J. Chem. Phys. 100, 5106 (1994).
  * (44) S. Jang and G. A. Voth, A derivation of centroid molecular dynamics and other approximate time evolution methods for path integral centroid variables, J. Chem. Phys. 111, 2371 (1999).
  * (45) R. RamíRez and T. LóPez-Ciudad, The Schrödinger formulation of the Feynman path centroid density, J. Chem. Phys. 111, 3339 (1999).
  * (46) E. A. Polyakov, A. P. Lyubartsev, and P. N. Vorontsov-Velyaminov, Centroid molecular dynamics: Comparison with exact results for model systems, J. Chem. Phys. 133, 194103 (2010).
  * (47) I. R. Craig and D. E. Manolopoulos, Quantum statistics and classical mechanics: Real time correlation functions from ring polymer molecular dynamics, J. Chem. Phys. 121, 3368 (2004).
  * (48) B. J. Braams and D. E. Manolopoulos, On the short-time limit of ring polymer molecular dynamics, J. Chem. Phys. 125 ,124105 (2006).
  * (49) S. Habershon, D. E. Manolopoulos, T. E. Markland, and T. F. Miller 3rd, Ring-polymer molecular dynamics: quantum effects in chemical dynamics from classical trajectories in an extended phase space, Annu. Rev. Phys. Chem. 64, 387 (2013).
  * (50) T. E. Markland and M. Ceriotti, Nuclear quantum effects enter the mainstream, Nat. Rev. Chem. 2, 0109 (2018).
  * (51) T. Dornheim, M. Invernizzi, J. Vorberger, and B. Hirshber, Attenuating the fermion sign problem in path integral Monte Carlo simulations using the Bogoliubov inequality and thermodynamic integration, J. Chem. Phys. 153, 234104 (2020).
  * (52) Y. N. Xiong and H. W. Xiong, A solution of fermion sign problem for large fermion systems, arXiv:2206.08341 (2022).