# Robust Semi-Supervised Learning for Histopathology Images through Self-
Supervision Guided Out-of-Distribution Scoring

Nikhil Cherian Kurian, Varsha S, Abhijit Patil, Shashikant Khade, and Amit
Sethi NC Kurian, S Varsha, A Patil, S Khade and A Sethi are with the
Department of Electrical Engineering, Indian Institute of Technology Bombay,
Mumbai, MH, 400076 India e-mail<EMAIL_ADDRESS>
###### Abstract

Semi-supervised learning (semi-SL) is a promising alternative to supervised
learning for medical image analysis when obtaining good quality supervision
for medical imaging is difficult. However, semi-SL assumes that the underlying
distribution of unaudited data matches that of the few labeled samples, which
is often violated in practical settings, particularly in medical images. The
presence of out-of-distribution (OOD) samples in the unlabeled training pool
of semi-SL is inevitable and can reduce the efficiency of the algorithm.
Common preprocessing methods to filter out outlier samples may not be suitable
for medical images that involve a wide range of anatomical structures and rare
morphologies. In this paper, we propose a novel pipeline for addressing open-
set supervised learning challenges in digital histology images. Our pipeline
efficiently estimates an OOD score for each unlabelled data point based on
self-supervised learning to calibrate the knowledge needed for a subsequent
semi-SL framework. The outlier score derived from the OOD detector is used to
modulate sample selection for the subsequent semi-SL stage, ensuring that
samples conforming to the distribution of the few labeled samples are more
frequently exposed to the subsequent semi-SL framework. Our framework is
compatible with any semi-SL framework, and we base our experiments on the
popular Mixmatch semi-SL framework. We conduct extensive studies on two
digital pathology datasets, Kather colorectal histology dataset and a dataset
derived from TCGA-BRCA whole slide images, and establish the effectiveness of
our method by comparing with popular methods and frameworks in semi-SL
algorithms through various experiments.

###### Index Terms:

Semi Supervised learning, open-set, label-noise, mixmatch

##  I Introduction

Medical image analysis requires large volumes of supervised data to train deep
learning models effectively, but obtaining good quality supervision for
medical imaging is inherently difficult due to the associated labor,
expertise, and time required [1, 2, 3]. In such scenarios, semi-supervised
learning (semi-SL) offers an efficient alternative, especially when there are
only a few labeled samples but plenty of unlabeled or unaudited data. Semi-SL
algorithms can leverage the vast pool of unaudited training data by extracting
discriminative information from the structure of unlabeled data that
complements the knowledge gained from a small number of supervisory data
samples. However, semi-SL assumes that the underlying distribution of the
unaudited data matches that of the few labeled samples [4, 5].

The idea of the data distributions matching is frequently challenged in
practical situations, particularly when the data is sourced from datasets with
high heterogeneity. An example of such datasets’ is medical images where the
manifestation of patient characteristics and diseases can exhibit significant
diversity.[6]. Furthermore, representing disease continuums as discrete
categories inherently introduces out-of-distribution (OOD) samples in medical
imaging data. Consequently, the presence of OOD samples in the unlabelled
training pool of semi-supervised learning (semi-SL) is inevitable. In a
supervised learning framework, such OOD training samples are referred to as
open-set samples, as they may not strictly conform to the distribution of any
predefined classes of interest in a classification problem [7].

In the case of deep semi-SL frameworks, the presence of OOD samples in the
unlabeled training data can have a negative impact on knowledge propagation,
thereby reducing the efficiency of the algorithm [4]. To address this issue,
common preprocessing methods aim to filter out the outlier samples [4].
However, such methods may not be suitable for semi-SL frameworks that deal
with medical images. Medical images often involve a wide range of anatomical
structures with significant differences in their visible features, making it
difficult to eliminate samples without losing valuable information.
Additionally, certain disease conditions in medical images can present as rare
morphologies, which may be missed by filtering out samples. Therefore, it is
important to develop tailored solutions for open-set semi-supervision in
medical imaging.

The focus of our work is to address the challenges of open-set supervised
learning in digital histology images, where the availability of labeled data
is limited but noise-free with respect to its supervision and unlabeled data
contains open-set samples. We propose a pipeline that combines an OOD detector
based on self-supervised learning with a semi-SL framework to effectively
utilize the large pool of unaudited data for training. Our pipeline is
designed specifically for medical images, which often contain novel samples
that should not be naively discarded. The OOD detector’s outlier score is used
to modulate the sampling selection for the subsequent semi-SL stage to ensure
that the samples conform to the distribution of the few labeled samples.

We conducted extensive experiments on two digital pathology datasets: Kather
colorectal histology dataset [8] and a dataset derived from TCGA-BRCA whole
slide images (WSIs) [9]. The Kather dataset was customized to contain varying
proportions of open-set data in unlabeled samples. Our experiments show that
our proposed method is effective in addressing open-set supervised learning
challenges in digital histology images compared to other popular methods and
frameworks in semi-SL algorithms.

Our method is compatible with any semi-SL framework, though we base our
experiments on the popular Mixmatch semi-SL framework. We demonstrate that our
method outperforms other semi-SL frameworks in terms of accuracy,
demonstrating the effectiveness of our method in digital histology images. The
pipeline can be used as a solution to address open-set supervised learning
challenges in digital histology images where labeled data is limited and
unlabeled data contains open-set samples.

Figure 1: Proposed self-supervision based open-set semi-SL framework.

##  II Related Work

Several methods have been reported in literature that developed novel
algorithms that can leverage the supervisory information from labeled data and
propagate it into unlabelled samples. The state-of-the-art methods in semi-SL
employ consistency regularization [10] to reliably extract and incorporate
unsupervised information coherent with the labeled data. In consistency
regularization [10], the classification results for the unlabeled sample are
regularized to remain unchanged with data augmentations ensuring that only
robust information encoded in the structure of unlabeled data is learned. The
major semi-SL algorithms that use this consistency regularization [10] are
enlisted and concisely explained below. We also summarize a multi-task
curriculum learning framework (MTL) [4] approach, which has recently been
reported to successfully address the open-set issues in semi-supervised
learning.

####  II-1 MixMatch

MixMatch [11] is a popular semi-SL algorithm that integrated three distinct
algorithmic subroutines: consistency regularization [10], entropy
minimization, and MixUp [12], to ensure efficient training and knowledge
abstraction from the labeled and unlabeled training samples. As mentioned
before, consistency regularization [10] is employed in these algorithms to
ensure that training data (including the unlabeled data) and its augmentations
have similar predictions throughout the training. Entropy minimization
sharpens the pseudo-labeled predictions on multiple input augmentations,
encouraging the model to output more confident predictions on the unlabelled
samples. Finally, by augmenting the input data with convex combinations of the
training samples and their corresponding labels, the concept of MixUp [12]
introduces a linear relationship between the training samples. The entire
framework uses a cross-entropy based loss function for the fully supervised
samples, whereas a mean square loss is adopted to optimize the unlabeled data.

####  II-2 FixMatch

FixMatch [13] is another semi-SL algorithm that works on the principles of
pseudo-labeling and consistency regularization [10] concepts. FixMatch [13]
modifies the original MixMatch [11] framework by removing the MixUp [12] stage
and introducing an enhanced version of consistency regularization [10] between
weakly and strongly augmented versions of training samples. FixMatch [13]
ensures superior pseudo-labeling by imposing an additional threshold value on
weak augmentation so that only highly confident predictions from these samples
are considered for consistency regularization.

####  II-3 ReMixMatch

ReMixMatch [14] enhances the MixMatch [11] framework by introducing the
principles of distribution alignment and augmentation anchoring. The notion of
distribution alignment replaces the sharpening step of the MixMatch [11]
framework in an attempt to match model aggregated class prediction to that of
the marginal distribution of the given ground truth. Further, ReMixMatch [14]
introduces another principle of augmentation alignment in place of the
consistency regularization step of MixMatch [11] to encourage each output to
be close to the prediction for a weakly-augmented version of the same input.
By incorporating these changes, ReMixMatch [14] has been reported to be data
efficient compared to MixMatch [11] frameworks.

####  II-4 Multi-task Curriculum learning (MTL) Framework Guided Semi-SL

Unlike previous methods that focused on developing strong semi-SL models, the
goal of MTL [4] guided semi-SL is to explicitly address the issues associated
with the presence of open-set samples in semi-SL. The proposed framework
addresses this issue by simultaneously training and optimizing OOD detection,
and a MixMatch [11]-based semi-SL in a multitask learning framework. The paper
proposes a novel OOD detector based on the model’s capacity to identify noisy
labeled training data. The OOD detector ensures that the subsequent semi-SL
framework, based on MixMatch [11], is only trained with the inliers samples by
selecting an appropriate threshold on the OOD score and filtering out the
outlier samples from the unlabelled data. The method significantly produced
superior results in settings where samples from a different distribution
contaminates the unlabelled data.

Along with these frameworks, a few studies have also re-examined the original
MixMatch [11] algorithm itself. One such study found that the performance
degradation caused by open-set samples in unlabelled data is primarily due to
the Pseudo-Labelling (PL) task of MixMatch. Additionally, they also report
that the incorporation of consistency regularization has the potential to
enhance the performance of MixMatch in semi-SL, even when unlabelled training
samples contain open-set data.

##  III Method

In this section, we present our robust semi-supervised learning (semi-SL)
framework which involves multiple stages of deep learning-oriented data
preprocessing. Our approach is designed to limit the exposure of open-set
samples to a MixMatch-based semi-supervised learning framework.

The first stage involves training a deep learning (DL) architecture based on
self-supervised learning to create a novel outlier detector. In the second
stage, we calculate an out-of-distribution (OOD) score based on Gaussian
mixture modeling. Finally, we develop a custom data sampler that modifies the
exposure of unlabeled samples based on the OOD scores we calculated.

The entire methodology is summarized in Figure 1. In the following sections,
we provide a detailed explanation of each of these steps.

####  III-1 Self-supervised learning stage

To develop our OOD detector, we utilized self-supervised learning (SSL) which
has been reported to discover the underlying structure of the data without
explicit supervision [15]. To achieve this, we employed the widely used
SimCLR-based self-supervised learning paradigm [16] and tailored it to train a
novel OOD detector for our purposes.

The SimCLR framework [16] is based on contrastive learning principles where a
sample and its augmented version are trained to produce similar embeddings in
a high-dimensional latent space of dimension $D$ by optimizing a loss function
based on cosine similarity. In order to train such networks, each mini-batch
must include multiple augmented views of a training sample, and each training
sample must serve as a positive anchor for all of its augmented versions
during loss calculations.

In our methodology, in addition to the augmented versions of a training
sample, we also encourage positive anchors to attract any samples that lie
within an $\epsilon$ neighbourhood in the high-dimensional embedding space
based on their cosine similarity. This ensures that samples with similar
underlying structures are clustered closer together at the end of training.

This means that a sample $x_{j}$ is included as a positive pair with an anchor
sample $x_{i}$ if it is an augmented version of $x_{i}$, or if the cosine
similarity of their latent space embeddings $z_{i}$ and $z_{j}$ are similar to
each other based on Equation 1.

$sim(z_{i},z_{j})\geq 1-\epsilon$ (1)

Here $sim$ refers to the cosine similarity function as defined below in
Equation 2.

$sim(z_{i},z_{j})=\frac{z_{i}^{T}z_{j}}{\left|z_{i}\right|_{2}\left|z_{j}\right|_{2}}$
(2)

where $T$, represents the vector transpose operation and the notation
$\left|z_{i}\right|_{2}$ represents the $L_{2}$ norm of the vector $z_{i}$.

For a minibatch of size $2N$ ($N$ samples collated with their augmented
versions), the loss function $\mathds{L}(z_{i},z_{j})$, for the latent space
embeddings $(z_{i},z_{j})$ corresponding to a positive pair of training
samples $(x_{i},x_{j})$ is given by the following negative log likelihood
formulation over softmax as described in Equation 3

$\mathds{L}(z_{i},z_{j})=-log\frac{exp(sim(z_{i},z_{j})/\tau)}{\sum_{k=1}^{2N}\mathbbm{1}_{[k\neq
i]}exp(sim(z_{i},z_{k})/\tau)}$ (3)

In this equation, $\mathbbm{1}_{[k\neq i]}\in\\{0,1\\}$ is an indicator
function evaluating to 1 iff $k\neq i$, $\tau$ denotes a temperature parameter
and $k$ denotes the index of all samples that forms the negative pair to the
anchor samples index $i$.

####  III-2 GMM-based outlier score calculation

Once we have trained the self-supervised OOD detector model, we proceed to
estimate outlier scores for each unlabeled training sample by modelling latent
space as a mixture of Gaussians[17]. The Gaussian mixture model is a
probabilistic clustering model that assumes all data points are generated from
a mixture of a finite number of Gaussian distributions with unknown parameters
[17].

In this case, for the number of cluster components $N_{cls}$, we assume that
for each cluster $c_{n}$, the distribution of the $D$ dimensional latent space
embedding variable $z$ is having a marginal probability distribution of the
form:

$P(z)=\sum_{n=1}^{N_{cls}}\pi_{n}P(z|C_{i}=c_{n})$ (4)

where $i\in\\{1,2,\cdots N_{cls}\\}$, $P(z|C_{i}=c_{n})$ is the likelihood of
the mixture component and $\pi_{n}$ is the prior probability associated with
the cluster component $c_{n}$.

In GMMs, the likelihood distribution $P(z|C_{i}=c_{n})$ is further assumed to
be a multivariate $D$ dimensional normal distribution,
$\mathcal{N}(z;\mu_{n},\Sigma_{n})$ with mean vector $\mu_{n}$ and covariance
matrix $\Sigma_{n}$.

Finally, we use an iterative Expectation-Maximization (EM) algorithm [18] to
estimate the $\mu_{n}$, $\Sigma_{n}$ and $\pi_{n}$ parameters of the GMM. The
GMM-based soft clustering provides a normalized soft-cluster membership score
of $\phi_{k}(c_{j})$ in the range of $(0,1)$ for the cluster component $c_{j}$
as a posterior distribution from the prior $\pi_{n}$ and the likelihood,
$\mathcal{N}(z;\mu_{n},\Sigma_{n})$, as defined in Equation 5

$\phi_{k}(c_{j})=\frac{\pi_{j}\mathcal{N}(z_{k};\mu_{j},\Sigma_{j})}{\sum_{n=1}^{N_{cls}}\pi_{n}\mathcal{N}(z_{k};\mu_{n},\Sigma_{n})}$
(5)

To estimate a cluster impurity score (CIS), $CIS(c_{j})$, for the GMM cluster
component $c_{j}$, we use the cumulative GMM membership of all labeled
training samples to the total GMM membership of all the samples in a given
cluster, as specified in Equation 6. The intuition behind this approach is
that a cluster is likely to have more inlier samples if they have more GMM
membership affiliations from the inlier labeled training samples.

$CIS(c_{j})=-\log\frac{\sum_{l=1}^{|\mathcal{L}|}\phi_{\mathcal{L}_{l}}(c_{j})}{\sum_{l=1}^{|\mathcal{L}|}\phi_{\mathcal{L}_{l}}(c_{j})+\sum_{u=1}^{|\mathcal{U}|}\phi_{\mathcal{U}_{u}}(c_{j})}$
(6)

In Equation 6, $\mathcal{L}$ and $\mathcal{U}$ represents the labeled and
unlabeled subsets of the training dataset. Finally, OOD score for an unlabeled
sample, $x_{k}$, is given by the weighted sum of cluster impurity scores with
sample membership scores over all the $N_{cls}$ as given in Equation 7

$OOD\text{ }Score(x_{k})=\sum_{j=1}^{N_{cls}}CIS(c_{j})\times\phi_{k}(c_{j})$
(7)

####  III-3 Design of OOD Score-Driven Data Sampler and MixMatch

To prevent outlier samples from affecting the performance of the MixMatch-
based semi-supervised learning framework, we employ the OOD scores computed
for each unlabeled sample to guide the data selection process in each mini-
batch. Specifically, we create a data sampler that prioritizes or samples more
from the inlier regions of the unlabeled training data. Below are the steps to
design this sampler:

  1. 1:

We cluster the unlabeled samples into $N_{cls}$ using a GMM-based approach, as
described before. Afterwards, we compute a cluster impurity score (CIS) using
Equation 6. To avoid having clusters with the same weights, we combine them to
form super clusters, resulting in $N$ clusters with unique CIS values.

  2. 2:

We use a weighted sampler to select one of the $N$ clusters defined in step 1.
The weighting factor is the inverse of the impurity score of each score, as
this ensures that more samples will be chosen from the less impure clusters.

  3. 3:

From the cluster selected in step 2, we sample an unlabeled data sample using
weights that are inversely proportional to the OOD scores as defined in
Equation 7.

We repeat steps 1 and 2 for all the samples in a mini-batch. For the final
semi-supervised learning stage, we use the MixMatch framework. The MixMatch
framework consists of two components in its loss function: a cross-entropy
loss for the labeled training data, as given by Equation 8, and a consistency
loss for the unlabeled data, as given by Equation 9,

$\mathds{L}_{\mathcal{L}}=mean_{|\mathcal{L}|}(-\sum_{c=1}^{C}p_{sharp_{c}}\log(p_{model_{c}}(y_{li}/x_{li})))$
(8)

where, $p_{sharp}$ is the sharpened prediction of labeled sample over $k$
augmentations of the input $x_{li}$ with labels $y_{li}$, $C$ denotes the
total number of classes and $p_{model}$ denotes the predictions from the
model.

The unlabeled loss is mean-square error defined in the Equation 9,

$\mathds{L}_{\mathcal{U}}=mean_{|\mathcal{U}|}(-||q-p_{model}(x_{u})||_{2})$
(9)

where, $q$ is the pseudo label guessed after sharpening operation by the model
over $k$ augmentations, $p_{models}$ is the model’s raw prediction over one
unlabeled data, $x_{u}$ denotes the unlabeled input. Relating to the MixMatch
[11] theory mentioned in section II-1, as per the Equations 8 and 9, the
sharpening operation incorporates entropy minimization and the loss
minimization ensures consistency regularization.

##  IV Experiments

Our study involved a series of experiments to evaluate the effectiveness of
our proposed method in comparison to existing semi-supervised learning
frameworks when the unlabeled training data contained open-set samples.
Additionally, we conducted experiments to determine optimal design parameters
for our approach.

For our experiments, we utilized two datasets primarily. The first dataset,
the Kather colorectal cancer (CRC) histopathology dataset [8], was publicly
available. To simulate an open-set semi-SL problem, we reserved a portion of
the classes in the dataset as OODs, which were solely present in the unlabeled
training data. To assess the robustness of our algorithm, we varied the amount
of OODs in the unlabeled training data. Our second dataset, the TCGA-BRCA
cohort [9], was intended to demonstrate the effectiveness of our method in
scenarios where open-set samples naturally exist in large unaudited data
pools.

In all of our experiments, we utilized an OOD detector based on a SimCLR
architecture using ResNet18 [16]. We made adjustments to the ResNet18 model by
eliminating the fully-connected layers after the global average pooling layer
and appended a linear layer with dimensions of $128\times 1$. The latent space
embeddings that perform the SimCLR self-supervision [16] were obtained from
this linear layer. Additionally, this ResNet model [19] had already undergone
pre-training on a vast array of histology images utilizing SimCLR-based self-
supervision [15], making it ideal for fine-tuning on our dataset.

For all our experiments, we selected $\epsilon=0.05$ based on empirical
observations. We employed the MixMatch framework for semi-supervised learning,
utilizing a wide-ResNet architecture [20]. Our data augmentations included
color jitter, random horizontal and vertical flips, as well as elastic
deformations.

Figure 2: Cluster wise distribution of inlier or OOD samples vs purity of
cluster in the experiment on Kather dataset for the best and worst case
settings in our experiments. We can see that clusters with higher CIS have
more outlier samples.

###  IV-A Experiments with colorectal cancer (CRC) Dataset

The Kather colorectal dataset [8] comprises 100,000 images of colorectal
histology patches from H&E-stained WSIs. These patches include nine types of
tissue regions, such as adipose, background, debris, lymphocytes, mucus,
normal colon mucosa, smooth muscles, cancer-associated stroma, and tumor
epithelial regions. To evaluate open-set semi-supervised learning performance,
we reorganized this dataset and constructed a four-class classification
problem, using lymphocytes, normal colon mucosa, cancer-associated stroma, and
tumor epithelial regions as the major classes. We tagged the remaining five
classes as open-set samples.

To conform with the typical settings of the semi-supervised learning problem,
we used only a small number of labeled samples from the four major classes for
our experiments. The training data consisted mainly of unlabeled data sampled
from both inlier and open-set classes. We manually varied the number of open-
set samples in the unlabeled data to assess the performance drop caused by
these samples. We used a labeled training data size of either 25 or 50 in each
of our experiments, while keeping the number of unlabeled inlier data constant
at 10,000. Additionally, we used three different sample sizes for the open-set
unlabeled data: 8000, 10000, and 15000 samples.

We held out a class-balanced testing set of 4800 samples for final evaluations
and also used a validation data size of 400 in our experiments.

The DNN models, including the OOD detector and MixMatch semi-SL framework,
were trained using the Adam optimizer with a learning rate of 0.0003. The
batch size for the OOD detector and MixMatch models was set to 32 and 128,
respectively. In addition, we used 12 clusters in GMM, which was determined
based on additional experiments.

###  IV-B Experiments with TCGA-BRCA Whole Slide Images

Our second set of experiments focused on classifying H&E stained whole slide
images (WSIs) based on their intrinsic subtype. We choose a Basal versus
Luminal A PAM50 subtype classification for this experiment as these represent
the most common histological subtype representing the worst case and best case
prognosis respectively. The experiment was performed on patches extracted from
WSIs and evaluated at a WSI level through patch-wise aggregation. The patches
were extracted from regions of interest (tumor regions) annotated by a
pathologist, which led to label corruptions due to intrinsic subtype
heterogeneity in certain patients. As all patches inherited the same tumor
subtype label, the label noise was more of an open-set or out-of-distribution
noise, as heterogeneous tumors can exhibit a wide range of spatial cues.

The genomic criteria used to filter WSI cases for the experiment was based on
a semi-supervised non-negative matrix factorization (SS-NMF) of PAM 50 genes
detailed in [21]. Further considering quality control factors we choose a
subset of 180 WSIs from the TCGA-BRCA dataset [9]. This dataset was then
divided into train, validation, and test sets. We selected samples with the
least heterogeneity for the labeled part of the training data and for the
held-out validation and test sets. As is typical in semi-SL frameworks, the
number of training samples needed to be significantly smaller than the
unlabeled training data. Therefore, we chose 10 WSIs (5 from Luminal A and 5
from Basal) with the most pure genomic tumor subtype signature for training,
10 WSIs for validation, and 70 WSIs for testing, all of which exhibited higher
genomic tumor subtype purity compared to remaining cases. The remaining 90
WSIs, which had an admixed tumor mutational landscape, were selected as
unlabeled data. In total, we extracted 50,000 patches that were approximately
class-balanced across all data splits for the analysis.

TABLE I: Classification test accuracies on Kather dataset [8] with varied
levels of labelled data and OOD noise. Results reported here are average
accuracy over three set of runs with different random seeds.

Labelled data | 25 |  | 50
---|---|---|---
OOD Noise | 8000 | 10000 | 15000 |  | 8000 | 10000 | 15000
MixMatch [11] | 92.37$\pm$1.27 | 91.39$\pm$0.67 | 90.31$\pm$0.76 |  | 94.20$\pm$1.01 | 91.77$\pm$0.15 | 91.09$\pm$0.54
MixMatch (w/o PL) | 93.68$\pm$0.66 | 91.42$\pm$0.54 | 92.09$\pm$0.41 |  | 93.36$\pm$0.48 | 92.99$\pm$0.55 | 91.17$\pm$0.08
FixMatch [13] | 94.13$\pm$0.26 | 90.05$\pm$0.47 | 88.46$\pm$0.61 |  | 94.59$\pm$0.39 | 92.50$\pm$0.20 | 91.22$\pm$0.26
ReMixMatch [14] | 95.10$\pm$0.22 | 93.10$\pm$0.31 | 90.50$\pm$0.27 |  | 95.20$\pm$0.38 | 93.30$\pm$0.42 | 90.90$\pm$0.31
MTL [4] | 94.55$\pm$0.02 | 94.75$\pm$0.15 | 92.63$\pm$0.14 |  | 96.79$\pm$0.25 | 94.27$\pm$0.47 | 93.43$\pm$0.11
Proposed Method | 97.22$\pm$0.29 | 96.55$\pm$0.20 | 95.58$\pm$0.67 |  | 97.82$\pm$0.23 | 97.90$\pm$0.06 | 97.58$\pm$0.46

##  V Results

###  V-A Results on Kather Dataset

After completing the training of our OOD detector, we evaluate its performance
by analyzing and visualizing the distribution of GMM cluster membership for
in-class and OOD samples in each of the 12 clusters. This provides us with a
measure of the effectiveness of our cluster impurity score and the final OOD
score. As shown in Figure 2, we observe that impure clusters with higher CIS
contain more OOD samples, while pure clusters have fewer OOD samples. Figure 3
demonstrates how our sampling technique regulates the level of exposure of OOD
samples to the following MixMatch [11] framework.

Figure 3: Data sampler used in the experiments with the Kather dataset. We can
see here that this sampler chose more samples from the unlabeled training data
that had less OOD scores. Figure 4: The distribution of test samples with the
OOD scores using the proposed method in Kather dataset experiments.

Table I presents the main results of our experiments on the Kather dataset
[8]. It can be observed that as the number of open-set samples in the
unlabeled data increases, there is a significant decrease in test accuracy.
However, our proposed method demonstrates robustness to such training data
corruption compared to all other baseline methods. The MTL [4] framework also
performs better in these experiments. Moreover, we notice that the overall
performance significantly improves when more labeled samples are included in
the training data. This increase in accuracy can be attributed to the
simplicity of the Kather dataset [8] that has less intra-class variance.
Therefore, more labeled samples can effectively overcome the issues associated
with open-set unlabeled samples in the Kather dataset [8]. Additionally, we
also ran the MixMatch [11] framework in a scenario where there was no
corruption in the unlabeled data, and we report an average accuracy over three
runs to be $98.04\%\pm 0.23$ when there were 25 labelled training data and an
average accuracy of $98.04\%\pm 0.23$, when the labeled training data was 50.
Furthermore, the effectiveness of our proposed sampling method in calibrating
the exposure of OOD samples to the subsequent MixMatch [11] framework is
illustrated in Figure 3. The sampler modifies the data selection distribution
in the unlabelled training data by choosing more samples from the low OOD
score region.

Figure 5: Summary of results from experiment with RCGA-BRCA WSI dataset. Tumor
region annotation and patch extraction in (A), and patch clusters with low and
high impurity scores in (B). High cellularity patches are found in low CIS
clusters, while high CIS clusters exhibit patches with lower cellularity or
reduced label correlation due to the presence of artifacts [6]. TABLE II:
Comparison of OOD detection AUC scores on the test dataset of Kather dataset
that had 4800 samples from inlier and OOD samples.

Labelled data | 25 |  | 50
---|---|---|---
OOD Noise | 8000 | 10000 | 15000 |  | 8000 | 10000 | 15000
MTL [4] | 0.91 | 0.91 | 0.94 |  | 0.94 | 0.89 | 0.87
Proposed OOD Detector | 0.92 | 0.94 | 0.94 |  | 0.97 | 0.93 | 0.96

Figure 6: The effect of $N_{cluster}$ on the overall performance in Kather
dataset experiments for the best case scenario corresponding to 50 labeled
training samples and 8000 unlabeled open-set data, and the worst case scenario
in which only 25 labeled train samples and 15000 unlabeled open-set were used.
Based on this study, we fixed $N_{cluster}$ to 12.

Additional experiments were conducted to evaluate the efficiency of our OOD
detector. We plotted the distribution of test samples with the OOD scores in
Figure 4. To further assess the OOD detector’s performance, we tabulated the
area under the receiver operating characteristic curve in Table II. This was
done on a separate test dataset consisting of 9600 samples of both OOD and
inlier images, evenly sampled across all classes.

We also examined the impact of the hyperparameters, specifically
$N_{clusters}$, on the overall performance, as depicted in Figure 6. We
conducted these evaluations on a dataset of 9600 instances, which consisted of
4800 OOD samples added to the inlier test data.

TABLE III: Patient-level accuracy and patch-level shannon entropy (averaged on
correctly classified patients) of proposed openset SSL method for TCGA dataset
[9]

Method | Patient Level Accuracy | Patch-Level Entropy
---|---|---
MixMatch [11] | 72.86% | 0.421
MixMatch (w/o PL) | 75.71% | 0.351
FixMatch [13] | 77.14% | 0.372
ReMixMatch [14] | 75.71% | 0.396
MTL [4] | 71.43% | 0.394
Proposed Method | 81.43% | 0.342
Proposed Method (weighted voting) | 81.43% | 0.305

Additionally, it should be noted that while our test samples were chosen to
have a relatively pure genomic signature, there is still a possibility of
label noise corruption in test patches. To address this issue, we implemented
a modified weighted majority voting strategy that incorporated an inlier
scores (1- normalized OOD scores) as weights in the patient level aggregation.
However, with this modified weighting the results in Table III does not
indicate any significant improvement in patient level accuracy. On the other
hand, we did observe a decrease in entropy associated with patch-wise
prediction when we calculated the Shannon entropy of the number of patches
that were correctly classified as Luminal A and non-Luminal A. This finding
indicates higher confidence in patient level predictions when the entropy is
lower, which is a promising outcome for our approach. Moreover, we conducted
experiments where we sampled patches from GMM clusters with low and high CIS,
and our observations provide valuable insights into the effectiveness of our
approach. The results of these experiments are summarized in Figure 5.
Overall, our study contributes to the growing body of research in this field
and highlights the importance of addressing label noise corruption in genomic
signature analysis. We also sampled some of the patches present in the GMM
clusters with low and high CIS and the observations are summarized in Figure
5.

##  VI Conclusion

In this paper, we proposed a novel multi-stage framework to tackle the open-
set semi-supervised learning problem in histopathology. Our approach was
designed to address the challenges of limited labeled data and the presence of
open-set samples in the unlabeled data. Unlike other methods that use an OOD
detector to discard open-set samples from the unlabeled training data, our
approach limited the exposure of such samples to the semi-supervised learning
framework. This was done to ensure that the training data was not artificially
limited, while still mitigating the impact of potentially noisy data.

Our proposed framework was tailored specifically for medical images, which
typically have a higher degree of novelty than other types of data. We
demonstrated that our framework not only preserved all the information in the
data but also resulted in more robust semi-supervised learning. Our
experiments showed that our approach outperformed other semi-supervised
learning frameworks on two sets of histopathology images, demonstrating the
effectiveness of our algorithms.

In conclusion, our multi-stage framework is a promising approach to address
the open-set semi-supervised learning problem in histopathology. By limiting
the exposure of open-set samples to the semi-supervised learning framework, we
were able to maintain the integrity of the training data while achieving
robust and effective learning. Our work has important implications for the
field of medical image analysis, where limited labeled data and high novelty
in unlabeled data are common challenges.

## Acknowledgment

The authors would like to thank Dr Stephanie McGregor, Assistant Professor at
the University of Wisconsin School of Medicine and Public Health for their
contributions in organizing the experiment data.

## References

  * [1] Mingyu Kim, Jihye Yun, Yongwon Cho, Keewon Shin, Ryoungwoo Jang, Hyun-jin Bae, and Namkug Kim,  “Deep learning in medical imaging,”  Neurospine, vol. 16, no. 4, pp. 657, 2019.
  * [2] Nikhil Cherian Kurian, Amit Sethi, Anil Reddy Konduru, Abhishek Mahajan, and Swapnil Ulhas Rane,  “A 2021 update on cancer image analytics with deep learning,”  Wiley Interdisciplinary Reviews: Data Mining and Knowledge Discovery, vol. 11, no. 4, pp. e1410, 2021.
  * [3] Amber L Simpson, Michela Antonelli, Spyridon Bakas, Michel Bilello, Keyvan Farahani, Bram Van Ginneken, Annette Kopp-Schneider, Bennett A Landman, Geert Litjens, Bjoern Menze, et al.,  “A large annotated medical image dataset for the development and evaluation of segmentation algorithms,”  arXiv preprint arXiv:1902.09063, 2019.
  * [4] Qing Yu, Daiki Ikami, Go Irie, and Kiyoharu Aizawa,  “Multi-task curriculum framework for open-set semi-supervised learning,”  in European Conference on Computer Vision. Springer, 2020, pp. 438–454.
  * [5] Jongjin Park, Sukmin Yun, Jongheon Jeong, and Jinwoo Shin,  “Opencos: Contrastive semi-supervised learning for handling open-set unlabeled data,”  in Computer Vision–ECCV 2022 Workshops: Tel Aviv, Israel, October 23–27, 2022, Proceedings, Part II. Springer, 2023, pp. 134–149.
  * [6] Nikhil Cherian Kurian, Pragati Shuddhodhan Meshram, Abhijeet Patil, Sunil Patel, and Amit Sethi,  “Sample specific generalized cross entropy for robust histology image classification,”  in 2021 IEEE 18th International Symposium on Biomedical Imaging (ISBI). IEEE, 2021, pp. 1934–1938.
  * [7] Huixiang Luo, Hao Cheng, Fanxu Meng, Yuting Gao, Ke Li, Mengdan Zhang, and Xing Sun,  “An empirical study and analysis on open-set semi-supervised learning,”  arXiv preprint arXiv:2101.08237, 2021.
  * [8] Jakob Nikolas Kather, Johannes Krisam, Pornpimol Charoentong, Tom Luedde, Esther Herpel, Cleo-Aron Weis, Timo Gaiser, Alexander Marx, Nektarios A Valous, Dyke Ferber, et al.,  “Predicting survival from colorectal cancer histology slides using deep learning: A retrospective multicenter study,”  PLoS medicine, vol. 16, no. 1, pp. e1002730, 2019.
  * [9] John N Weinstein, Eric A Collisson, Gordon B Mills, Kenna R Shaw, Brad A Ozenberger, Kyle Ellrott, Ilya Shmulevich, Chris Sander, and Joshua M Stuart,  “The cancer genome atlas pan-cancer analysis project,”  Nature genetics, vol. 45, no. 10, pp. 1113–1120, 2013.
  * [10] Han Zhang, Zizhao Zhang, Augustus Odena, and Honglak Lee,  “Consistency regularization for generative adversarial networks,”  arXiv preprint arXiv:1910.12027, 2019.
  * [11] David Berthelot, Nicholas Carlini, Ian Goodfellow, Nicolas Papernot, Avital Oliver, and Colin A Raffel,  “Mixmatch: A holistic approach to semi-supervised learning,”  Advances in neural information processing systems, vol. 32, 2019\.
  * [12] Hongyi Zhang, Moustapha Cisse, Yann N Dauphin, and David Lopez-Paz,  “mixup: Beyond empirical risk minimization,”  arXiv preprint arXiv:1710.09412, 2017.
  * [13] Kihyuk Sohn, David Berthelot, Nicholas Carlini, Zizhao Zhang, Han Zhang, Colin A Raffel, Ekin Dogus Cubuk, Alexey Kurakin, and Chun-Liang Li,  “Fixmatch: Simplifying semi-supervised learning with consistency and confidence,”  Advances in neural information processing systems, vol. 33, pp. 596–608, 2020.
  * [14] David Berthelot, Nicholas Carlini, Ekin D Cubuk, Alex Kurakin, Kihyuk Sohn, Han Zhang, and Colin Raffel,  “Remixmatch: Semi-supervised learning with distribution alignment and augmentation anchoring,”  arXiv preprint arXiv:1911.09785, 2019.
  * [15] Ozan Ciga, Tony Xu, and Anne Louise Martel,  “Self supervised contrastive learning for digital histopathology,”  Machine Learning with Applications, vol. 7, pp. 100198, 2022.
  * [16] Ting Chen, Simon Kornblith, Mohammad Norouzi, and Geoffrey Hinton,  “A simple framework for contrastive learning of visual representations,”  in Proceedings of the 37th International Conference on Machine Learning, Hal Daumé III and Aarti Singh, Eds. 13–18 Jul 2020, vol. 119 of Proceedings of Machine Learning Research, pp. 1597–1607, PMLR.
  * [17] Douglas A Reynolds,  “Gaussian mixture models.,”  Encyclopedia of biometrics, vol. 741, no. 659-663, 2009.
  * [18] Todd K Moon,  “The expectation-maximization algorithm,”  IEEE Signal processing magazine, vol. 13, no. 6, pp. 47–60, 1996\.
  * [19] Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun,  “Deep residual learning for image recognition,”  in Proceedings of the IEEE conference on computer vision and pattern recognition, 2016, pp. 770–778.
  * [20] Sergey Zagoruyko and Nikos Komodakis,  “Wide residual networks,”  arXiv preprint arXiv:1605.07146, 2016.
  * [21] Neeraj Kumar, Peter H Gann, Stephanie McGregor, and Amit Sethi,  “Quantification of subtype purity in luminal a breast cancer predicts clinical characteristics and survival,”  medRxiv, pp. 2023–02, 2023.


This manuscript has been authored by UT-Battelle, LLC under Contract No. DE-
AC05-00OR22725 with the U.S. Department of Energy. The United States
Government retains and the publisher, by accepting the article for
publication, acknowledges that the United States Government retains a non-
exclusive, paid-up, irrevocable, world-wide license to publish or reproduce
the published form of this manuscript, or allow others to do so, for United
States Government purposes. The Department of Energy will provide public
access to these results of federally sponsored research in accordance with the
DOE Public Access Plan(http://energy.gov/downloads/doe-public-access-plan).

# Microscopic view of heat capacity of matter: solid, liquid, and gas

Jaeyun Moon  To whom correspondence should be addressed; E-mail:
<EMAIL_ADDRESS>Materials Science and Technology Division,
Oak Ridge National Laboratory, Oak Ridge, Tennessee 37831, USA  Simon Thébaud
Materials Science and Technology Division,
Oak Ridge National Laboratory, Oak Ridge, Tennessee 37831, USA  Lucas Lindsay
Materials Science and Technology Division,
Oak Ridge National Laboratory, Oak Ridge, Tennessee 37831, USA  Takeshi Egami
Materials Science and Technology Division,
Oak Ridge National Laboratory, Oak Ridge, Tennessee 37831, USA
Department of Materials Science and Engineering,
University of Tennessee, Knoxville, Tennessee 37996, USA
Department of Physics and Astronomy,
University of Tennessee, Knoxville, Tennessee 37996, USA

(August 27, 2024)

###### Abstract

Understanding thermodynamics in liquids at the atomic level is challenging
because of strong atomic interactions and lack of symmetry. Recent prior
theoretical works have focused on describing heat capacity of liquids in terms
of phonon-like excitations but often rely on fitting parameters and ad hoc
assumptions. In this work, we perform microscopic analysis on instantaneous
normal modes and velocity autocorrelations on molecular dynamics simulations
of single element systems over wide ranges of temperature (up to 108 K) and
pressure (up to 1 TPa). Our results demonstrate that heat capacity of liquids
can be described by a combination of both solid-like and gas-like degrees of
freedom, leading to a unified framework to describe heat capacity of all three
phases of matter: solid, liquid, and gas.

In various applications of liquids from nuclear reactors
[williams_assessment_2006-1] to thermal energy storage devices
[li_sensible_2016], performance and efficiency of such device are directly
controlled by heat capacity. However, engineering and customization of heat
capacity of liquids are challenging due to a general lack of microscopic
understanding of the thermodynamics of liquids. Liquids have disordered
structures lacking spatial periodicity, yet their physical densities are
similar to solids. In gases, atoms are similarly disordered; however, liquids
are strongly interacting and dynamically correlated. These characteristics
make rigorous theoretical treatment of liquid thermodynamics a daunting
challenge. As such, the physics of liquids has been historically studied
starting from more established gas and solid perspectives.

Considering liquid and gas states as a unified concept dates back two
centuries to Cagniard de la Tour [de_la_tour_expose_1822], Faraday, Whewell,
Dumas [faraday_xvii_1823, faraday_selected_1971], Mendeleev
[mendelejeff_ueber_1861], Andrews [andrews_bakerian_1869,
tait_scientific_1889], and van der Waals [van_der_waals_over_1873]. The work
by van der Waals also motivated development of hard-sphere paradigm that is
widely used to study soft matter, granular materials, gases, and liquids
[chandler_van_1983, dyre_simple_2016]. On the other hand, many have also
suggested that liquids can be modeled by atoms as organized in crystals
[mie_zur_1903, mott_resistance_1934, lennard-jones_critical_1937,
granato_specific_2002]. Frenkel considered liquids and solids as a continuity
[frenkel_uber_1926, frenkel_continuity_1935, frenkel_liquid_1937,
frenkel_kinetic_1947] and unified these phases of matter as “condensed
bodies”. Frenkel proposed that when the characteristic time of interest is
less than $\tau=\eta/\mathcal{N}$ where $\eta$ is viscosity and $\mathcal{N}$
is rigidity modulus, materials behave like solid, and otherwise like liquid
[frenkel_liquid_1937]. The use of the term “condensed matter” has been since
popularized by Anderson and Heine to include liquid physics in the 1960s.
Frenkel’s seminal ideas have led to many recent works considering
thermodynamics of liquids from phonon quasi-particles as in solids
[trachenko_heat_2011, bolmatov_phonon_2012, andritsos_heat_2013,
trachenko_collective_2016, wang_direct_2017, yang_emergence_2017,
brazhkin_liquid-like_2018, tomiyoshi_heat_2019, khusnutdinoff_collective_2020,
Kryuchkov_universal_2020, zaccone_universal_2021, baggioli_explaining_2021].

Motivated by these historical accounts studying liquids from either gas or
solid perspectives, we consider here liquid thermodynamics and thermal
properties as a subset of a more general, unified framework that includes all
classical phases of matter simultaneously: solids, liquids, and gases. To
obtain microscopic understanding of thermodynamics and thermal properties of
matter, it is imperative to characterize the effective heat carriers and
understand their nature. We build conclusions of effective heat carriers in
liquids from systems over a very broad range of thermodynamic states rather
than from a limited temperature and pressure ranges typically reported in the
literature [bolmatov_phonon_2012, baggioli_explaining_2021,
zaccone_universal_2021, tomiyoshi_heat_2019, Kryuchkov_universal_2020].

Here we focus our discussion on heat carriers related to atomic motion. In
solids, they are phonons characterized by normal modes of the lattice
vibration. Phonon energies ($\hbar\omega$) and their spectral distribution,
$g(\omega)$, i.e., the density of states, are typically calculated by
diagonalizing the dynamical matrices built from atomic positions and force
constants or using velocity autocorrelation functions (VACF($\omega$)) from
molecular dynamics - both methods yield equivalent results under the harmonic
approximation. In this approximation, total energy and constant volume heat
capacity can be expressed in terms of phonons by $E=\int
d\omega\hbar\omega(n+\frac{1}{2})g(\omega)$ and $C_{V}=(dE/dT)_{V}=\int
d\omega
k_{B}(\frac{\hbar\omega}{2k_{B}T})^{2}\sinh^{-2}(\hbar\omega/2k_{B}T)g(\omega)$,
where $n$ is the Bose-Einstein distribution, and corrections can be further
made to include phonon anharmonicity [togo_first_2015]. In the classical
limit, these become $E=3Nk_{B}T$ and $C_{V}=3Nk_{B}$, also known as the
Dulong-Petit law of specific heat. Quadratic potential and kinetic energies
lead to equipartition of energy with $\frac{1}{2}k_{B}T$ per degree of
freedom. For gases where potential interactions are weak compared to thermal
energy, an atomic particle picture is used to describe thermodynamics and
thermal properties. Diagonalizing dynamical matrices of a gas state is rarely
used even for dense gases where potential interactions are relevant, but
VACF($\omega$) yields well-known Lorentzian lineshapes characterizing
independent uncorrelated atomic collisions. For monatomic dilute gases, in
which classical approximations are typically acceptable, total energy is
equivalent to the total kinetic energy,
$E=E_{KE}=\sum_{i}\frac{1}{2}m|\bm{v}_{i}|^{2}=\frac{3}{2}Nk_{B}T$, where the
sum is over all atoms and the corresponding specific heat is
$C_{V}=\frac{3}{2}Nk_{B}$.

These solid and gas limits ($C_{V}=3Nk_{B}$ vs. $\frac{3}{2}Nk_{B}$) suggest
that the effective heat carriers in liquids may potentially be described as an
intermediate state with a varying degree of importance of potential
interactions. In contrast, widely used heat capacity theories of liquids
typically rely on phonon quasi-particles alone [trachenko_heat_2011,
bolmatov_phonon_2012, andritsos_heat_2013, trachenko_collective_2016,
wang_direct_2017, yang_emergence_2017, brazhkin_liquid-like_2018,
tomiyoshi_heat_2019, khusnutdinoff_collective_2020, Kryuchkov_universal_2020,
zaccone_universal_2021, baggioli_explaining_2021] and describe the decrease in
heat capacity of liquids with increase in temperature solely due to a
disappearance of particular phonon populations [zaccone_universal_2021,
baggioli_explaining_2021]. Thus, the full translational atomic degrees of
freedom ($3N$) are not taken into account [zaccone_universal_2021,
baggioli_explaining_2021]. In addition, these theories are based on
assumptions of Debye densities of states [bolmatov_phonon_2012] or Gaussian
densities of states of normal modes at high frequencies with free fitting
parameters [zaccone_universal_2021, baggioli_explaining_2021].

In this work, we systematically characterize effective heat carriers of
monatomic systems (argon, silicon, and iron) under constant volume from solid
to gas through both instantaneous normal mode analysis and spectral velocity
autocorrelation functions over a wide range of temperatures from 1 to 108 K.
Based on these microscopic calculations, we define general instability
parameters (IP1 and IP2) to describe the ‘gasness’ of a system, and we draw
connection between seemingly different instantaneous normal mode distributions
and spectral velocity autocorrelation functions in liquid and gas states. We
further demonstrate agreement between constant volume specific heat predicted
from IPs and independent calculations using molecular dynamics (MD). These
results pave the way towards a unified approach to thermodynamics of matter
and thermal properties from solid to gas and provide new insights into a
microscopic view of liquid thermodynamics.

Interatomic interactions are described by Lennard-Jones (argon)
[jones_determination_1924, jones_determination_1924-1, rahman_molecular_1976],
Stillinger-Weber (silicon) [stillinger_computer_1985], and modified Johnson
potentials (iron) [srolovitz_structural_1981, levashov_equipartition_2008]. We
focus our scope on specific heat over a wide range of temperatures maintaining
a constant volume determined by the equilibrium density, $\rho_{0}$, at 1 K
and 1 bar for each system to minimize the effect of anharmonicity. Effects of
volumetric expansion and the role of anharmonicity may be addressed in future
works. Detailed simulation procedures are discussed in Supplementary
Materials.

We perform molecular dynamics (MD) simulations using Large-scale Atomic/
Molecular Massively Parallel Simulator (LAMMPS) [plimpton_fast_1995] to (1)
generate equilibrated atomic structures at a given temperature, (2) compute
spectral velocity autocorrelation functions, and (3) calculate specific heats.
We further use the equilibrated atomic structures generated from MD to perform
lattice dynamics calculations (GULP [gale_gulp:_1997] and in-house codes) to
obtain eigenfrequencies of the $\Gamma$ point dynamical matrices with the
entire domain considered as a unit cell. For better statistics, 10 structures
at each temperature in MD simulations were used for the lattice dynamics
calculations. Three independent MD simulations using different initial
velocities at each temperature were done for specific heat calculations.

The same procedures were additionally applied to structures at different
densities to characterize the sensitivity to density and anharmonicity:
$0.9\rho_{0}$ and $0.8\rho_{0}$ for argon, $1.1\rho_{0}$ for silicon, and
$0.9\rho_{0}$ for iron. For silicon, due to a large relative increase in
coordination numbers from 4 to more than 6 upon melting, reducing density to
$0.9\rho_{0}$ and $0.8\rho_{0}$ led to segregation of atoms and large empty
space in the simulation domains for some liquid temperatures. Thus, we instead
chose a more dense $\rho=1.1\rho_{0}$ for silicon.

To demonstrate that we are sampling different liquid and gas states, we
examine pair distribution functions (PDF), $g(r)$, of all $\rho_{0}$ systems.
PDF is defined as $g(r)=\frac{1}{4\pi
Nnr^{2}}\sum_{i,j}\langle\delta(r-|\bm{r}_{i}-\bm{r}_{j}|)\rangle$ where $N$
is the number of atoms, $n$ is the number density, $\bm{r}_{i}$ is the atomic
position of the $i$th atom, and the angled bracket denotes an ensemble
average. PDFs above melting temperatures are shown in Fig. 1. As we increase
the temperature, we observe progressive disappearance of well-defined peaks
and valleys in the $g(r)$ highlighting a transition from liquid to dense gas
states. The first distance point at which $g(r)$ becomes finite also decreases
with increase in temperature as the distance where potential energy $\sim
k_{B}T$ becomes smaller, as expected. Thus, the effective atomic diameter
becomes smaller with increase in temperature and the effective density
decreases. The compressibility factor, $Z=\frac{pV}{Nk_{B}T}$, which is a
measure of ideal gasness, is around 1.05 for the highest temperatures for all
systems, confirming that our systems are indeed in gas states despite the high
nominal density.

We now characterize their effective heat carriers. We first examine their
instantaneous normal modes by diagonalizing their dynamical matrices.
Resulting instantaneous normal mode densities of states, INM($\omega$), for
all $\rho_{0}$ systems beyond the melting temperatures are shown in Fig. 2. We
have done lattice dynamics calculations on small systems consisting of $\sim$
600 atoms and larger systems of $\sim$ 8000 atoms at select temperatures for
which we did not observe qualitative differences in the spectral
distributions. We expect slight size effects at low frequencies
$\omega\lesssim 1$ THz; however, low frequency modes constitute only a small
portion of the overall mode population and should not affect our results. In
addition, specific heat does not have strong size effects when normalized by
number of atoms so calculations on systems with $\sim$ 2000 atoms are
sufficient here. Negative frequencies denote modes with imaginary frequencies
arising from instabilities of the structure. Lightest blue shades represent
lowest temperatures considered above the melting temperature and the shades
become progressively red with increase in temperature. At the lowest
temperatures, INM($\omega$) for all systems are dominated by real modes as
measured by the areas under the curves. As temperature is increased, however,
imaginary mode populations become more prominent as observed previously
[keyes_instantaneous_1997].

We note that relaxed solids have only real modes and we observe that numbers
of real and imaginary modes become nearly equal at high temperatures in the
gas limit from our calculations. The square of the instantaneous normal mode
frequency reflects the curvature of the local potential energy landscape (PEL)
that atoms participating in the normal mode see at that instant. The concept
of the PEL [debenedetti_supercooled_2001] is usually applied to the whole
system, but here we expand the concept and define the PEL for the subsystem of
the normal mode. Such a limited view of the PEL is not new. For instance, when
a relaxation of glass is considered, only the atoms that are involved in the
relaxation phenomenon are taken into account in depicting the PEL
[goldstein_viscous_1969, debenedetti_supercooled_2001]. At low temperatures
the system is largely trapped in the valleys of the PEL, so the instantaneous
normal mode frequencies are mostly real. In contrast at very high temperatures
the system samples the positions with positive curvature (valleys) as well as
the positions with negative curvature (hills) equally. For symmetric sparse
random matrices, the eigenvalue distribution results in Wigner’s semi-circle
law with even number of positive and negative eigenvalues
[wigner_distribution_1958]. Thus, it is reasonable that the number of positive
and negative eigenvalues are equal for dynamical matrices of high temperature
gases. To describe this transition of instantaneous normal mode spectra from
solid to gas, we propose two phenomenological parameters called instability
parameters (IP1 and IP2) as a measure of how unstable the system is in the
configurational space. We define $\text{IP}_{1}=\frac{2N_{i}}{3N}$, where
$N_{i}$ is the total number of imaginary modes in our quantized systems, as a
measure of the ‘gasness’ of the system. IP1 is linearly related to the
fraction of the PEL with negative curvature that the system sees. We also
define $\text{IP}_{2}=\frac{N_{i}}{3N-N_{i}}$, which represents the ratio of
the portion of the PELs with negative curvature to those with positive
curvature. In both cases, IP1,2 $=0$ represents solid, IP1,2 $=1$ represents
gas, and liquids are described by IP1,2 values in-between. Temperature
dependent IP1,2 values for all $\rho_{0}$ systems are shown in Fig. S2 in
Supplementary Materials. IP1,2 increases with temperature for all systems as
expected and continues to increase slightly in the gas phase, possibly due to
finite size effects.

With these observations in INM($\omega$), we turn to velocity autocorrelation
spectra (calculated via MD simulations) of liquids and gases and see if we can
obtain a consistent picture of effective heat carriers existing in solid,
liquid, and gas phases. The spectral velocity autocorrelation function,
VACF($\omega$), is a useful tool to study atomic dynamics. VACF($\omega$)
describes phonon density of states for solids under harmonic approximations
and are equivalent to INM($\omega$) [dove_introduction_1993] as demonstrated
in Fig. S3 in Supplementary Materials for crystalline silicon at 1 K.
VACF($\omega$) can also describe non-phononic dynamics as VACF($\omega$ = 0)
represents the self-diffusion coefficient in the system (in our normalization,
$D=\frac{k_{B}T}{12mN}\text{VACF}(\omega=0)$ where $D$ and $m$ represent self-
diffusion coefficient and atomic mass, respectively). As such, there have been
prior attempts to decompose VACF($\omega$) into diffusion and phonon
contributions from both MD simulations and experiments utilizing Langevin-
Brownian motion and hard-sphere diffusion processes with varying degrees of
success [rahman_correlations_1964, verkerk_velocity_1989]. We take a different
approach and examine INM($\omega$) and VACF($\omega$) together to characterize
the heat carriers. At high temperatures, we observe Lorentzian lineshapes for
all gas phases demonstrated in Fig. 3 (C, I, F) as expected. Similar to
INM($\omega$), integration of VACF($\omega$) over frequency leads to the total
translational atomic degrees of freedom, $3N$. This property is where we make
connections between VACF($\omega$) and INM($\omega$). In the case of a gas
(IP1,2 $=1$ for our analysis), VACF($\omega$) is Lorentzian with a peak given
by $D$ and integral describing all $3N$ atomic degrees of freedom, no phonon
quasi-particles. For a solid (IP1,2 $=0$), VACF($\omega$) strictly describes
phonon density of states, with no free diffusing atoms. Motivated by this
contrast, we partition VACF($\omega$) into a gas-like portion,
$\text{VACF}_{\text{gas}}(\omega)$, with negligible attractive potential
interactions and a solid-like portion, $\text{VACF}_{\text{solid}}(\omega)$,
with strong potential interactions such that

$\frac{1}{3N}\int\text{VACF}_{\text{gas}}(\omega)d\omega=\text{IP}_{1,2}$ (1)
$\frac{1}{3N}\int\text{VACF}_{\text{solid}}(\omega)d\omega=1-\text{IP}_{1,2}$
(2)

VACFgas($\omega$) is assumed to have a Lorentzian lineshape with the height
and width determined by the diffusion coefficient and IP1,2, respectively.
VACFsolid($\omega$) is then determined by subtracting VACF($\omega$) by
VACFgas($\omega$). Calculated VACF($\omega$) and corresponding gas-like and
solid-like decompositions for all systems at select low, intermediate, and
high temperatures are shown in Fig. 3. If the ‘gasness’ in the system is
severely overpredicted, hence affecting the Lorenzian linewidth,
$\text{VACF}_{\text{solid}}(\omega)$ will be negative over a wide range of
frequencies, which is unphysical. The results here satisfy this test.

We first note that solid-like populations continuously decrease with
increasing temperature as expected from other prior phonon works
[bolmatov_phonon_2012, zaccone_universal_2021, baggioli_explaining_2021] on
liquids. However, rather than phononic degrees of freedom simply disappearing
[zaccone_universal_2021, baggioli_explaining_2021], they transition to gas-
like degrees of freedom and are eventually all absorbed into the gas-like
population at high temperatures as shown in Fig. 3 (C, F, I). This observation
in combination with our pair distribution function analysis at high
temperatures (see Fig. 1) leads to an important insight into INM($\omega$) of
liquid and gas states: real frequency modes do not have the same meaning as
conventional harmonic oscillators (phonons) in solids. Prior instantaneous
normal mode works have generally focused on the origin of imaginary modes
alone and simply considered real modes to derive from harmonic oscillators in
solids [stassen_instantaneous_1994, melzer_instantaneous_2012]. However,
average collision frequencies for high temperature liquids and gases,
conservatively estimated by inverse of the time required for atoms to travel
one interatomic distance ($(V/N)^{1/3}$) is well above 100 THz, higher than
most of the real normal mode frequencies in Fig. 2. There are nearly
$\frac{3}{2}N$ modes with real frequencies in INM($\omega$) for high
temperature gases and interpreting them as having conventional phonon harmonic
oscillators is highly questionable. Our work connecting INM($\omega$) and
VACF($\omega$) calls for a new interpretation of these real frequency modes
for liquid and gas states in terms of the local curvature in the PEL, and
using the instability parameter as a measure of ‘gasness’ to interpret the
correlated INM and VACF behaviors.

With VACF($\omega$) decomposed into VACFsolid($\omega$) and VACFgas($\omega$),
we now assess heat capacity contributions from the separate solid-like and
gas-like pictures. We assume that quadratic terms in the time-dependent
potential energy Taylor expansion is dominant for the constant volume systems
with $\rho_{0}$ for VACFsolid($\omega$). The total energy considered here
under classical approximation is given by

$E=\int d\omega\text{VACF}_{solid}(\omega)k_{B}T+\int
d\omega\text{VACF}_{gas}(\omega)\frac{1}{2}k_{B}T$ (3)

Eq. 3 can be further simplified in terms of IP1,2 by
$E=(1-\text{IP}_{1,2})(3Nk_{B}T)+\text{IP}_{1,2}\big{(}\frac{3}{2}Nk_{B}T)$.
The corresponding constant volume heat capacity for solids, liquids, and gases
under the classical and harmonic approximations is then simply

$C_{V,\text{IP}_{1,2}}=(1-\text{IP}_{1,2})(3Nk_{B})+\text{IP}_{1,2}\Big{(}\frac{3}{2}Nk_{B}\Big{)}-\frac{d\text{IP}_{1,2}}{dT}\Big{(}\frac{3}{2}Nk_{B}T\Big{)}$
(4)

Instability parameter derivatives at each $T$ were found from additional
instability parameters calculated at adjacent temperatures ($0.9T$ and
$1.1T$), i.e. numerical derivatives.

We can now compare our predictions from Eq. 4 with heat capacity independently
calculated by $C_{V,MD}=\frac{\langle E\rangle^{2}-\langle
E^{2}\rangle}{k_{B}T^{2}}$ from molecular dynamics (see Fig. 4). We see good
agreement between our predictions from Eq. 4 and MD specific heat values for
all systems. Clear convergence to the gas limit of $C_{V}=1.5Nk_{B}$ at IP
$=1$ is also observed as predicted (see Fig. S5 in Supplementary Materials).
$C_{V,\text{IP}_{1}}$ tends to underestimate $C_{V,MD}$ slightly, whereas
$C_{V,\text{IP}_{2}}$ shows a better overall agreement. However, in spite of
small differences, these observations give strong evidence supporting our
interpretations of effective heat carriers in INM($\omega$) and VACF($\omega$)
of solid, liquid, and gas phases from a unified perspective via simply defined
microscopic instability parameters.

We anticipate two possible sources of error when comparing the specific heat
values. The first one is the theory itself, i.e., IP1,2 does not represent
‘gasness’ correctly. The other source of error is neglecting anharmonicity in
potential interactions. To minimize the effect of volume expansion/contraction
on the potential interactions, we have focused our discussion primarily to
$\rho_{0}$ systems. Temperature dependent MD specific heat values are shown in
Fig. S4 in Supplementary Materials. We have also computed specific heat and
IP1,2 for densities slightly different from $\rho_{0}$ to examine the
sensitivity of volume effect near $\rho_{0}$ as shown in Fig. 4. We note that
our work does not assume Debye or Gaussian densities of states that are often
used in the literature [bolmatov_phonon_2012, zaccone_universal_2021,
baggioli_explaining_2021]. Rather, we examined the actual INM($\omega$) and
VACF($\omega$) from realistic potentials and make connections between the two
spectra of liquid and gas systems. Direct comparison with experiments is
challenging as most heat capacity measurements are done under constant
pressure conditions. Further, it is expected that anharmonicity will be
important near melting and glass transition temperatures. Pressure and volume
dependence and the role of anharmonicity in heat capacity merit further
investigations.

We expect that our work will also be useful in studying thermodynamics and
thermal properties of non-conventional materials including liquid crystals and
solid ionic conductors. There has been a lot of recent research interests for
thermoelectric power generators in certain solid ionic conductors dubbed
‘phonon-liquid, electron crystal’ where atoms at sub-lattice sites are fixed
while others are diffusing [liu_copper_2012, voneshen_hopping_2017,
niedziela_selective_2019]. This leads to desirable low heat capacity and non-
electronic thermal conductivity but the origin of these is not clear. It is
possible that our work helps identify the mechanism behind these complex
phenomena.

In summary, we have addressed thermodynamics of liquids from both solid and
gas perspectives. We propose to characterize effective heat carriers in
liquids via instability parameters describing ‘gasness’ in both instantaneous
normal mode and velocity autocorrelation spectra. In our approach we interpret
instantaneous normal modes reflecting the local curvature of the potential
energy landscape, rather than considering the instantaneous normal modes with
real frequencies equivalent to harmonic oscillators in solids as is often done
in the literature. We provide strong evidence in support of our proposal by
good agreement between predicted specific heat and specific heat values
calculated from molecular dynamics. Our work provides some insights into the
long-standing problem of thermodynamics of liquids and suggests pathways to a
unified framework in studying thermodynamics of solid, liquid, and gas phases.

This research was supported by the U.S. Department of Energy, Office of
Science, Basic Energy Sciences, Materials Sciences and Engineering Division.
This work used the Extreme Science and Engineering Discovery Environment
(XSEDE) Expanse under Allocation No. TG-MAT200012. This research used
resources of the National Energy Research Scientific Computing Center (NERSC),
a U.S. Department of Energy Office of Science User Facility located at
Lawrence Berkeley National Laboratory, operated under Contract No. DE-
AC02-05CH11231 using NERSC award BES-ERCAPERCAP0020503.

Figure 1: Pair distribution function, $g(r)$, of (A) argon from 103 to 108 K,
(B) silicon from $5\times 10^{3}$ to 107 K, and (C) iron from $5\times 10^{3}$
to 106 K. All systems were evaluated at $\rho_{0}$. Lightest blue and darkest
red curves represent lowest and highest temperatures, respectively. Well-
defined peaks and valleys in $g(r)$ progressively disappear as temperature
increases highlighting transition from liquid states at lower T to gas states
at higher T.

Figure 2: Instantaneous normal mode density of states (INM($\omega$)) of (A)
argon from 103 to 108 K, (B) silicon from $5\times 10^{3}$ to 107 K, and (C)
iron from $5\times 10^{3}$ to 106 K. All systems were evaluated at $\rho_{0}$.
Each INM($\omega$) is normalized such that the integral over the entire
frequency range is unity. Lightest blue and darkest red represent lowest and
highest temperatures, respectively. Negative frequency in the figure means
imaginary frequency. At low temperatures, modes with real frequency modes
dominate the spectra. As temperature increases, however, we observe that
imaginary mode populations become more significant.

Figure 3: Decomposition of spectral velocity autocorrelation function,
VACF($\omega$) at select temperatures of (A-C) argon from 103 to 108 K, (D-F)
silicon from $5\times 10^{3}$ to 107 K, and (G-I) iron from $5\times 10^{3}$
to 106 K. All systems were evaluated at $\rho_{0}$. Solid and dashed curves
represent results based on IP1 and IP2, respectively. Black curves represent
gas-like degrees of freedom spectra, VACFreal($\omega$), from instability
parameter and orange curves represent solid-like degrees of freedom spectra,
VACFphonon($\omega$), obtained from subtraction of VACF($\omega$) by the black
curves.

Figure 4: Predicted constant volume heat capacity of various monatomic systems
with different densities from instability parameters (Eq. 4), compared to
independently calculated heat capacity from molecular dynamics. Dashed black
lines represent one to one correspondence. For $C_{V,\text{IP}_{2}}$,
deviations are within 5%.


networks themselves are not critical due to the heterogeneity of the neurons.
The red curves are the predictions of the hierarchy of nonlinearities (Eqs.
(15), truncated at fourth order for subcritical and critical cases and first
order for the supercritical case, using the results of loh2011accurate  to
compute the eigenvalue distributions. Blue data points are simulated data,
using networks of $N=10^{3}$ neurons. Networks are created by generating a
symmetric Gaussian random matrix, diagonalizing it, and drawing new
eigenvalues from a Beta distribution SuppInfo .

## References

  * [1] R. E. Kass, S.-I. Amari, K. Arai, E. N. Brown, C. O. Diekman, M. Diesmann, B. Doiron, U. T. Eden, A. L. Fairhall, G. M. Fiddyment, T. Fukai, So. Grün, M. T. Harrison, M. Helias, H. Nakahara, J.-N. Teramae, P. J. Thomas, M. Reimers, J. Rodu, H. G. Rotstein, E. Shea-Brown, H. Shimazaki, S. Shinomoto, B. M. Yu, and M. A. Kramer.  Computational neuroscience: Mathematical and statistical perspectives.  Annual Review of Statistics and Its Application, 5(1):183–214, 2018\.
  * [2] A. A. Faisal, L.P. J. Selen, and D. M. Wolpert.  Noise in the nervous system.  Nat. Rev. Neurosci., 9(4):292–303, 2008.
  * [3] Mikhail I Rabinovich, Pablo Varona, Allen I Selverston, and Henry DI Abarbanel.  Dynamical principles in neuroscience.  Reviews of modern physics, 78(4):1213, 2006.
  * [4] Kechen Zhang.  Representation of spatial orientation by the intrinsic dynamics of the head-direction cell ensemble: a theory.  Journal of Neuroscience, 16(6):2112–2126, 1996.
  * [5] Carlo R Laing and Carson C Chow.  Stationary bumps in networks of spiking neurons.  Neural computation, 13(7):1473–1494, 2001.
  * [6] Klaus Wimmer, Duane Q Nykamp, Christos Constantinidis, and Albert Compte.  Bump attractor dynamics in prefrontal cortex explains behavioral precision in spatial working memory.  Nature neuroscience, 17(3):431–439, 2014.
  * [7] Sung Soo Kim, Hervé Rouault, Shaul Druckmann, and Vivek Jayaraman.  Ring attractor dynamics in the drosophila central brain.  Science, 356(6340):849–853, 2017.
  * [8] G Bard Ermentrout and Jack D Cowan.  A mathematical theory of visual hallucination patterns.  Biological cybernetics, 34(3):137–150, 1979.
  * [9] Paul C Bressloff.  Metastable states and quasicycles in a stochastic wilson-cowan model of neuronal population dynamics.  Physical Review E, 82(5):051903, 2010.
  * [10] Thomas Charles Butler, Marc Benayoun, Edward Wallace, Wim van Drongelen, Nigel Goldenfeld, and Jack Cowan.  Evolutionary constraints on visual cortex architecture from the dynamics of hallucinations.  Proceedings of the National Academy of Sciences, 109(2):606–609, 2012.
  * [11] H. Sompolinsky, A. Crisanti, and H. J. Sommers.  Chaos in random neural networks.  Phys. Rev. Lett., 61:259–262, Jul 1988.
  * [12] David Dahmen, Sonja Grün, Markus Diesmann, and Moritz Helias.  Second type of criticality in the brain uncovers rich multiple-neuron dynamics.  Proceedings of the National Academy of Sciences, 116(26):13051–13060, 2019.
  * [13] John M Beggs and Dietmar Plenz.  Neuronal avalanches in neocortical circuits.  Journal of neuroscience, 23(35):11167–11177, 2003.
  * [14] Michael A Buice and Jack D Cowan.  Field-theoretic approach to fluctuation effects in neural networks.  Physical Review E, 75(5):051919, 2007.
  * [15] N. Friedman, S. Ito, B. A. W. Brinkman, M. Shimono, R. E. L. DeVille, K. A. Dahmen, J. M. Beggs, and T. C. Butler.  Universal critical dynamics in high resolution neuronal avalanche data.  Phys. Rev. Lett., 108:208102, May 2012.
  * [16] Jonathan Kadmon and Haim Sompolinsky.  Transition to chaos in random neuronal networks.  Physical Review X, 5(4):041030, 2015.
  * [17] Woodrow L Shew, Hongdian Yang, Shan Yu, Rajarshi Roy, and Dietmar Plenz.  Information capacity and transmission are maximized in balanced cortical networks with neuronal avalanches.  Journal of neuroscience, 31(1):55–63, 2011.
  * [18] J. Beggs and N. Timme.  Being critical of criticality in the brain.  Frontiers in Physiology, 3:163, 2012.
  * [19] W. L. Shew and D. Plenz.  The functional benefits of criticality in the cortex.  The Neuroscientist, 19(1):88–100, 2013.  PMID: 22627091.
  * [20] Jonathan Touboul and Alain Destexhe.  Can power-law scaling and neuronal avalanches arise from stochastic dynamics?  PloS one, 5(2):e8982, 2010.
  * [21] Jonathan Touboul and Alain Destexhe.  Power-law statistics and universal scaling in the absence of criticality.  Physical Review E, 95(1):012413, 2017.
  * [22] Marcel Nonnenmacher, Christian Behrens, Philipp Berens, Matthias Bethge, and Jakob H Macke.  Signatures of criticality arise from random subsampling in simple population models.  PLoS computational biology, 13(10):e1005718, 2017.
  * [23] Anna Levina and Viola Priesemann.  Subsampling scaling.  Nature communications, 8(1):1–9, 2017.
  * [24] G. K. Ocker, K. Josić, E. Shea-Brown, and M. A. Buice.  Linking structure and activity in nonlinear spiking networks.  PLOS Computational Biology, 13(6):1–47, 06 2017.
  * [25] Jonas Stapmanns, Tobias Kühn, David Dahmen, Thomas Luu, Carsten Honerkamp, and Moritz Helias.  Self-consistent formulations for stochastic nonlinear neuronal dynamics.  Physical Review E, 101(4):042124, 2020.
  * [26] M. Kordovan and S. Rotter.  Spike train cumulants for linear-nonlinear poisson cascade models, 2020\.
  * [27] Gabriel Koch Ocker.  Dynamics of stochastic integrate-and-fire networks.  Phys. Rev. X, 12:041007, Oct 2022.
  * [28] L. Canet, B. Delamotte, O. Deloubrière, and N. Wschebor.  Nonperturbative renormalization-group study of reaction-diffusion processes.  Phys. Rev. Lett., 92:195703, May 2004.
  * [29] Léonie Canet, Hugues Chaté, Bertrand Delamotte, Ivan Dornic, and Miguel A Munoz.  Nonperturbative fixed point in a nonequilibrium phase transition.  Physical review letters, 95(10):100601, 2005.
  * [30] T. Machado and N. Dupuis.  From local to critical fluctuations in lattice models: A nonperturbative renormalization-group approach.  Phys. Rev. E, 82:041128, Oct 2010.
  * [31] L. Canet, H. Chaté, and B. Delamotte.  General framework of the non-perturbative renormalization group for non-equilibrium steady states.  Journal of Physics A: Mathematical and Theoretical, 44(49):495001, Nov 2011.
  * [32] A. Rançon and N. Dupuis.  Nonperturbative renormalization group approach to strongly correlated lattice bosons.  Phys. Rev. B, 84:174513, Nov 2011.
  * [33] A. A. Winkler and E. Frey.  Long-range and many-body effects in coagulation processes.  Phys. Rev. E, 87:022136, Feb 2013.
  * [34] Ingo Homrighausen, Anton A. Winkler, and Erwin Frey.  Fluctuation effects in the pair-annihilation process with lévy dynamics.  Phys. Rev. E, 88:012111, Jul 2013.
  * [35] T. Kloss, L. Canet, B. Delamotte, and N. Wschebor.  Kardar-parisi-zhang equation with spatially correlated noise: A unified picture from nonperturbative renormalization group.  Phys. Rev. E, 89:022108, Feb 2014.
  * [36] Ivan Balog and Gilles Tarjus.  Activated dynamic scaling in the random-field ising model: A nonperturbative functional renormalization group approach.  Physical Review B, 91(21):214201, 2015.
  * [37] Léonie Canet, Bertrand Delamotte, and Nicolás Wschebor.  Fully developed isotropic turbulence: Symmetries and exact identities.  Physical Review E, 91(5):053004, 2015.
  * [38] P. Jakubczyk and A. Eberlein.  Thermodynamics of the two-dimensional $\mathit{XY}$ model from functional renormalization.  Phys. Rev. E, 93:062145, Jun 2016.
  * [39] C. Duclut and B. Delamotte.  Nonuniversality in the erosion of tilted landscapes.  Phys. Rev. E, 96:012149, Jul 2017.
  * [40] Leenoy Meshulam, Jeffrey L Gauthier, Carlos D Brody, David W Tank, and William Bialek.  Coarse graining, fixed points, and scaling in a large population of neurons.  Physical review letters, 123(17):178103, 2019.
  * [41] Serena Bradde and William Bialek.  PCA meets RG.  Journal of statistical physics, 167(3):462–475, 2017.
  * [42] Lorenzo Tiberi, Jonas Stapmanns, Tobias Kühn, Thomas Luu, David Dahmen, and Moritz Helias.  Gell-mann–low criticality in neural networks.  Physical Review Letters, 128(16):168301, 2022.
  * [43] Daniel Martí, Nicolas Brunel, and Srdjan Ostojic.  Correlations between synapses in pairs of neurons slow down dynamics in randomly connected neural networks.  Physical Review E, 97(6):062314, 2018.
  * [44] Braden AW Brinkman, Fred Rieke, Eric Shea-Brown, and Michael A Buice.  Predicting how and when hidden neurons skew measured synaptic interactions.  PLoS computational biology, 14(10):e1006490, 2018.
  * [45] Bertrand Delamotte.  An introduction to the nonperturbative renormalization group.  In Renormalization group and effective field theory approaches to many-body systems, pages 49–132. Springer, 2012.
  * [46] Léonie Canet.  Reaction–diffusion processes and non-perturbative renormalization group.  Journal of Physics A: Mathematical and General, 39(25):7901, 2006\.
  * [47] L. Canet and H. Chaté.  A non-perturbative approach to critical dynamics.  Journal of Physics A: Mathematical and Theoretical, 40(9):1937–1949, feb 2007.
  * [48] Nicolas Dupuis, L Canet, Astrid Eichhorn, W Metzner, Jan M Pawlowski, M Tissier, and N Wschebor.  The nonperturbative functional renormalization group and its applications.  Physics Reports, 910:1–114, 2021.
  * [49] Braden A. W. Brinkman.  Supplemental material.
  * [50] C. Wetterich.  Exact evolution equation for the effective potential.  Physics Letters B, 301(1):90 – 94, 1993.
  * [51] Jean-Michel Caillol.  The non-perturbative renormalization group in the ordered phase.  Nuclear Physics B, 855(3):854–884, 2012.
  * [52] James Cooper Robinson and C Pierre.  Infinite-dimensional dynamical systems: An introduction to dissipative parabolic pdes and the theory of global attractors. cambridge texts in applied mathematics.  Appl. Mech. Rev., 56(4):B54–B55, 2003.
  * [53] Yen Lee Loh.  Accurate calculation of green functions on the d-dimensional hypercubic lattice.  Journal of Physics A: Mathematical and Theoretical, 44(27):275201, 2011.
  * [54] Juergen Berges, Nikolaos Tetradis, and Christof Wetterich.  Non-perturbative renormalization flow in quantum field theory and statistical physics.  Physics Reports, 363(4-6):223–386, 2002.
  * [55] Nigel Goldenfeld.  Lectures on Phase Transitions and the Renormalization Group.  Westview Press, 1992.
  * [56] Solving Burger’s equation with NDSolve at large time.  StackExchange; accessed January 18 2023.
  * [57] Haim Sompolinsky and Annette Zippelius.  Relaxational dynamics of the edwards-anderson model and the mean-field theory of spin-glasses.  Physical Review B, 25(11):6860, 1982.
  * [58] Jannis Schuecker, Sven Goedeke, and Moritz Helias.  Optimal sequence memory in driven random networks.  Physical Review X, 8(4):041029, 2018.
  * [59] Aslı Tuncer and Ayşe Erzan.  Spectral renormalization group for the gaussian model and $\psi$ 4 theory on nonspatial networks.  Physical Review E, 92(2):022106, 2015.
  * [60] Hans-Karl Janssen, Frédéric van Wijland, Olivier Deloubrière, and Uwe C. Täuber.  Pair contact process with diffusion: Failure of master equation field theory.  Phys. Rev. E, 70:056114, Nov 2004.
  * [61] Hans-Karl Janssen and Uwe C Täuber.  The field theory approach to percolation processes.  Annals of Physics, 315(1):147–192, 2005.
  * [62] Malo Tarpin, Federico Benitez, Léonie Canet, and Nicolás Wschebor.  Nonperturbative renormalization group for the diffusive epidemic process.  Physical Review E, 96(2):022137, 2017.
  * [63] Malte Henkel, Haye Hinrichsen, and Sven Lübeck.  Universality classes different from directed percolation.  Non-Equilibrium Phase Transitions: Volume I: Absorbing Phase Transitions, pages 197–259, 2008.
  * [64] Fan Zhong.  Imaginary fixed points can be physical.  Physical Review E, 86(2):022104, 2012.
  * [65] Fan Zhong.  Renormalization-group theory of first-order phase transition dynamics in field-driven scalar model.  Frontiers of Physics, 12(5):1–31, 2017.
  * [66] Xin An, David Mesterházy, and Mikhail A Stephanov.  Functional renormalization group approach to the yang-lee edge singularity.  Journal of High Energy Physics, 2016(7):1–19, 2016.
  * [67] Tim R Morris.  Properties of derivative expansion approximations to the renormalization group.  International Journal of Modern Physics B, 12(12n13):1343–1354, 1998.
  * [68] Léonie Canet, Bertrand Delamotte, Dominique Mouhanna, and Julien Vidal.  Nonperturbative renormalization group approach to the ising model: a derivative expansion at order $\partial^{4}$.  Physical Review B, 68(6):064421, 2003.
  * [69] Tim R Morris.  Derivative expansion of the exact renormalization group.  Physics Letters B, 329(2-3):241–248, 1994.
  * [70] J-P Blaizot, R Mendez Galain, and N Wschebor.  Non-perturbative renormalization group calculation of the transition temperature of the weakly interacting bose gas.  EPL (Europhysics Letters), 72(5):705, 2005.
  * [71] Uwe C Täuber.  Field-theory approaches to nonequilibrium dynamics.  In Ageing and the Glass Transition, pages 295–348. Springer, 2007\.
  * [72] Uwe C Täuber.  Renormalization group: applications in statistical physics.  Nuclear Physics B-Proceedings Supplements, 228:7–34, 2012.
  * [73] Uwe C Täuber.  Critical dynamics: a field theory approach to equilibrium and non-equilibrium scaling behavior.  Cambridge University Press, 2014.
  * [74] Wulfram Gerstner, Werner M Kistler, Richard Naud, and Liam Paninski.  Neuronal dynamics: From single neurons to networks and models of cognition.  Cambridge University Press, 2014.
  * [75] John J Hopfield.  Neural networks and physical systems with emergent collective computational abilities.  Proceedings of the national academy of sciences, 79(8):2554–2558, 1982.
  * [76] Sen Song, Per Jesper Sjöström, Markus Reigl, Sacha Nelson, and Dmitri B Chklovskii.  Highly nonrandom features of synaptic connectivity in local cortical circuits.  PLoS biology, 3(3):e68, 2005.
  * [77] Kevin Berlemont and Gianluigi Mongillo.  Glassy phase in dynamically-balanced neural networks.  bioRxiv, 2022.
  * [78] Jef Hooyberghs, Ferenc Iglói, and Carlo Vanderzande.  Absorbing state phase transitions with quenched disorder.  Physical Review E, 69(6):066140, 2004.
  * [79] Gilles Tarjus and Matthieu Tissier.  Nonperturbative functional renormalization group for random-field models: The way out of dimensional reduction.  Physical review letters, 93(26):267008, 2004.
  * [80] Gilles Tarjus and Matthieu Tissier.  Nonperturbative functional renormalization group for random field models and related disordered systems. i. effective average action formalism.  Physical Review B, 78(2):024203, 2008.
  * [81] Matthieu Tissier and Gilles Tarjus.  Nonperturbative functional renormalization group for random field models and related disordered systems. ii. results for the random field o (n) model.  Physical Review B, 78(2):024204, 2008.
  * [82] Gašper Tkačik, Thierry Mora, Olivier Marre, Dario Amodei, Stephanie E Palmer, Michael J Berry, and William Bialek.  Thermodynamics and signatures of criticality in a network of neurons.  Proceedings of the National Academy of Sciences, 112(37):11508–11513, 2015.
  * [83] Rashid V Williams-García, Mark Moore, John M Beggs, and Gerardo Ortiz.  Quasicritical brain dynamics on a nonequilibrium widom line.  Physical Review E, 90(6):062714, 2014.
  * [84] Adam G Kline and Stephanie E Palmer.  Gaussian information bottleneck and the non-perturbative renormalization group.  New Journal of Physics, 24(3):033007, 2022.
  * [85] A. Skrondal and S. Rabe-Hesketh.  Some applications of generalized linear latent and mixed models in epidemiology: Repeated measures, measurement error and multilevel modeling.  Norsk Epidemiologi, 13(2), Oct. 2003.
  * [86] D. I. Warton and L. C. Shepherd.  Poisson point process models solve the “pseudo-absence problem” for presence-only data in ecology.  Ann. Appl. Stat., 4(3):1383–1402, 09 2010.
  * [87] A. Bray and F. P. Schoenberg.  Assessment of point process models for earthquake forecasting.  Statist. Sci., 28(4):510–520, 11 2013.
  * [88] Braden AW Brinkman, Michael LeBlanc, Yehuda Ben-Zion, Jonathan T Uhl, and Karin A Dahmen.  Probing failure susceptibilities of earthquake faults using small-quake tidal correlations.  Nature communications, 6(1):1–7, 2015.
  * [89] Braden AW Brinkman, Michael P LeBlanc, Jonathan T Uhl, Yehuda Ben-Zion, and Karin A Dahmen.  Probabilistic model of waiting times between large failures in sheared media.  Physical Review E, 93(1):013003, 2016.