# VQSynergy: Robust Drug Synergy Prediction with Vector Quantization Mechanism

Jiawei Wu*
School of Data Science and Engineering
East China Normal University
Shanghai, China
&Mingyuan Yan
School of Data Science and Engineering
East China Normal University
Shanghai, China
&Dianbo Liu
School of Medicine and College of Design&Engineering
National University of Singapore
Singapore
* equal contribution

###### Abstract

The advancement in optimizing cancer therapies has been significantly boosted
by the ability to accurately predict drug synergies. While traditional
approaches like clinical trials offer reliability, they are hampered by
considerable time and financial requirements. The advent of high-throughput
screening and computational breakthroughs has marked a shift towards more
efficient methods for investigating drug interactions. In this study, we
introduce VQSynergy, an innovative machine learning framework that leverages
the Vector Quantization (VQ) mechanism, combined with gated residuals and a
custom attention mechanism, to improve the accuracygeneralizability of drug
synergy predictions. Our results show that VQSynergy outperforms existing
models in robustness, especially in conditions with Gaussian noise,
underscoring its superior efficacy and value in the intricate and often
unpredictable field of drug synergy research. This study highlights
VQSynergy’s potential to transform the field through its advanced predictive
capabilities, thus contributing to the refinement of cancer treatment
strategies.

_K_ eywords machine learning $\cdot$ drug synergy $\cdot$ vector quantization
$\cdot$ generalization $\cdot$ robustness

##  1 Introduction

Figure 1: Overview of VQSynergy

The treatment of cancer faces challenges due to its complex nature. While the
traditional ‘one gene, one drug’ approach proposed by Ehrlich [1, 2] remains
common for many cancer types, it targets actively proliferating cells
indiscriminately, causing damage to healthy tissue alongside cancerous cells.
This oversimplified view of cancer biology overlooks the intricate interplay
of factors involved. Further, monotherapy treatment is prone to drug
resistance as continuous exposure to a single compound prompts cancer cells to
activate alternative survival pathways [3, 4, 5]. In response to these
limitations, there has been a shift towards combination therapy, where drug
synergy offers the potential for increased efficacy at lower doses, thereby
reducing toxicity and the likelihood of drug resistance [6, 7]. This paradigm
shift emphasizes the need to reconsider our approach to cancer treatment,
moving beyond the simplistic ‘one tumor, one gene, one drug’ strategy towards
a more holistic understanding of cancer biology [8].

The landscape of drug synergy prediction has been revolutionized by advances
in high-throughput screening methods and computational modeling [9]. Despite
these advancements, the scarcity of synergistic drug pairs and the vast search
space pose challenges, necessitating the development of computational models
to prioritize viable combinations [9]. These technologies have led to the
construction of extensive datasets like O’Neil [10], NCI-ALMANAC [11], and
others, capturing over a million experimentally measured synergy scores. These
resources also introduce the complexity of handling vast amounts of data,
underscoring the need for efficient computational models.

In addressing this need, recent advancements in machine learning and deep
learning methods have shown promise in predicting drug synergy. DeepSynergy
[12] uses deep learning approach for predicting drug synergy, utilizing
chemical and genomics data to outperform traditional machine learning methods
in identifying novel synergistic drug combinations. The DTF [13] model
combines tensor factorization with a deep neural network (DNN) to predict drug
synergy. Furthermore, HypergraphSynergy [14] predicts anti-cancer drug synergy
by formulating synergistic drug combinations over cancer cell lines as a
hypergraph, which often grapples with the intricate multi-way relations
existing between drugs and cell lines. This study underscores the necessity to
employ models that can capture and interpret these multi-dimensional
relationships, leading to the development of hypergraph neural networks
(HGNNs) [15], which have shown promise in various fields, including biomedical
applications.

However, the message passing in HGNNs can hinder the model’s generalization,
especially with new node features absent in the training data. After message
passing, these new features might inadvertently degrade predictive
performance. To address this, a more meticulous approach to message passing is
required. Implementing Vector Quantization (VQ) [16] in the message passing
process presents a promising solution. By conceptualizing message passing as a
form of communication between different entities [17], the application of a
simpler ‘lingua franca’ can facilitate more efficient and effective
interactions within the hypergraph framework.

Previous studies in drug synergy prediction have been significantly hindered
by a heavy reliance on the quality of data pertaining to drug-associated
biological entities. This dependency has often limited their performance,
particularly in scenarios involving novel drugs or diseases where prior
information is sparse. Recognizing this limitation, we developed VQSynergy, a
robust framework for drug synergy prediction, employing a Vector Quantization
(VQ) mechanism and gated residuals, coupled with a specially designed
attention mechanism, VQSynergy not only achieves high predictive accuracy, but
also enhances the model’s generalization. This innovative approach ensures
that our model is not constrained by the availability of prior knowledge,
making it particularly effective for the early screening stages of drug
development.

##  2 Methods

###  2.1 Overview of VQSynergy

Our methodology reconceptualizes drug synergy prediction within an advanced
hypergraph framework, treating it as a hyperedge prediction problem. This
innovative approach unfolds through three interconnected phases 1:
initialization, refinement, and consolidation. Each phase plays a pivotal role
in processing and integrating varied data types, thereby synthesizing
comprehensive information crucial for accurate drug synergy predictions.
Besides, our model introduces novel methodologies and integrates state-of-the-
art techniques to substantially enhance the prediction accuracy and
reliability in each phase of the analysis.

###  2.2 Initialization: Feature Acquisition

The initialization stage is essential in our approach to predicting synergy
involving collecting initial features for drugs and cell lines, and aligning
them to a uniform dimensionality, which is critical for later analysis in our
hypergraph model.

Molecular Structure Data for Drugs: For each drug, we construct molecular
graphs $\mathbf{G}=(\mathbf{X}_{\text{atoms}},\mathbf{A})$ using the
Simplified Molecular Input Line Entry System (SMILES) [18]. In these graphs,
the attribute matrix $\mathbf{X}_{\text{atoms}}$, with dimensions
$\mathbb{R}^{n\times F_{d}}$, captures atomic attributes, and the adjacency
matrix $\mathbf{A}$, sized $\mathbb{R}^{n\times n}$, represents atomic
connections. To fit our model’s requirements for one-dimensional data
representations, we convert these multidimensional matrices into a singular
feature vector $\mathbf{x_{d}}$ for each drug, encapsulating its molecular
structure. This conversion relies primarily on pooling techniques. Before
pooling, we apply TransformerConv [19] to more effectively extract structural
details from the molecular graphs.

Genomic Data for Cell Lines: For cell lines, we refine gene expression
profiles into feature vectors $\mathbf{x_{c}}$ through a process of
normalization and transformation. This step is pivotal in capturing the
distinct genomic characteristics of each cell line. To ensure compatibility
with our model’s uniform dimensional requirements, we utilize a Multi-Layer
Perceptron (MLP). The MLP functions to project these genomic features into a
dimensionally consistent space. This projection is essential for integrating
the genomic features of cell lines with the molecular structure data of drugs
in the hypergraph model, ensuring a cohesive and comprehensive dataset for
analysis.

###  2.3 Refinement: Enhancing Feature Representation

In the refinement phase of our model, we transform raw data into enhanced
feature representations by constructing a drug synergy hypergraph. This
process involves utilizing Hypergraph Neural Network Convolution, integrating
gated residual connections, and employing Vector Quantization (VQ) techniques.
The combination of these methods, particularly the inclusion of VQ,
constitutes a meticulous message passing approach. This approach is crucial
for accurately identifying complex interactions between drugs and cell lines,
which are fundamental to effective drug synergy prediction. The incorporation
of VQ allows for more precise feature discretization, further refining the
model’s ability to understand and predict intricate interrelationships in drug
synergy applications. This phase, with its focus on detail and precision,
significantly enhances the model’s analytical capabilities, ensuring robust
and reliable predictions in the complex field of drug synergy.

####  2.3.1 Representing Multi-way Relations

In this phase, we create a hypergraph to represent synergistic drug-cell line
combinations. This hypergraph is essential for understanding the complex
interactions of drugs and their effects on various cell lines, aiding in
assessing the potential efficacy of drug combinations.

Nodes $v$ in the hypergraph, initialized as detailed in Section 2.2, are
described by the feature matrix $\mathbf{X}=\begin{bmatrix}\mathbf{X}_{D}\\\
\mathbf{X}_{C}\end{bmatrix}$. Here, $\mathbf{X}_{D}$ denotes molecular
structure features of drugs, and $\mathbf{X}_{C}$ indicates gene expression
features of cell lines.

The hypergraph structure is captured through an incidence matrix
$\mathbf{H}\in\mathbb{R}^{|V|\times|E|}$, focusing on the Drug-Cell Line
Synergy Relationship. This relationship is central to the hypergraph, as it
reveals the potential synergistic effects of specific drugs on different cell
lines, which is crucial for predicting the effectiveness of drug combinations.

####  2.3.2 Meticulous Message Passing in Graph

In this section, we present a novel method for refining features in Graph
Neural Networks (GNNs), emphasizing meticulous message passing. Our approach
integrates hypergraph convolution, quantization strategies, and gated residual
connections, tailored for modeling object interactions and predicting their
future states, while embedding relational information.

A major challenge in GNNs is “over-smoothing” [20], a phenomenon where
repeated convolutions cause node features to become overly similar and
indistinct. To combat this, we incorporate a residual connection to preserve
feature uniqueness. The node feature update at $l^{\text{th}}$ layer is
formulated as:

$\mathbf{X}^{(l+1)}=\mathbf{X}^{(l)}+\Delta\mathbf{X}^{(l)},$ (1)

where $\Delta\mathbf{X}^{(l)}$ represents the feature update. This update is
derived from message updates and gating functions, subsequently discretized
into a neural communication format [17], drawing inspiration from VQ-VAE [16].
For each layer, we utilize a discrete space vector $e\in\mathbb{R}^{L\times
m}$, denoted by $e^{(l)}$, with $L$ as the size of the discrete latent space
and $m$ representing the dimension of each latent vector $e_{j}$. We segment
each target vector into $G$ parts, or discretization heads, quantize each head
individually using the layer-specific codebook, and concatenate the outcomes:

$\Delta\mathbf{X}^{(l)}=q^{(l)}(\mathbf{m}^{(l)}\odot\mathbf{i}^{(l)},L,G),$
(2)

To manage multi-way relationships characteristic of hypergraphs, where one
hyperedge can connect multiple nodes, our model employs Hypergraph Neural
Networks (HGNN) [15, 21]. The message passing in HGNN is described as:

$\displaystyle\mathbf{m}^{(l)}=\text{ReLU}(\mathbf{D}^{-1}\mathbf{H}\mathbf{B}^{-1}\mathbf{H}^{\top}\mathbf{X}^{(l)}\Theta_{m}^{(l)}),$
(3) $\displaystyle\mathbf{D}=\sum_{e}Y_{ve},\quad\mathbf{B}=\sum_{v}Y_{ve}.$

where $\mathbf{D}$ and $\mathbf{B}$ denote diagonal matrices representing node
and hyperedge degrees, respectively.

A key element of our model is its detailed and carefully designed message
passing mechanism. This mechanism notably incorporates a gating function to
modulate the influence of adjacent nodes. We strategically initialize the
biases $b$ with negative values. This initial setting enables the gating
function to start functioning close to an identity operation, which is crucial
for the efficient integration of features.

For the $l^{\text{th}}$ layer, the gating function is defined as follows:

$\mathbf{i}^{(l)}=\sigma(\mathbf{X}^{(l)}\Theta_{i}^{(l)}),$ (4)

where $\sigma$ is the sigmoid function, crucial for preserving the uniqueness
of node features and , thereby contributing significantly to the model’s
learning dynamics.

The overall graph refinement update rule in our model is formulated as:

$\mathbf{X}^{(l+1)}=\mathbf{X}^{(l)}+q^{(l)}(\text{ReLU}(\mathbf{D}^{-1}\mathbf{H}\mathbf{B}^{-1}\mathbf{H}^{\top}\mathbf{X}^{(l)}\Theta_{m}^{(l)})\odot\sigma(\mathbf{X}^{(l)}\Theta_{i}^{(l)}),L,G).$
(5)

This meticulous message passing approach is instrumental in ensuring the
effectiveness of our model. By accurately capturing and refining complex
features, it significantly enhances the model’s ability to make precise
predictions in the context of drug synergy.

###  2.4 Consolidation: Integration and Prediction

In the culminating phase of our methodology, known as the consolidation phase,
we synthesize the previously refined features into predictive insights. This
is achieved through a binary classification model that leverages the enriched
features from the hypergraph. It distinguishes between synergistic and non-
synergistic drug-cell line interactions, providing a binary outcome for each
drug synergy scenario. This phase is crucial for accurately identifying
effective drug combinations.

####  2.4.1 Predictive Model

Our predictive model processes data formatted as tuples
$(\mathbf{x}_{d_{1}},\mathbf{x}_{d_{2}},...,\mathbf{x}_{d_{n}},\mathbf{x}_{c})$,
where $\mathbf{x}_{d_{i}}$ represents a drug and $\mathbf{x}_{c}$ a cell line.
Each $\mathbf{x}$ is an output from the $l^{th}$ layer of a prior Meticulous
Message Passing (see Section 2.3). A key feature of our model is its
capability to interpret these tuples irrespective of the sequence of drugs, as
their order is arbitrary and should not influence the outcome.

To accommodate the unordered nature of the drugs, we utilizes an attention
mechanism, as outlined by Vaswani [22]. This mechanism, prior to position
encoding, inherently does not assume any specific order among inputs, aligning
well with our model’s requirements. To further enhance the model’s
expressiveness, we introduce learnable modality embeddings for both drugs
($\mathbf{m}_{d}$) and cell lines ($\mathbf{m}_{c}$). In our attention-based
model, we define the components for each attention head $j$ as query
($\mathbf{Q}_{j}$), key ($\mathbf{K}_{j}$), and value ($\mathbf{V}_{j}$):

$\mathbf{Q}_{j}=(\mathbf{x}_{c}+\mathbf{m}_{c})^{T}W_{Q}^{j},$ (6)

$\mathbf{K}_{j}=\begin{bmatrix}(\mathbf{x}_{d_{1}}+\mathbf{m}_{d})^{T}W_{K}^{j}\\\
(\mathbf{x}_{d_{2}}+\mathbf{m}_{d})^{T}W_{K}^{j}\\\ \vdots\\\
(\mathbf{x}_{d_{n}}+\mathbf{m}_{d})^{T}W_{K}^{j}\\\
(\mathbf{x}_{c}+\mathbf{m}_{c})^{T}W_{K}^{j}\\\ \end{bmatrix},$ (7)

$\mathbf{V}_{j}=\begin{bmatrix}(\mathbf{x}_{d_{1}}+\mathbf{m}_{d})^{T}W_{V}^{j}\\\
(\mathbf{x}_{d_{2}}+\mathbf{m}_{d})^{T}W_{V}^{j}\\\ \vdots\\\
(\mathbf{x}_{d_{n}}+\mathbf{m}_{d})^{T}W_{V}^{j}\\\
(\mathbf{x}_{c}+\mathbf{m}_{c})^{T}W_{V}^{j}\\\ \end{bmatrix},$ (8)

where $W_{Q}^{j}$, $W_{K}^{j}$, and $W_{V}^{j}$ are respective weight matrices
for the query, key, and value vectors in the $j^{th}$ attention head.

The computation in each attention head is executed as follows:

$\displaystyle\mathbf{head}_{j}$
$\displaystyle=\mathrm{Attention}(\mathbf{Q}_{j},\mathbf{K}_{j},\mathbf{V}_{j})$
$\displaystyle=\mathrm{Softmax}\left(\frac{\mathbf{Q}_{j}\mathbf{K}_{j}^{T}}{\sqrt{d_{k}}}\right)\mathbf{V}_{j},$
(9)

Following this, the model forms a comprehensive representation by
concatenating the outputs from all attention heads:

$\mathbf{h}_{d_{1},d_{2},\ldots,d_{n},c}=\mathrm{Concat}(\mathbf{head}_{1},\mathbf{head}_{2},\ldots,\mathbf{head}_{h})W^{O},$
(10)

Ultimately, a Multi-layer Perceptron (MLP) is applied to this concatenated
output to generate the final predictive output. The MLP processes the
representation and delivers a scalar value, which is normalized to a range
between 0 and 1 using a sigmoid activation function. This normalized output
represents the probability of potential synergistic effects among the drugs:

$\mathbf{s}_{d_{1},d_{2},\cdots,d_{n},c}=\text{MLP}(\mathbf{h}_{d_{1},d_{2},\cdots,d_{n},c}),$
(11)

This structure enables our model to adeptly discerns intricate relationships
between drugs and cell lines, without being influenced by the order in which
the drugs are presented.

####  2.4.2 Training objective

Given the binary nature of our prediction task, we employ the cross-entropy
loss function, a standard and effective choice for binary classification
models. The cross-entropy loss is defined as:

$L_{\mathrm{cross-
entropy}}=-\frac{1}{N}\sum_{i=1}^{N}[y_{i}\log(\hat{y}_{i})+(1-y_{i})\log(1-\hat{y}_{i})],$
(12)

Here, $y_{i}$ represents the actual label (synergistic or non-synergistic) of
the $i$-th sample, and $\hat{y}_{i}$ is the predicted probability generated by
the model for that sample.

##  3 Experiments and Results

###  3.1 Datasets

We collected four categories of data for our research, which include drug
synergy data, molecular information for drugs, genomic characteristics of
cancer cell lines, and indications for drug therapy in diseases from multiple
publicly available databases. Here are the specific details:

  * •

Drug synergy datasets

We gathered data on the synergy of anti-cancer drugs from two prominent large-
scale tumor screening datasets – the O’Neil dataset [10] and the NCI-ALMANAC
dataset [11]. The O’Neil dataset comprises 23,062 samples involving 38 unique
drugs and 39 distinct human cancer cell lines. Each sample measures Loewe
synergy scores for two drugs in combination with a specific cell line. The
NCI-ALMANAC dataset contains 304,549 samples, including ComboScores for 104
FDA-approved drugs in pairings across the NCI-60 cell line panel.

  * •

Drug molecular structures

Information on the simplified molecular input line entry system (SMILES) of
drugs was obtained from the PubChem database [18].

  * •

Gene expression in cancer cell lines

Data on gene expression in cancer cell lines were sourced from the Cell Lines
Project within the COSMIC database [23]. In this context, we specifically
considered the expression values of 651 genes related to the COSMIC cancer
gene census. These expression values were subjected to logarithmic (log2)
transformation and z-score normalization.

We conducted data preprocessing on both the NCI-ALMANAC and O’Neil datasets,
excluding cell lines lacking gene expression data and drugs without SMILES
information. Following this refinement, the NCI-ALMANAC dataset comprises
74,139 measurement samples of ComboScores for 87 drugs across 55 cancer cell
lines, while the O’Neil dataset encompasses 18,950 samples of Loewe synergy
scores for 38 drugs and 39 cancer cell lines. Subsequently, we removed drugs
from the ’drug-disease dataset’ that were not present in the aforementioned
drug synergy datasets. This led to the extraction of indications for 111
diseases corresponding to 111 drugs within the NCI-ALMANAC dataset and 222
diseases corresponding to 222 drugs within the O’Neil dataset. These datasets
are employed in our research study.

Table 1: Comparison of Different Datasets

Dataset | #Drugs | #Cell lines | #Samples
---|---|---|---
NCI-ALMANAC | 87 | 55 | 75,493
O’Neil | 38 | 39 | 18,950

###  3.2 Baselines

In this study, we conducted a comparative analysis of our approach with
several contemporary drug synergy prediction techniques. Here is a brief
overview of each of the baseline methods:

  * •

DeepSynergy [12]: It utilizes a three-layer feedforward neural network to
predict synergy scores, incorporating gene expression as cell line features
and three types of chemical descriptors as drug features.

  * •

DTF [13]: DTF extracts latent features from the drug synergy matrix through
tensor factorization and employs them to train a deep neural network model for
predicting drug synergy.

  * •

HypergraphSynergy [14]: This method formulates synergistic drug combinations
across cancer cell lines as a hypergraph. In this hypergraph, drugs and cell
lines are represented by nodes, while synergistic drug–drug–cell line triplets
are represented by hyperedges. It leverages the biochemical features of drugs
and cell lines as node attributes. Additionally, a hypergraph neural network
is designed to learn drug and cell line embeddings from the hypergraph and
predict drug synergy.

###  3.3 Experiment Setup

To assess the predictive performance of our method in comparison to the
baseline models, we adopted a cross-validation (CV) approach. This involved
randomly partitioning each dataset into a CV set comprising 90% of the samples
and an independent test set containing the remaining 10%. Within the CV set,
we considered three distinct practical scenarios:

1\. Random CV: In this scenario, a five-fold CV was executed by randomly
dividing the samples into five equal folds. The aim here was to rediscover
established anti-cancer drug synergies.

2\. Stratified CV for Cell Lines: We implemented a five-fold CV by randomly
splitting the samples at the cell line level, ensuring that the test set
exclusively contained previously unseen cell lines that were not part of the
training set.

3\. Stratified CV for Drug Combinations: Similar to the previous scenario, we
conducted a five-fold CV by partitioning the samples at the level of drug
combinations. This was done to ensure that the test set exclusively comprised
previously unencountered drug combinations not present in the training set.

Additionally, we employed an independent test set to evaluate the performance
of our model and the baseline models. All methods were trained on the entire
CV set and then used for making predictions on the independent test set,
facilitating an objective evaluation.

For classification tasks, we transformed synergy scores into binary outcomes
using a threshold established in prior studies [12, 13]. Samples with synergy
scores exceeding 30 were categorized as positive, while those below this
threshold were considered negative.

###  3.4 Performance Comparison and Analysis

For the sake of evaluating the performance of VQSynergy, we compare it with
baselines mentioned above. The comparative results of different methods for
drug synergy prediction, as summarized in Table 3, distinctly highlights the
superiority of the VQSynergy model across various testing scenarios. In the
Random (Rand.) scenario, VQSynergy achieves the highest performance with a
score of 0.8566, which is a noticeable improvement over its closest
competitor, HypergraphSynergy, at 0.8447. This represents a performance
increase of approximately 1.41%. In the Stratified Cell scenario, VQSynergy
again leads with a score of 0.7808, surpassing HypergraphSynergy’s 0.7752 by a
narrower margin of 0.56%. The most significant performance differential is
observed in the Stratified DrugComb scenario, where VQSynergy’s score of
0.8002 outperforms HypergraphSynergy’s 0.7745 by 2.57%.

The overall analysis indicates that VQSynergy consistently outperforms other
methods like DFT, DeepSynergy, and HypergraphSynergy. This consistent lead
across different evaluation metrics underscores the robustness and
effectiveness of the VQSynergy model. It also suggests that the incorporation
of the Vector Quantization mechanism provides a substantial advantage in
accurately predicting drug synergies, especially in complex scenarios
represented by Stratified DrugComb. The performance of VQSynergy in these
various scenarios not only demonstrates its efficacy but also highlights its
potential as a reliable tool in the field of drug synergy prediction, capable
of handling diverse and complex datasets with higher accuracy compared to
existing methods.

Table 2: Comparison of Different Methods (NCI-ALMANAC)

Method | Rand. | Strat. CellLine | Strat. DrugComb
---|---|---|---
DFT | 0.8238 | 0.7551 | 0.7502
DeepSynergy | 0.8350 | 0.7035 | 0.7739
HypergraphSynergy | 0.8447 | 0.7752 | 0.7745
VQSynergy | 0.8566 | 0.7808 | 0.8002

Table 3: Comparison of Different Methods (O’Neil)

Method | Rand. | Strat. CellLine | Strat. DrugComb
---|---|---|---
DFT | 0.9138 | 0.8387 | 0.8424
DeepSynergy | 0.9060 | 0.8188 | 0.8583
HypergraphSynergy | 0.9230 | 0.8554 | 0.8621
VQSynergy | 0.9232 | 0.8611 | 0.8715

###  3.5 Robustness Validation

In this part, our experiment is structured to compare the performance of
VQSynergy under two main scenarios: without Gaussian noise and with Gaussian
noise (standard deviation, $\sigma$ = 0.1). We evaluate the model in three
different settings. The aim is to assess the effectiveness of the VQ mechanism
in maintaining the model’s predictive accuracy under noisy conditions.

Initially, the model was tested without the introduction of Gaussian noise. In
this setting, as shown in Table 4, the performance of VQSynergy without the VQ
mechanism is marginally better than when the VQ mechanism is incorporated. The
difference in performance is minimal, suggesting that under normal conditions
without external noise, the addition of the VQ mechanism does not
significantly enhance the model’s capability.

The robustness of the VQSynergy model becomes particularly evident with the
introduction of Gaussian noise ($\sigma$ = 0.1), as shown in Table 5. When the
VQ mechanism is employed, there is a noticeable improvement in performance
across all three scenarios. Specifically, the use of VQ results in an increase
of 2.52% for Random (82.19% vs. 79.67%), 0.64% for stratified cell line
(75.90% vs. 75.26%), and 1.60% for stratified drug combination (76.77% vs.
75.17%) compared to the model without VQ mechanism. These results clearly
illustrate the effectiveness of the VQ mechanism in enhancing the model’s
performance, particularly under conditions with data imperfections like
Gaussian noise.

Table 4: Results w/o Gaussian Noise

Method | Rand. | Strat. CellLine | Strat. DrugComb
---|---|---|---
VQSynergy (w/o VQ) | 0.8591 | 0.7836 | 0.8014
VQSynergy (w/ VQ) | 0.8566 | 0.7808 | 0.8002

Table 5: Results in Gaussian Noise ($\sigma$ = 0.1)

Method | Rand. | Strat. CellLine | Strat. DrugComb
---|---|---|---
VQSynergy (w/o VQ) | 0.7967 | 0.7526 | 0.7517
VQSynergy (w/ VQ) | 0.8219 | 0.7590 | 0.7677

##  4 Conclusion and Future Direction

In this study, we introduce VQSynergy, a robust framework employing the Vector
Quantization (VQ) mechanism for enhanced drug synergy prediction. This
innovative model combines VQ with gated residuals and a bespoke attention
mechanism, setting a new standard in accuracy and generalizability in
predictive modeling. Notably, robustness tests reveal that VQSynergy maintains
superior predictive performance even under Gaussian noise, a testament to its
resilience. This demonstrates the distinct advantage of incorporating the VQ
mechanism in handling complex, noisy datasets typical in drug synergy studies.

The emerging trend towards higher-order combination therapies, as evidenced by
the shift from traditional dual-drug approaches to more complex multi-drug
regimens, presents an exciting frontier for future research [24, 25, 26]. This
evolution, marked by the approval and investigation of therapies involving
three or more drugs in diseases like cancer, HIV, and tuberculosis, signifies
a rich area for exploration [27]. The success of regimens like the five-drug
R-CHOP therapy for Diffuse Large B-Cell Lymphoma underscores the potential of
such approaches [27]. In line with this progression, our future work will
attempt to extend the VQSynergy model to analyze and predict the efficacy of
these complex multi-drug combinations. The burgeoning repository of
dose–response data for numerous drug combinations will be instrumental in
advancing this research, paving the way for more effective and tailored
therapeutic strategies [28].

Besides, the conventional VQ methodology, while foundational, is somewhat
rudimentary and not without potential drawbacks. Specifically, the traditional
’arg min’ approach in VQ could introduce issues in optimization and
representation. To address this, an enhancement inspired by recent
advancements, such as the method proposed in [29], can be considered. This
enhancement transforms the distance calculation into a probability
distribution function and incorporates entropy as part of the loss function in
this setup. Furthermore, the straight-through [30] method commonly used in VQ
has been gradually superseded by the Gumbel-Softmax [31, 32] approach.
However, even more advanced gradient estimation techniques, as suggested by
recent research [33, 34, 35], remain unexplored and warrant further
investigation. These methods can be attempt to offer improvements in
optimization and representation fidelity, thus enhancing the overall
performance of VQ-based models in various applications.

## References

  * [1] Jonathan B Fitzgerald, Birgit Schoeberl, Ulrik B Nielsen, and Peter K Sorger.  Systems biology and combination therapy in the quest for clinical efficacy.  Nature chemical biology, 2(9):458–466, 2006.
  * [2] Andrew L Hopkins.  Network pharmacology: the next paradigm in drug discovery.  Nature chemical biology, 4(11):682–690, 2008.
  * [3] Michael M Gottesman, Tito Fojo, and Susan E Bates.  Multidrug resistance in cancer: role of atp–dependent transporters.  Nature reviews cancer, 2(1):48–58, 2002.
  * [4] Ayman Khdair, Di Chen, Yogesh Patil, Linan Ma, Q Ping Dou, Malathy PV Shekhar, and Jayanth Panyam.  Nanoparticle-mediated combination chemotherapy and photodynamic therapy overcomes tumor drug resistance.  Journal of Controlled Release, 141(2):137–144, 2010.
  * [5] Caitriona Holohan, Sandra Van Schaeybroeck, Daniel B Longley, and Patrick G Johnston.  Cancer drug resistance: an evolving paradigm.  Nature Reviews Cancer, 13(10):714–726, 2013.
  * [6] Bissan Al-Lazikani, Udai Banerji, and Paul Workman.  Combinatorial drug therapy for cancer in the post-genomic era.  Nature biotechnology, 30(7):679–692, 2012.
  * [7] Reza Bayat Mokhtari, Sushil Kumar, Syed S Islam, Mehrdad Yazdanpanah, Khosrow Adeli, Ernest Cutz, and Herman Yeger.  Combination of carbonic anhydrase inhibitor, acetazolamide, and sulforaphane, reduces the viability and growth of bronchial carcinoid cell lines.  BMC cancer, 13(1):1–18, 2013.
  * [8] Xiaochen Sun, Santiago Vilar, and Nicholas P Tatonetti.  High-throughput methods for combinatorial drug discovery.  Science translational medicine, 5(205):205rv1–205rv1, 2013.
  * [9] Liye He, Evgeny Kulesskiy, Jani Saarela, Laura Turunen, Krister Wennerberg, Tero Aittokallio, and Jing Tang.  Methods for high-throughput drug combination screening and synergy scoring.  Cancer systems biology: methods and protocols, pages 351–398, 2018.
  * [10] Jennifer O’Neil, Yair Benita, Igor Feldman, Melissa Chenard, Brian Roberts, Yaping Liu, Jing Li, Astrid Kral, Serguei Lejnine, Andrey Loboda, et al.  An unbiased oncology compound screen to identify novel combination strategies.  Molecular cancer therapeutics, 15(6):1155–1162, 2016.
  * [11] Susan L Holbeck, Richard Camalier, James A Crowell, Jeevan Prasaad Govindharajulu, Melinda Hollingshead, Lawrence W Anderson, Eric Polley, Larry Rubinstein, Apurva Srivastava, Deborah Wilsker, et al.  The national cancer institute almanac: a comprehensive screening resource for the detection of anticancer drug pairs with enhanced therapeutic activity.  Cancer research, 77(13):3564–3576, 2017.
  * [12] Kristina Preuer, Richard PI Lewis, Sepp Hochreiter, Andreas Bender, Krishna C Bulusu, and Günter Klambauer.  Deepsynergy: predicting anti-cancer drug synergy with deep learning.  Bioinformatics, 34(9):1538–1546, 2018.
  * [13] Zexuan Sun, Shujun Huang, Peiran Jiang, and Pingzhao Hu.  Dtf: deep tensor factorization for predicting anticancer drug synergy.  Bioinformatics, 36(16):4483–4489, 2020.
  * [14] Xuan Liu, Congzhi Song, Shichao Liu, Menglu Li, Xionghui Zhou, and Wen Zhang.  Multi-way relation-enhanced hypergraph representation learning for anti-cancer drug synergy prediction.  Bioinformatics, 38(20):4782–4789, 2022.
  * [15] Yifan Feng, Haoxuan You, Zizhao Zhang, Rongrong Ji, and Yue Gao.  Hypergraph neural networks.  In Proceedings of the AAAI conference on artificial intelligence, volume 33, pages 3558–3565, 2019.
  * [16] Aaron Van Den Oord, Oriol Vinyals, et al.  Neural discrete representation learning.  Advances in neural information processing systems, 30, 2017.
  * [17] Dianbo Liu, Alex M Lamb, Kenji Kawaguchi, Anirudh Goyal ALIAS PARTH GOYAL, Chen Sun, Michael C Mozer, and Yoshua Bengio.  Discrete-valued neural communication.  Advances in Neural Information Processing Systems, 34:2109–2121, 2021.
  * [18] Sunghwan Kim, Jie Chen, Tiejun Cheng, Asta Gindulyte, Jia He, Siqian He, Qingliang Li, Benjamin A Shoemaker, Paul A Thiessen, Bo Yu, et al.  Pubchem 2019 update: improved access to chemical data.  Nucleic acids research, 47(D1):D1102–D1109, 2019.
  * [19] Yunsheng Shi, Zhengjie Huang, Shikun Feng, Hui Zhong, Wenjin Wang, and Yu Sun.  Masked label prediction: Unified message passing model for semi-supervised classification.  arXiv preprint arXiv:2009.03509, 2020.
  * [20] Qimai Li, Zhichao Han, and Xiao-Ming Wu.  Deeper insights into graph convolutional networks for semi-supervised learning.  In Proceedings of the AAAI conference on artificial intelligence, volume 32, 2018.
  * [21] Song Bai, Feihu Zhang, and Philip HS Torr.  Hypergraph convolution and hypergraph attention.  Pattern Recognition, 110:107637, 2021.
  * [22] Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N Gomez, Łukasz Kaiser, and Illia Polosukhin.  Attention is all you need.  Advances in neural information processing systems, 30, 2017.
  * [23] Simon A Forbes, David Beare, Prasad Gunasekaran, Kenric Leung, Nidhi Bindal, Harry Boutselakis, Minjie Ding, Sally Bamford, Charlotte Cole, Sari Ward, et al.  Cosmic: exploring the world’s knowledge of somatic mutations in human cancer.  Nucleic acids research, 43(D1):D805–D811, 2015.
  * [24] Geraint Davies, Martin Boeree, Dave Hermann, and Michael Hoelscher.  Accelerating the transition of new tuberculosis drug combinations from phase ii to phase iii trials: New technologies and innovative designs.  PLoS medicine, 16(7):e1002851, 2019.
  * [25] Martina Properzi, Paola Magro, Francesco Castelli, and Eugenia Quiros-Roldan.  Dolutegravir-rilpivirine: first 2-drug regimen for hiv-positive adults.  Expert Review of Anti-Infective Therapy, 16(12):877–887, 2018.
  * [26] Philip Gotwals, Scott Cameron, Daniela Cipolletta, Viviana Cremasco, Adam Crystal, Becker Hewes, Britta Mueller, Sonia Quaratino, Catherine Sabatos-Peyton, Lilli Petruzzelli, et al.  Prospects for combining targeted and conventional cancer therapy with immunotherapy.  Nature Reviews Cancer, 17(5):286–301, 2017.
  * [27] Ritu Lakhtakia and Ikram Burney.  A historical tale of two lymphomas: part ii: non-hodgkin lymphoma.  Sultan Qaboos University Medical Journal, 15(3):e317, 2015.
  * [28] Elif Tekin, Cynthia White, Tina Manzhu Kang, Nina Singh, Mauricio Cruz-Loya, Robert Damoiseaux, Van M Savage, and Pamela J Yeh.  Prevalence and patterns of higher-order drug interactions in escherichia coli.  NPJ systems biology and applications, 4(1):31, 2018.
  * [29] Aditya Ramesh, Mikhail Pavlov, Gabriel Goh, Scott Gray, Chelsea Voss, Alec Radford, Mark Chen, and Ilya Sutskever.  Zero-shot text-to-image generation.  In International Conference on Machine Learning, pages 8821–8831. PMLR, 2021.
  * [30] Yoshua Bengio, Nicholas Léonard, and Aaron Courville.  Estimating or propagating gradients through stochastic neurons for conditional computation.  arXiv preprint arXiv:1308.3432, 2013.
  * [31] Eric Jang, Shixiang Gu, and Ben Poole.  Categorical reparameterization with gumbel-softmax.  arXiv preprint arXiv:1611.01144, 2016.
  * [32] Chris J Maddison, Andriy Mnih, and Yee Whye Teh.  The concrete distribution: A continuous relaxation of discrete random variables.  arXiv preprint arXiv:1611.00712, 2016.
  * [33] Max B Paulus, Chris J Maddison, and Andreas Krause.  Rao-blackwellizing the straight-through gumbel-softmax gradient estimator.  arXiv preprint arXiv:2010.04838, 2020.
  * [34] Ting-Han Fan, Ta-Chung Chi, Alexander I Rudnicky, and Peter J Ramadge.  Training discrete deep generative models via gapped straight-through estimator.  In International Conference on Machine Learning, pages 6059–6073. PMLR, 2022.
  * [35] Liyuan Liu, Chengyu Dong, Xiaodong Liu, Bin Yu, and Jianfeng Gao.  Bridging discrete and backpropagation: Straight-through and beyond.  arXiv preprint arXiv:2304.08612, 2023.


# Dynamical phase transition in the occupation fraction statistics for non-
crossing Brownian particles

Soheli Mukherjee<EMAIL_ADDRESS>Department of Solar Energy and
Environmental Physics, Blaustein Institutes for Desert Research, Ben-Gurion
University of the Negev, Sede Boqer Campus, 8499000, Israel  Naftali R. Smith
<EMAIL_ADDRESS>Department of Solar Energy and Environmental Physics,
Blaustein Institutes for Desert Research, Ben-Gurion University of the Negev,
Sede Boqer Campus, 8499000, Israel

###### Abstract

We consider a system of $N$ non-crossing Brownian particles in one dimension.
We find the exact rate function that describes the long-time large deviation
statistics of their occupation fraction in a finite interval in space.
Remarkably, we find that, for any general $N\geq 2$, the system undergoes
$N-1$ dynamical phase transitions of second order. The $N-1$ transitions are
the boundaries of $N$ phases that correspond to different numbers of particles
which are in the vicinity of the interval throughout the dynamics. We achieve
this by mapping the problem to that of finding the ground-state energy for $N$
noninteracting spinless fermions in a square-well potential. The phases
correspond to different numbers of single-body bound states for the quantum
problem. We also study the process conditioned on a given occupation fraction
and the large-$N$ limiting behavior.

##  I Introduction

Fluctuations in stochastic systems are of central importance in statistical
mechanics and other fields fluctuations ; fluctuations1 ; OvMe2010 ;
MeersonAssaf2017 ; neqpt ; currentfluc ; currentfluc1 ; condensation .
Examples include population dynamics OvMe2010 ; MeersonAssaf2017 , non-
equilibrium phase transitions neqpt , current statistics currentfluc ;
currentfluc1 , condensation phenomena condensation , etc to name a few. In
particular large deviations (or rare events) have been a central theme of
interest over the past few decades ellis ; hugo2009 ; Varadhan ; O1989 ; DZ ;
Hollander ; Jack20 . The fluctuations of the quantities of interest are
encoded in the large deviation functions (LDFs), which are analogous to
thermodynamic potentials in equilibrium systems hugo2009 . One of the most
remarkable phenomena in the context of large deviations is dynamical phase
transitions (DPTs), which correspond to singularities in LDFs exclusion ;
exclusion1 ; baek ; baek1 ; kafri ; glass ; singularities .

The occupation time of a particle is the time that it spends in some specified
spatial domain. Fluctuations of the occupation time have been of interest in
many random systems, including coarsening and phase ordering dynamics in
magnetic systems phaseorder ; phaseorder1 ; phaseorder2 , financial time
series finance , random walks on graphs randomwalk ; randomwalk1 . Since the
seminal work of Lévy levy , who computed the exact probability distribution of
the occupation time for an ordinary Brownian motion, there have been studies
of the fluctuations of occupation time for various non-equilibrium processes
kac ; kac2 ; godreche ; majumdar1 ; TouchetteMinimalModel ;
Touchetteoneparticle ; AKM19 ; KA22 , because of their potential applications
in many physical systems. For example, to analyze the morphological dynamics
of interfaces morphology , the fluorescence intermittency emitting from
colloidal semiconductor dots semidots , optical imaging imaging  etc.

Recently, Touchette and Nyawo studied fluctuations of occupation fraction in a
closed interval for a Brownian motion with and without drift
TouchetteMinimalModel ; Touchetteoneparticle . They demonstrated that this
relatively simple model exhibits a DPT in the fluctuations of the occupation
time in presence of a drift in the long-time limit. Their analysis of this
model used the Donsker-Varadhan (DV) large-deviation formalism hugo2009 ;
donsker ; donsker1 ; majumdar4 ; majumdar5 ; Jack20 , which maps the problem
onto that of finding the ground-state energy of a quantum problem of a single
particle inside an effective potential well. The transition occurs between the
escape and confinement of the Brownian motion and is first order in nature.
Similar DPTs have been observed in many systems in the limit of low noise
and/or large system size exclusion ; exclusion1 ; baek ; baek1 ; kafri ; glass
; singularities ; fluctuations1 ; lownoise ; naftali2 . To observe the DPT in
drifted Brownian motion occupation fraction it is sufficient to consider the
long-time limit without taking any additional limits of small or large
parameters TouchetteMinimalModel ; Touchetteoneparticle .

It is natural to ask how the occupation fraction fluctuates for $N>1$ Brownian
particles, conditioned not to cross each other. The non-crossing condition
introduces correlations between the particles, making this a nontrivial, many-
body problem. Since the pioneering works genne  and fishersurvival , non-
crossing random walkers have been studied in the context of wetting and
melting fishersurvival , networks of polymers polymers , persistence
properties in nonequilibrium systems viciouswalker1  and more. Non-crossing
Brownian motions are also known to be closely related to random matrix theory
because the joint distribution of their positions (which has been studied in
many contexts schehr1 ; schehr ; schehr2 ; NIBM ; GLMS19 ; fermions )
coincides, in some special cases, to that of eigenvalues of certain types of
random matrices Dyson62 ; rmt ; rmt1 .

In this work, we extend the study of occupation fraction statistics to $N>1$
non-crossing Brownian particles in an interval $[-l,l]$ for a given time
interval $[0,T]$. By extending the DV formalism to non-crossing Brownian
motions, we study the large deviation function for all $N\geq 1$. We find that
the DV formalism maps this problem to $N$ noninteracting, spinless _fermions_
in a square-well potential, cf. GLMS19 ; fermions . Interestingly, we find
that the large deviation functions show multiple singularities: The system
undergoes $N-1$ second order DPTs. These transitions are of very different
nature to the DPT found for a single Brownian particle in
TouchetteMinimalModel ; Touchetteoneparticle . In particular, they occur in
the absence of a drift, i.e., the dynamics obey time-reversal symmetry. In
each of the $N$ different phases, separated by the DPTs, a different number of
particles remains in the vicinity of the interval for the entire dynamics,
while the other particles wander away from the interval.

The paper is organized as follows. In Sec. II we define the model and present
the scaling behavior of the fluctuations of the occupation fraction. In Sec.
III we show how to extend the Donsker-Varadhan (DV) formalism to our system.
In Sec. IV, we solve the DV problem and show that there are $N-1$ second order
phase transitions. We give explicit results for $N=2$. We also discuss the
process conditioned on a given value of the occupation fraction, and some
limiting behaviors that emerge in the limit $N\gg 1$. Finally, we summarize
and conclude in Sec. V. Some details of the calculations are given in the
Appendices.

##  II Model

We consider $N$ one-dimensional Brownian particles, which are defined by the
following stochastic differential equations

$\dot{X}_{i}(t)=\sqrt{2D}\;\xi_{i}(t),\qquad i=1,2,\dots,N,$ (1)

where $X_{i}$ is the position of the $i_{th}$ particle, $D$ is the diffusion
constant and $\xi_{i}(t)$ are Gaussian white noises with
$\langle\xi_{i}(t)\rangle=0$ and
$\langle\xi_{i}(t)\xi_{j}(t^{\prime})\rangle=\delta_{ij}\delta(t-t^{\prime})$.
Here $\langle.\rangle$ denotes the ensemble average over realizations of the
noise. We condition the particles to be non-crossing, i.e,
$X_{1}(t)<X_{2}(t)<\dots<X_{N}(t)$ for all $t\in[0,T]$.

The occupation fraction is defined as

$\rho_{T}=\frac{1}{T}\int_{0}^{T}\Big{(}\mathds{1}_{[-l,l]}(X_{1}(t))+\dots...+\mathds{1}_{[-l,l]}(X_{N}(t))\Big{)}\,dt,$
(2)

where $\mathds{1}_{[\ -l,l]}(x)$ is the indicator function

$\mathds{1}_{[\
-l,l]}(x)=\begin{cases}1,&\text{for}\,\,x\in[-l,l],\\\\[5.69054pt]
0,&\text{otherwise.}\end{cases}$ (3)

$\rho_{T}$ is a random variable that takes values on $[0,N]$, where
$\rho_{T}=0$ ($\rho_{T}=N$) corresponds to realizations in which none (all) of
the Brownian particles stay in the interval $[-l,l]$ for the entire duration
of the dynamics. For $T\longrightarrow\infty$, it becomes very unlikely for
the Brownian particles to spend an extensive time $O(T)$ in the interval
$[-l,l]$, so we expect that the probability distribution $P(\rho_{T}=\rho)$
concentrates around $\rho=0$.

As we show below, in the long-$T$ limit, the fluctuations of $\rho_{T}$ follow
the large deviation principle ellis ; Varadhan ; O1989 ; DZ ; Hollander ;
hugo2009 ; Jack20 ; donsker ; donsker1 ; gartner ; ellis1 , which states that
the probability distribution $P(\rho_{T}=\rho)$ scales as

$P(\rho_{T}=\rho)\sim e^{-TI(\rho)},$ (4)

describing an exponential decay with $T$ at a rate described by the LDF or the
rate function $I(\rho)\geq 0$, which is defined as

$I(\rho)=-\lim_{T\longrightarrow\infty}\frac{1}{T}\ln P(\rho_{T}=\rho)$ (5)

In general, the rate function is difficult to calculate directly. However,
according to the Gärtner-Ellis theorem gartner ; ellis1  the rate function
$I(\rho)$ can be calculated via the Legendre-Fenchel transform of the scaled
cumulant generating function (SCGF), which is defined as hugo2009 ; gartner ;
ellis1

$\lambda(k)=\lim_{T\longrightarrow\infty}\frac{1}{T}\,\,\ln\langle
e^{Tk\rho_{T}}\rangle$ (6)

The rate function is then expressed as the Legendre-Fenchel transformation of
$\lambda(k)$

$I(\rho)=\sup_{k}\\{k\,\rho-\lambda(k)\\}$ (7)

provided that $\lambda(k)$ exists and differentiable hugo2009 . If
$\lambda(k)$ is convex, then the transform reduces to the Legendre transform,
so that the conjugate variable of $k$, $\rho$ is given by LegendreNutshell

$\rho=\frac{d\lambda(k)}{dk}\,.$ (8)

Hence, the problem reduces to that of calculating the SCGF of the occupation
fraction of $N$ non-crossing Brownian particles. In the next section we
explain how this can be done, by extending the DV formalism to noncrossing
Brownian particles.

##  III Donsker-Varadhan formalism

An equivalent description of the dynamics of the $N$ Brownian particles is
given in terms of the Fokker-Planck equation

$\frac{\partial Q(X_{1},\dots,X_{N},t)}{\partial
t}=L^{\dagger}\,\,Q(X_{1},\dots,X_{N},t)\,.$ (9)

Here $Q(X_{1},\dots,X_{N},t)$ is the unnormalized time dependent joint
probability density function of the $X_{i}$’s, $L^{\dagger}=D\nabla^{2}$ is
called the Fokker-Planck generator, with
$\nabla^{2}=\frac{\partial^{2}}{\partial
X_{1}^{2}}+\dots+\frac{\partial^{2}}{\partial X_{N}^{2}}$ . $L^{\dagger}$ is a
Hermitian operator ($L^{\dagger}=L$) as the dynamics have time-reversal
symmetry LegendreNutshell . Due to the non-crossing condition, we consider
only the domain $X_{1}<\dots<X_{N}$, with the Dirichlet boundary conditions

$Q(X_{1},\dots,X_{N},t)\Big{|}_{X_{i}=X_{i+1}}\\!\\!\\!=0\,\,\,\text{for
any}\;i\\!\in\\!\\{1,\dots,N-1\\}.$ (10)

The _normalized_ joint probability density function of $X_{1},\dots,X_{N}$ is

$P(X_{1},\dots,X_{N},t)=\frac{Q(X_{1},\dots,X_{N},t)}{\mathcal{N}(t)}$ (11)

where $\mathcal{N}(t)$ is the (time-dependent) normalization constant

$\mathcal{N}(t)=\int_{-\infty}^{\infty}\\!\\!dX_{1}\dots\int_{-\infty}^{\infty}\\!\\!dX_{N}Q(X_{1},\dots,X_{N},t).$
(12)

At long times $t\longrightarrow\infty$, $\mathcal{N}(t)$ goes as a power law
$\mathcal{N}(t)\sim t^{-N(N-1)/4}$ fishersurvival ; fishersurvival1 ;
viciouswalker ; viciouswalker1 ; GLMS19 .

A mathematical framework for calculating large deviations of time-averaged
quantities in stochastic processes of the form Eq. (2) (which are often
referred to as dynamical or additive observables) was formulated by Donsker
and Varadhan in donsker ; donsker1 . We will now extend this formalism to the
case of noncrossing Brownian particles. In order to calculate $\lambda(k)$,
the Feynman-Kac formula is used kac ; Majumdar2007 , which states that the
evolution of the generating function of $\rho_{T}$ involves a linear operator,
called tilted generator $\mathcal{L}_{k}^{\dagger}$. In our case the
expression of the $\mathcal{L}_{k}^{\dagger}$ is

$\mathcal{L}_{k}^{\dagger}=D\,\,\nabla^{2}+k\Big{(}\mathds{1}_{[-l,l]}(X_{1})+\dots+\mathds{1}_{[-l,l]}(X_{N})\Big{)}\,.$
(13)

The SCGF of $\rho_{T}$ is evaluated as the largest eigenvalue of
$\mathcal{L}_{k}$ Touchette2018  (see Appendix A for more detailed
calculation)

The eigenvalue problem is defined by the equation

$\mathcal{L}_{k}\,\,r_{k}(X_{1},\dots,X_{N})=\lambda(k)\,\,r_{k}(X_{1},\dots,X_{N})$
(14)

where $r_{k}(X_{1},\dots,X_{N})$ is the right eigenfunction of
$\mathcal{L}_{k}$, together with the same boundary conditions as for $Q$,
i.e.,

$r(X_{1},\dots,X_{N},t)\Big{|}_{X_{i}=X_{i+1}}\\!\\!=0\,\,\,\text{for
any}\,\,\,i\in\\{1,\dots,N-1\\}.$ (15)

Note that as the operator $\mathcal{L}_{k}$ is Hermitian, its left
eigenfunctions are the same as $r_{k}(X_{1},\dots,X_{N})$. We are generally
interested in normalizable solutions, and we (arbitrarily) choose the
normalization constant to be 1

$\int_{-\infty}^{\infty}\\!dX_{1}\dots\int_{-\infty}^{\infty}\\!dX_{N}\,|r_{k}(X_{1},\dots,X_{N})|^{2}=1\,.$
(16)

Figure 1: (Color online) The effective quantum problem of $N$ fermions in a
square potential well (17) of depth $-k$, whose solution yields the behavior
of fluctuations of the occupation fraction in the interval $[-l,l]$ of $N$
non-crossing Brownian particles. The solid red lines represent the single-
particle energy levels occupied by the fermions. $-\lambda(k)$ is the many-
body ground-state energy.

Eq. (14) with the tilted operator (13) and the boundary condition (15) can be
recognized as the time-independent Schrödinger equation of $N$ noninteracting,
spinless fermions up to an overall minus sign with $\frac{\hslash^{2}}{2m}=D$
in a finite potential well

$V_{k}(x)=-k\mathds{1}_{[-l,l]}(x)\,.$ (17)

Fig. 1 shows a schematic plot of the equivalent quantum problem. This implies
that the $\lambda(k)$ equals minus the (many-body) ground state energy of
$-\mathcal{L}_{k}$. For $N=1$ the problem reduces to that of a single particle
in the potential well $V_{k}$, which was solved in TouchetteMinimalModel ;
Touchetteoneparticle .

##  IV Results

The SCGF $\lambda(k)$ for general $N$ is calculated from the ground state
energy of $N$ noninteracting, spinless fermions in a square well potential of
depth $k$, which is given by the sum of the single-body energy levels. For a
single Brownian particle ($N=1$), there is always a bound state for any $k>0$
given by cosine inside the well and two decaying exponentials outside the
well. The single-body energy levels of a quantum particle inside a square-well
potential can be found in any standard textbook of quantum physics hall ; qm ,
and each of them is given by the solution to one of two transcendental
equations (see below). For any general $N$, an interesting effect occurs: As
$k$ is increased, the number of solutions to the transcendental equations
increases, and, as shown below, this results in $N-1$ dynamical phase
transitions which correspond to critical values of $k$.

In this section, we first briefly recall the single-body energy levels of a
quantum particle in a finite potential well. Then, we give a brief review of
the results for the case $N=1$ TouchetteMinimalModel ; Touchetteoneparticle .
Next, we study the problem for general $N$, giving explicit results for the
case $N=2$ which is already sufficient to observe the dynamical phase
transition. We then consider the large-$N$ limit where we uncover a universal
behaviour of the system. We also discuss the conditioned process, both for
general $N$ and in the large-$N$ limit. We choose the units of space and time
such that $D=\frac{1}{2}$, $l=1$ throughout the entire section.

###  IV.1 Single-body energy levels

Due to the symmetry in the potential Eq. (17), the single-body wave functions
of bound states are either symmetric or antisymmetric. For any $k>0$, there
always exists at least one bound state with a symmetric wave function hall .
As the depth of the potential $k$ increases, the number of bound energy levels
increases, with the wave functions alternating between symmetric and
antisymmetric. The energy levels $-\lambda_{1},-\lambda_{2},\dots$ satisfy the
following transcendental equations hall ; qm

$A+A\tan^{2}\left(\sqrt{2A}\right)=k\qquad\text{for symmetric case}$ (18)

and

$A+A\cot^{2}\left(\sqrt{2A}\right)=k\qquad\text{for antisymmetric case}$ (19)

where

$A=k-\lambda_{i}(k).$ (20)

The $i^{\text{th}}$ energy level ($i=1,2,\dots$) is the solution of Eq. (18)
(for odd $i$) or Eq. (19) (for even $i$) for the following range of $A$’s

$\frac{\pi^{2}}{8}(i-1)^{2}<A_{i}<\frac{\pi^{2}}{8}i^{2}\,.$ (21)

###  IV.2 $N=1$

For a single ($N=1$) Brownian motion the problem was solved in
TouchetteMinimalModel ; Touchetteoneparticle . The calculation of
$\lambda(k)=\lambda_{1}(k)$ boils down to calculating the ground state energy
of one quantum particle trapped in a potential well with depth $k$. The ground
state wave function is symmetric satisfying Eq. (18), where $0\leq k<\infty$
and $0\leq A\leq\frac{\pi^{2}}{8}$. Eqs. (18) and (20) thus give $\lambda(k)$
in a parametric form (with $k$ and $\lambda$ both given explicitly as
functions of the parameter $A$).

The conjugate variable $\rho$ and the rate function $I(\rho)$ of the problem
can also be expressed in terms of the variable $A$, thus giving the rate
function in a parametric form:

$\displaystyle\rho$ $\displaystyle=$
$\displaystyle\frac{d\lambda}{dk}=\frac{d\lambda/dA}{dk/dA}=1-\frac{\cos^{2}\left(\sqrt{2A}\right)}{1+\sqrt{2A}\tan\left(\sqrt{2A}\right)}\,,$
(22) $\displaystyle I$ $\displaystyle=$
$\displaystyle\rho\,k-\lambda=\frac{A\sqrt{2A}\tan\left(\sqrt{2A}\right)}{\sqrt{2A}\tan\left(\sqrt{2A}\right)+1}\,.$
(23)

The asymptotic behaviours of $\lambda(k)$ and $I(\rho)$ are given by
Touchetteoneparticle  (see also Appendix B for a derivation)

$\displaystyle\lambda\left(k\right)$ $\displaystyle\simeq$
$\displaystyle\begin{cases}2k^{2},&k\longrightarrow 0,\\\\[2.84526pt]
k-\frac{\pi^{2}}{8}+\frac{\pi^{2}}{4\sqrt{2k}},&k\longrightarrow\infty,\end{cases}$
(24) $\displaystyle I\left(\rho\right)$ $\displaystyle\simeq$
$\displaystyle\begin{cases}\frac{\rho^{2}}{8},&\rho\longrightarrow
0,\\\\[2.84526pt]
\frac{\pi^{2}}{8}-3\left(\frac{\pi^{2}}{8\sqrt{2}}\right)^{2/3}\left(1-\rho\right)^{1/3},&\rho\longrightarrow
1.\end{cases}$ (25)

The SCGF and rate function for $N=1$ are shown in Figs. 2 and 3 respectively
($\lambda_{1}(k)$ and $I_{1}(\rho)$ respectively).

Figure 2: (Color online) Solid line: The plot of the SCGF $\lambda(k)$ as a
function of $k$ for $N=2$ non-crossing Brownian particle occupation fraction.
The solid circle represents the critical point $k_{c}=\pi^{2}/8$. For
$k<k_{c}$, $\lambda(k)$ coincides with the SCGF $\lambda_{1}(k)$ for $N=1$,
and for $k>k_{c}$, the SCGF is $\lambda(k)=\lambda_{1}(k)+\lambda_{2}(k)$
where $-\lambda_{2}(k)$ is the energy of the first excited state for the
potential (17). The dot-dashed line represents the continuation of the single
particle SCGF for $k>k_{c}$. The red dashed line depicts the asymptotic
behavior of $\lambda(k)$ for $k\gg 1$. The second derivative of $\lambda(k)$
has a jump discontinuity at $k=k_{c}$. This can be seen as a corner
singularity in the inset which shows $d\lambda/dk$ as a function of $k$. Units
are chosen such that $D=\frac{1}{2}$ and $l=1$.

###  IV.3 General $N$

For $N>1$, the problem is mapped to $N$ noninteracting, spinless fermions in
the potential well (17). As $k$ is increased, the number of single-body bound
states increases and the fermions alternately occupy the odd and even energy
levels which are found from Eqs. (18) and (19), respectively. These energy
levels correspond to multiple values of $A$ with ranges given by Eq. (21).

Importantly, there are $N-1$ critical values of $k=k_{c,i}\equiv
i^{2}\frac{\pi^{2}}{8}$ for $i\in\\{1,\dots,N-1\\}$, at which the number of
bound states increases. The critical values correspond to singularities of
$\lambda(k)$ and these correspond to dynamical phase transitions, i.e.,
singularities of $I(\rho)$. Indeed, at $k$ slightly larger than $k_{c,i}$,
$\lambda_{i+1}(k)$ behaves asymptotically as (for details see Appendix B)

$\lambda_{i+1}(k)\simeq\frac{1}{2}\left(k-k_{c,i}\right)^{2}$ (26)

and as a result, the SCGF $\lambda(k)$ shows a jump in the second derivative
at $k=k_{c,i}$. As we show below, this leads to a jump in the second
derivative of $I(\rho)$ which we interpret as a second order dynamical phase
transition (see next subsection for details for the particular case $N=2$).

For $k<k_{c,N-1}$, $\lambda(k)=\lambda_{1}(k)+\dots+\lambda_{M}(k)$, where $M$
is the number of single-body bound states. The eigenfunction
$r_{k}(X_{1},\dots,X_{N})$ in Eq. (14) is now given by the following Slater
determinant footnote:SlaterNormalization

$\Psi\left(X_{1},\dots,X_{N}\right)=\left|\begin{array}[]{ccc}\psi_{1}\left(X_{1}\right)&\dots&\psi_{N}\left(X_{1}\right)\\\
\vdots&\ddots\\\
\psi_{1}\left(X_{N}\right)&&\psi_{N}\left(X_{N}\right)\end{array}\right|$ (27)

where $\psi_{j}(X_{i})$ is the single particle state function. For
$k<k_{c,N-1}$, $M<N$ so the many-body energy spectrum is continuous, and $N-M$
of the fermions will not be localized around the interval $[-1,1]$. As we show
below, the interpretation of this in terms of the original problem is that $M$
of the Brownian particles remain in the vicinity of the interval $[-1,1]$ and
the remaining $N-M$ particles do not.

In the next subsection, we will analyze the results for a fixed $N$. We
consider $N=2$ non-crossing Brownian particles occupation fraction, as it
displays all of the interesting features for any general $N$.

###  IV.4 $N=2$

For $N=2$, the two lowest single-body energy levels are found from Eqs. (18)
and Eq. (19) which together take the form

$k=A_{2}\left(1+\cot^{2}\sqrt{2A_{2}}\right)=A_{1}\left(1+\tan^{2}\sqrt{2A_{1}}\right)$
(28)

There is a critical value of $k=k_{c}\equiv k_{c,1}=\frac{\pi^{2}}{8}$
(corresponds to $A_{1}=0.436191\dots$ and $A_{2}=\frac{\pi^{2}}{8}$), below
which only one physically possible solution $A_{1}$ contributes.

Whereas for $k>k_{c}$, the $\lambda_{2}(k)$ (corresponding to the solution
$A_{2}$) contributes and the total SCGF of the problem is

$\lambda(k)=\lambda_{1}(k)+\lambda_{2}(k).$ (29)

Fig. 2 shows the plot of the SCGF for the two non-crossing Brownian particle
occupation fraction (shown by the solid line). For $k<k_{c}$ the SCGF
coincides with $\lambda_{1}(k)$. At $k=k_{c}$, $\lambda_{2}(k)$ starts to
contribute and continues to exist with $\lambda_{1}(k)$ for all $k>k_{c}$.

Similarly, for the rate function, there is a critical value of the occupation
fraction $\rho=\rho_{c}=\lambda_{1}^{\prime}(k_{c})=0.84375\dots$, below which
the rate function coincides with the rate function $I_{1}(\rho)$ of the case
$N=1$. For $\rho>\rho_{c}$ the rate function is given by the Legendre
transform of the SCGF $\lambda_{1}(k)+\lambda_{2}(k)$. Fig. 3 shows the plot
of the rate function for $N=2$.

Next, we study the asymptotic behaviors of the SCGF, with particular emphasis
on the behavior near the critical point $k\longrightarrow k_{c}$
($\rho\longrightarrow\rho_{c}$) where we study the dynamical phase transition
in detail.

In the $k\longrightarrow 0$ limit which corresponds to $\rho\longrightarrow
0$, the SCGF and the rate function coincide with their counterparts for $N=1$,
and their asymptotic behaviors in these limits are given by the first lines of
Eqs. (24) and (25), respectively.

Figure 3:  (Color online) Solid line: The plot of the rate function $I(\rho)$
as a function of the occupation fraction $\rho$ for $N=2$ non-crossing
Brownian particles. The solid circle represents the critical point
$\rho_{c}=0.84375\dots$. For $\rho<\rho_{c}$ the rate function coincides with
the rate function for the corresponding single particle problem ($N=1$).The
dot-dashed line represents the continuation of the single particle rate
function for $\rho>\rho_{c}$. The red dashed line depicts the asymptotic
behavior of $I(\rho)$ for $\rho\longrightarrow 2$. The second derivative of
$I(\rho)$ has a jump discontinuity at $\rho=\rho_{c}$, which we interpret as a
second order dynamical phase transition. This can be seen as a corner
singularity in the inset which shows $dI/d\rho$. Units are chosen such that
$D=\frac{1}{2}$ and $l=1$.

In the critical regime, the SCGF behaves as

$\lambda\left(k\right)\simeq\begin{cases}\lambda_{1}\left(k\right),&k<k_{c},\\\\[5.69054pt]
\lambda_{1}\left(k\right)+\frac{1}{2}\left(k-k_{c}\right)^{2},&k-k_{c}\ll
k_{c}.\end{cases}$ (30)

Similarly the rate function near $\rho_{c}$ behaves as

$\\!I(\rho)\\!\simeq\\!\begin{cases}\\!\\!I_{1}(\rho),&\rho<\rho_{c},\\\\[5.69054pt]
\\!\\!I_{1}(\rho_{c})+k_{c}(\rho-\rho_{c})+B\\!\left(\rho-\rho_{c}\right)^{2}\\!,&\rho-\rho_{c}\\!\ll\\!\rho_{c}.\end{cases}$
(31)

where $B=\frac{1}{2\left[\lambda_{1}^{\prime\prime}(k_{c})+1\right]}\simeq
0.447705\dots$. The detailed calculations are presented in Appendix B.

Clearly the second derivative of the SCGF jumps at the critical point. As we
now show, the same is true for the rate function. The first derivative of the
rate function $I(\rho)$ is continuous because $I_{1}^{\prime}(\rho_{c})=k_{c}$
but $I^{\prime}(\rho)$ shows a corner singularity as shown in the bottom inset
of Fig. 3. Indeed, the second derivative shows a jump discontinuity (see
Appendix B)

$I^{\prime\prime}(\rho)\simeq\begin{cases}I_{1}^{\prime\prime}(\rho)=\frac{1}{\lambda_{1}^{\prime\prime}(k)}&\rho<\rho_{c}\\\
2B&\rho-\rho_{c}\ll\rho_{c}\end{cases}$ (32)

We interpret the jump in the 2nd derivative of $I$ as a dynamical phase
transition of second order baek ; baek1 .

In the $k\longrightarrow\infty$ limit the SCGF behaves as

$\lambda(k)\simeq
2k-\frac{5\pi^{2}}{8}+\frac{5\pi^{2}}{4}\frac{1}{\sqrt{2k}}\,.$ (33)

Taking the Legendre transform of this expression, we find that the
corresponding limit $\rho\longrightarrow 2$ of the rate function behaves as

$I(\rho)\simeq\frac{5\pi^{2}}{8}-3\left(\frac{5\pi^{2}}{8\sqrt{2}}\right)^{2/3}\left(2-\rho\right)^{1/3}$
(34)

The details of the calculations are given in Appendix B. The asymptotic
behaviors of the SCGF and the rate function are shown in Figs. 2 and 3
respectively by the dashed red lines.

###  IV.5 Large $N$ behavior

Let us consider the large-$N$ limit, where, as we show below, a universal
behavior emerges at $\rho\gg 1$. This corresponds to $k\longrightarrow\infty$,
where the single-body energy levels can be approximated as the energy levels
of particle-in-a-box i.e. $E_{i}\simeq-k+\frac{\pi^{2}}{8}i^{2}$.
$-\lambda(k)$ for large $N$ limit is thus the sum of the energy levels up to
$M$, where $M$ is the energy level at which $E_{M}\simeq 0$, which gives

$M\simeq\frac{2\sqrt{2k}}{\pi}$ (35)

Thus the SCGF can be written as

$\displaystyle\lambda(k)$ $\displaystyle=$
$\displaystyle-(E_{1}+E_{2}+....+E_{M})$ (36) $\displaystyle\simeq$
$\displaystyle
Mk-\frac{\pi^{2}}{24}M^{3}\simeq\frac{4\sqrt{2}}{3\pi}\,\,k^{3/2}\,.$

Taking the Legendre transform, we obtain

$I(\rho)\simeq\frac{\pi^{2}}{24}\,\,\rho^{3}$ (37)

Eq. (37) is valid in the limit $1\ll\rho<N$. Despite the fact that $I(\rho)$
has $N-1$ singularities, we find that the limiting behavior (37) is
nevertheless smooth. Fig. 4 shows the exact SCGF and rate function together
with the limiting behaviors for $N=5$, which turns out to be large enough to
observe excellent agreement. There are four critical points at which the
dynamical phase transitions of second order occur.

(a)

(b)

Figure 4:  The plot of the (a) SCGF as a function of $k$ and (b) rate function
as a function of $\rho$ for $N=5$ Brownian particles. There are four critical
points where the second-order dynamical phase transitions happen, shown by
black dots. The red dashed lines depict the asymptotic large-$k$ and
large-$\rho$ behaviours, see Eqs. (36) and (37). The exact and asymptotic
expressions agree surprisingly well blue and are indistinguishable at the
scale of the figure.

###  IV.6 Conditioned process

$|\Psi(X_{1},\dots,X_{N})|^{2}$ is the joint distribution of the positions of
the particles at some intermediate time $0<t<T$ conditioned on a given value
of $\rho$ Touchetteoneparticle . Therefore, the $N$ phases correspond to
conditioned processes in which $M$ particles stay in the vicinity of the
interval $[-l,l]$ (where $M=1,\dots,N$). It is then intriguing to study the
properties of the conditioned process in more detail. This can be done by
exploiting the interpretation of the joint distribution as that of the
positions of $M$ noninteracting fermions confined by the potential $V_{k}(x)$.

Quantum fluctuations of interacting or noninteracting fermions, possibly
confined by an external potential, have attracted much recent attention castin
; CMV2012 ; eisler_prl ; marino_prl ; DLMS16PRA , in particular due to their
importance in applications such as experiments in cold atoms Fermicro1 ;
Fermicro2 ; Fermicro3 ; Pauli ; BDZ08 ; flattrap . Using a variety of
analytical methods, such as determinantal processes and surprising connections
to random matrix theory eisler_prl ; marino_prl ; DLMS16PRA ; Calabrese_RMT ;
DLMS2019review ; DLMSS21 ; SLMS21 , fundamental properties such as the spatial
density of fermions and its correlations have been studied. For $M$
noninteracting fermions in the presence of a trapping potential at zero
temperature (i.e., in the many-body ground state), the quantum correlations
are encoded in a fundamental object called the kernel, that is given by

$K_{M}\left(x,y\right)=\sum_{j=1}^{M}\psi_{j}^{*}\left(x\right)\psi_{j}\left(y\right),$
(38)

where $\psi_{1}(x),\dots,\psi_{M}(x)$ are the wave functions of the $M$
single-body lowest energy states. The quantum correlations of the fermions’
positions can be expressed using the kernel. For instance, the density of
particles at point $x$ is given by
$\mathcal{R}_{M}\left(x\right)=K_{M}\left(x,x\right)$.

The above expressions for $\mathcal{R}_{M}$ and $K_{M}$ are valid for
arbitrary $M$, and (for the potential $V(x)=V_{k}(x)$) they describe the
density and density-density correlations of the $M$ particles that remain near
the interval $[-1,1]$, conditioned on a given value of the occupation
fraction. However, in the large-$M$ limit they approach universal behaviors as
we now describe. Over large scales, the density is in general given by the
local density (or Thomas-Fermi) approximation,
$\mathcal{R}_{M}\left(x\right)\simeq k_{F}\left(x\right)/\pi$ where
$k_{F}\left(x\right)\simeq\sqrt{2\left(\mu-V\left(x\right)\right)}$ is the
local Fermi wave vector and $\mu$ is the Fermi energy (we consider
$\hbar=m=1$). The density is normalized such that
$\int_{-\infty}^{\infty}\mathcal{R}_{M}\left(x\right)dx=M$. For the square-
well potential $V_{k}(x)$, as we saw above, $M$ is determined by the condition
$\mu=E_{M}\simeq 0$, so we obtain

$\mathcal{R}_{M}\left(x\right)\simeq\begin{cases}\sqrt{2k}/\pi\simeq
M/2,&x\in\left[-1,1\right],\\\\[5.69054pt]
0,&x\notin\left[-1,1\right],\end{cases}$ (39)

which thus determines $M\simeq 2\sqrt{2k}/\pi$ in agreement with Eq. (35)
above. This gives the mean inter-particle distance between particles inside
the interval, $d=1/\mathcal{R}_{M}\left(x\right)=2/M$. Moreover, the density
correlations inside the interval are given in terms of the celebrated sine
kernel eisler_prl ; DLMS16PRA ; DLMS2019review

$K_{M}\left(x,y\right)\simeq\frac{\sin\left[k_{F}\left(x\right)\left(x-y\right)\right]}{\pi\left(x-y\right)}.$
(40)

These expressions for the density and the kernel break down near the edges of
the interval, $x=\pm l$. In fact, the behavior of a noninteracting Fermi gas
in the presence of a step-function potential was analyzed in Ref. DLMSS21  in
some detail. A universal form of the density and kernel were obtained. For
$\mu\simeq 0$ (the so-called critical case), they are given by

$\mathcal{R}_{M}\left(x\right)\simeq\frac{1}{d}n_{c}\left(\frac{x}{d}\right),\quad
K_{M}\left(x,y\right)\simeq\frac{1}{d}\kappa_{c}\left(\frac{x}{d},\frac{y}{d}\right),$
(41)

respectively, where the universal functions $n_{c}$ and $\kappa_{c}$ are given
explicitly in DLMSS21 . They provide an interpolation of $\mathcal{R}_{M}$ and
$K_{M}$ between the regime inside the interval ($x\in[-l,l]$) and outside it
($x\notin[-l,l]$).

##  V Discussion

We have studied the fluctuations of the occupation fraction of $N$ non-
crossing Brownian particles in a finite interval. We found that the rate
function that describes these fluctuations at long times undergoes a sequence
of $N-1$ DPTs of second order. The $N$ different phases correspond to
different numbers of particles in the vicinity of the interval. We achieved
this by extending the DV formalism to noncrossing Brownian particles and found
that this maps the problem to that of $N$ noninteracting spinless fermions in
the presence of an effective potential. Under this mapping, the $N$ phases
correspond to different numbers of bound single-body states of the effective
potential.

This DPT observed here for any $N\geq 2$ is of very different nature to the
DPT that has been previously found for a single particle ($N=1$). The latter
DPT only occurs if the particle experiences an additional external drift
(which leads to breaking of time-reversal symmetry) TouchetteMinimalModel ;
Touchetteoneparticle , and its properties are quite different to those of the
DPTs that we found here. For instance, it involves coexistence between two
different phases, which correspond to the particle being confined (or not) to
the vicinity of the interval. However, in the present work, we have found
$N-1$ DPTs that separate between $N$ temporally-homogeneous phases, each
corresponding to a different number of particles which are in the vicinity of
the given interval for the entire dynamics.

It would be interesting to study fluctuations of the occupation time in
presence of other interactions between the Brownian particles (instead of or
in addition to the non-crossing condition). Another natural extension of the
present problem would be to study the effect of a drift in the system, or to
study the occupation time for other processes such as fractional Brownian
motion or for active particles. Moreover, one could study more detailed
quantities such as the empirical distribution of the positions of the
particles.

### Acknowledgments

We acknowledge several useful discussions with Pierre Le Doussal on related
topics.

##  Appendix A Donsker-Varadhan formalism for noncrossing particles

We first consider the unnormalized joint probability distribution
$Q(X_{1},\dots,X_{N},\mathcal{R},t)$ of the $X_{i}$’s and
$\mathcal{R}=t\rho_{t}$. $Q$ satisfies the Fokker-Planck equation

$\frac{\partial Q}{\partial
t}=L^{\dagger}\,Q-\mathds{1}_{\\{-l,l\\}}(x)\,\frac{\partial
Q}{\partial\mathcal{R}}$ (42)

together with the non-crossing boundary conditions given in Eq. (10) at
$X_{i}=X_{i+1}$. It is convenient to define the generating function

$G_{k}(X_{1},\dots,X_{N},t)=\int
e^{k\,\mathcal{R}}Q(X_{1},\dots,X_{N},\mathcal{R},t)d\mathcal{R}.$ (43)

According to Feynman-Kac formula kac ; Majumdar2007 , $G_{k}$ evolves with
time as

$\frac{\partial G_{k}(X_{1},\dots,X_{N},t)}{\partial
t}=\mathcal{L}_{k}^{\dagger}\,G_{k}(X_{1},\dots,X_{N},t)$ (44)

where

$\mathcal{L}_{k}^{\dagger}=L^{\dagger}+k\mathds{1}_{\\{-l,l\\}}(x)$ (45)

known as the tilted generator. Note that $G_{k}$ also satisfies the boundary
conditions (10).

Eq. (44) can be solved by expanding $G_{k}$ over the eigenbasis of
$\mathcal{L}_{k}$ footnote:SumIntegral

$G_{k}(X_{1},\dots,X_{N},t)=\sum_{i}\,c_{i}\,e^{t\zeta_{i}}\,r_{k}^{i}(X_{1},\dots,X_{N})$
(46)

where $\zeta_{i}$ and $r_{k}^{i}(X_{1},\dots,X_{N})$ are the eigenvalues and
eigenfunctions of $\mathcal{L}_{k}$ respectively, and the $c_{i}$’s are
coefficients that depend on the initial condition.

For $T\longrightarrow\infty$, the sum in Eq. (46) is dominated by the largest
eigenvalue $\zeta_{\max}$ (and if $\zeta_{\max}$ doesn’t exist then instead of
it we must use the supremum of the $\zeta_{i}$’s). In this long time limit the
nomalization factor $\mathcal{N}(t)$ goes as a power law $\mathcal{N}(t)\sim
t^{-N(N-1)/4}$ fishersurvival ; fishersurvival1 ; viciouswalker ;
viciouswalker1 ; GLMS19 . Comparing Eq. (46) to the definition of the SCGF in
Eq. (6) it can be concluded that $\lambda(k)=\zeta_{\max}$, as the logarithm
of $\mathcal{N}(t)$ can be neglected.

##  Appendix B Asymptotic behaviors of $\lambda(k)$ and of $I(\rho)$

In this section, we show the details of the calculation of some of the
asymptotic behaviors of the SCGF and the rate function. We mostly focus on the
case $N=2$, and consider the limits $k\longrightarrow 0$ ($\rho\longrightarrow
0$), $k\rightarrow\infty$ ($\rho\longrightarrow 2$) and near the critical
point, $k\simeq k_{c}$. We use Eqs. (18)-(20) to study the asymtotic
behaviours for different limits of $A$ (or $k$).

###  B.1 $k\simeq k_{c}$

For any general $N$, for $k$, that is near the critical point
$k_{c,i}=\frac{i^{2}\pi^{2}}{8}$ but slightly larger (i.e., for $k>k_{c,i}$
with $k-k_{c,i}\ll k_{c,i}$), using $\tan x=\tan\left(x-n\pi\right)\simeq
x-n\pi$ (for even $i$) and $\cot x=\cot\left(x-n\pi\right)\simeq x-n\pi$ (for
odd $i$) for $x\simeq n\pi$, the $k$ and $\lambda_{i+1}(k)$ behaves for
$A_{i}\longrightarrow\frac{i^{2}\pi^{2}}{8}$ as Eq. (18) becomes

$k\simeq\frac{1}{2}\left(\frac{i^{2}\pi^{2}}{8}-A_{i}\right)^{2}+A_{i}$ (47)

leading to

$\displaystyle\lambda_{i+1}(k)$ $\displaystyle=$ $\displaystyle
k-A_{i}\simeq\frac{1}{2}\left(\frac{i^{2}\pi^{2}}{8}-A_{i}\right)^{2}$ (48)
$\displaystyle\simeq$
$\displaystyle\frac{1}{2}\left(k-\frac{i^{2}\pi^{2}}{8}\right)^{2}$

which coincides with Eq. (26) of the main text. Using the above two equations
the behaviour of the SCGF $\lambda(k)$ for $N=2$ can now be written as

$\lambda(k)\simeq\lambda_{1}(k)+\frac{1}{2}\left(k-k_{c}\right)^{2}$ (49)

which is Eq. (30) in Sec. IV. Taking the derivative we have:

$\rho=\frac{d\lambda}{dk}\simeq\lambda_{1}^{\prime}(k)+(k-k_{c})\simeq\rho_{c}+\left[\lambda_{1}^{\prime\prime}(k_{c})+1\right](k-k_{c})$
(50)

where we used $\rho_{c}=\lambda_{1}^{\prime}(k_{c})$. Inverting this relation,
we obtain

$k\simeq k_{c}+\frac{\rho-\rho_{c}}{\lambda_{1}^{\prime\prime}(k_{c})+1}.$
(51)

Taking the Legendre transform of the Eq. (30) while using Eqs. (50) and (51),
the expression of the rate function at $\rho$ that is slightly larger than
$\rho_{c}$ can now be expressed as

$I(\rho)\simeq
I_{1}(\rho_{c})+k_{c}(\rho-\rho_{c})+\frac{1}{2\left[\lambda_{1}^{\prime\prime}(k_{c})+1\right]}\left(\rho-\rho_{c}\right)^{2}$
(52)

which coincides with the second line of Eq. (31) in Sec. IV.

###  B.2 $k\to 0$

In the $k\rightarrow 0$ limit, the SCGF $\lambda(k)=\lambda_{1}(k)\approx 0$,
which implies $A\longrightarrow 0$. Around $A=0$, using Eqs. (18) and (20),
the asymptotic behaviors of $\lambda(k)$ and $k$ are

$\lambda(k)\simeq 2A^{2},\qquad k\simeq A+2A^{2}$ (53)

respectively. This implies

$\lambda(k)\simeq 2k^{2},$ (54)

coinciding with the first line of Eq. (24) in Sec. IV. The Legendre transform
of $\lambda(k)$ in this limit is thus

$I(\rho)\simeq\frac{\rho^{2}}{8}$ (55)

which coincides the first line of Eq. (25) in Sec. IV.

###  B.3 $k\to\infty$

In the limit $k\longrightarrow\infty$, $A_{1}$ goes as
$A_{1}\rightarrow\frac{\pi^{2}}{8}$. Considering Eq. (18) the asymptotic
behavior of the left hand side of Eq. (18) near this value of $A_{1}$ we
obtain

$\displaystyle k$ $\displaystyle\simeq$
$\displaystyle\frac{\pi^{4}}{32\left(\frac{\pi^{2}}{8}-A_{1}\right)^{2}}-\frac{3\pi^{2}}{8\left(\frac{\pi^{2}}{8}-A_{1}\right)}+\frac{1}{24}\left(21+\pi^{2}\right)$
(56) $\displaystyle-$
$\displaystyle\left(\frac{1}{3}-\frac{1}{2\pi^{2}}\right)\left(\frac{\pi^{2}}{8}-A_{1}\right)$

which yields

$\displaystyle\lambda_{1}(k)$ $\displaystyle=$ $\displaystyle k-A_{1}$ (57)
$\displaystyle\simeq$
$\displaystyle\frac{\pi^{4}}{32\left(\frac{\pi^{2}}{8}-A_{1}\right)^{2}}-\frac{3\pi^{2}}{8\left(\frac{\pi^{2}}{8}-A_{1}\right)}-\frac{\pi^{2}}{12}+\frac{7}{8}$
$\displaystyle+$
$\displaystyle\left(\frac{2}{3}+\frac{1}{2\pi^{2}}\right)\left(\frac{\pi^{2}}{8}-A_{1}\right).$

From Eqs. (56) and (57), we now express $\lambda_{1}(k)$ as a function of $k$
in the limit $k\gg 1$

$\lambda_{1}(k)\simeq k-\frac{\pi^{2}}{8}+\frac{\pi^{2}}{4\sqrt{2k}},$ (58)

which coincides with the second line of Eq. (24) in the main text. This result
gives the asymptotic behavior of the SCGF for the case $N=1$.

For $N=2$, we must add the contribution of $\lambda_{2}(k)$, whose asymptotic
behavior in the large-$k$ limit we now calculate. For the second energy level,
$A_{2}$ goes as $A_{2}\rightarrow\frac{\pi^{2}}{2}$ and by using the
corresponding behavior of Eq. (19) and following similar steps as above, we
find

$\displaystyle\lambda_{2}(k)\simeq
k-\frac{\pi^{2}}{2}+\frac{\pi^{2}}{\sqrt{2k}}\,.$ (59)

Now using Eqs. (58) and (59), we find the asymptotic behavior of the SCGF
$\lambda(k)=\lambda_{1}(k)+\lambda_{2}(k)$ at $k\longrightarrow\infty$

$\lambda(k)\simeq 2k-\frac{5\pi^{2}}{8}+\frac{5\pi^{2}}{4}\frac{1}{\sqrt{2k}}$
(60)

which coincides with Eq. (33) in Sec. IV. We now calculate the Legendre
transform of this expression, in order to find the corresponding asymptotic
behavior of the rate function. Taking the derivative of Eq. (60) w.r.t $k$, we
get the limiting behaviour of $\rho\longrightarrow 2$ as

$\rho=\frac{d\lambda(k)}{dk}=2-\frac{5\pi^{2}}{8\sqrt{2}k^{3/2}}.$ (61)

By combining Eqs. (60) and (61), we find the behavior of the rate function in
the limit $\rho\longrightarrow 2$,

$I(\rho)=k\rho-\lambda(k)\simeq\frac{5\pi^{2}}{8}-3\left(\frac{5\pi^{2}}{8\sqrt{2}}\right)^{2/3}\left(2-\rho\right)^{1/3},$
(62)

which coincides with Eq. (34) in Sec. IV. The second line of Eq. (25) (for the
case $N=1$) is obtained similarly, by calculating the Legendre transform of
Eq. (58).

## References

  * (1) T. Bodineau, B. Derrida, Current Fluctuations in Nonequilibrium Diffusive Systems: An Additivity Principle, Phys. Rev. Lett. 92, 180601 (2004).
  * (2) R. J. Harris, A. Rákos, G. M. Schütz Current fluctuations in the zero-range process with open boundaries, J. Stat. Mech. 2005, P08003 (2005).
  * (3) B. Derrida, Non-equilibrium steady states: fluctuations and large deviations of the density and of the current, J. Stat. Mech. P07023 (2007).
  * (4) O. Ovaskainen and B. Meerson, Stochastic models of population extinction, Trends Ecol. Evol. 25, 643 (2010).
  * (5) N. Merhav, Y. Kafri, Bose–Einstein condensation in large deviations with applications to information systems, J. Stat. Mech. 2010, P02011 (2010).
  * (6) O. Cohen, D. Mukamel, Phase Diagram and Density Large Deviations of a Nonconserving ABC Model, Phys. Rev. Lett. 108, 060602 (2012)
  * (7) L. Bertini, A. De Sole, D. Gabrielli, G. Jona-Lasinio, C. Landim, Macroscopic fluctuation theory, Rev. Mod. Phys. 87, 593 (2015).
  * (8) M. Assaf and B. Meerson, WKB theory of large deviations in stochastic populations, J. Phys. A: Math. Theor. 50, 263001 (2017).
  * (9) S. S. Varadhan, Large Deviations and Applications, CBMS-NSF Regional Conference Series in Applied Mathematics, No. 46 (SIAM, Philadelphia, 1984).
  * (10) R. S. Ellis, Entropy, Large Deviations, and Statistical Mechanics, (Springer, New York) 1985.
  * (11) Y. Oono, Large Deviation and Statistical Physics, Prog. Theor. Phys. Suppl. 99, 165 (1989).
  * (12) A. Dembo and O. Zeitouni, Large Deviations Techniques and Applications, 2nd ed. (Springer, New York, 1998).
  * (13) F. den Hollander, Large Deviations, Fields Institute Monographs, vol. 14 (AMS, Providence, Rhode Island, 2000).
  * (14) H. Touchette, The large deviation approach to statistical mechanics, Phys. Rep. 478, 1 (2009).
  * (15) R. L. Jack, Ergodicity and large deviations in physical systems with stochastic dynamics, Eur. Phys. J. B 93, 1 (2020).
  * (16) T. Bodineau, B. Derrida, Distribution of current in nonequilibrium diffusive systems and phase transitions, Phys. Rev. E 72, 066110 (2005).
  * (17) J. P. Garrahan, R. L. Jack, V. Lecomte, E. Pitard, K. van Duijvendijk, F. van Wijland, Dynamical First-Order Phase Transition in Kinetically Constrained Models of Glasses, Phys. Rev. Lett. 98, 195702 (2007).
  * (18) G. Bunin, Y. Kafri, D. Podolosky, Non-differentiable large-deviation functionals in boundary-driven diffusive systems, J. Stat. Mech. 2012, L10001, (2012).
  * (19) C. P. Espigares, P. L. Garrido, P. I. Hurtado, Dynamical phase transition for current statistics in a simple driven diffusive system, Phys. Rev. E 87, 032115 (2013).
  * (20) Y. Baek, Y. Kafri, Singularities in large deviation functions, J. Stat. Mech. 2015, P08026 (2015).
  * (21) Y. Baek, Y. Kafri, V. Lecomte, Dynamical Symmetry Breaking and Phase Transitions in Driven Diffusive Systems,  Phys. Rev. Lett. 118, 030604 (2017) .
  * (22) Y. Baek, Y. Kafri, V. Lecomte, Dynamical phase transitions in the current distribution of driven diffusive channels,  J. Phys. A: Math. Theor. 51 105001 (2018) .
  * (23) I. Dornic and C. Godréche, Large deviations and nontrivial exponents in coarsening systems, J. Phys. A: Math. Gen. 31, 5413 (1998).
  * (24) J.-M. Drouffe and C. Godréche, Stationary definition of persistence for finite-temperature phase ordering, J. Phys. A: Math. Gen. 31, 9801 (1998).
  * (25) A. Baldassarri, J. P. Bouchaud, I. Dornic, C. GodrécheStatistics of persistent events: An exactly soluble model, Phys. Rev. E, 59, 1, 1999.
  * (26) J.-P. Bouchaud, M. Potters Theory of Financial Risk and Derivative Pricing: From Statistical Physics to Risk Management, 2nd ed. (Cambridge University Press, Cambridge, 2000).
  * (27) A. Montanari, R. Zecchina, Optimizing searches via rare events, Phys. Rev. Lett. 88, 178701 (2002).
  * (28) V. Kishore, M. S. Santhanam, R. E. Amritkar, Extreme events and event size fluctuations in biased random walks on networks, Phys. Rev. E 85, 056120 (2012).
  * (29) P. Lévy, Sur certains processus stochastiques homogènes. , Compositio Mathematica 7 (1940): 283-339.
  * (30) M. Kac, On Distributions of Certain Wiener Functionals, Trans. Am. Math. Soc. 65, 1, 1949.
  * (31) D. A. Darling, M. Kac, On occupation times for Markoff processes, Trans. Amer. Math. Soc. 84, 444-458, 1957.
  * (32) C. Godréche, J.M. Luck, Statistics of the Occupation Time of Renewal Processes, Journal of Statistical Physics 104, 489–524 (2001).
  * (33) S. N. Majumdar, A. Comtet, Local and Occupation Time of a Particle Diffusing in a Random Medium, Phys. Rev. Lett. 89, 6 2002.
  * (34) P. Tsobgni Nyawo, H. Touchette, A minimal model of dynamical phase transition, Europhys. Lett. 116, 50009 (2016).
  * (35) P. Tsobgni Nyawo, H. Touchette, Dynamical phase transition in drifted Brownian motion, Phys. Rev. E 98, 052103 (2018).
  * (36) T. Agranov, P. L. Krapivsky, and B. Meerson, Occupation-time statistics of a gas of interacting diffusing particles, Phys. Rev. E 99, 052102 (2019)
  * (37) M. Kimura and T. Akimoto, Occupation time statistics of the fractional Brownian motion in a finite domain, Phys. Rev. E 106, 064132 (2022).
  * (38) Z. Toroczkai, T.J. Newman, S. Das Sarma, Sign-time distributions for interface growth, Phys. Rev. E 60, R1115 (1999).
  * (39) G. Messin, J.P. Hermier, E. Giacobino, P. Desbiolles, M. Dahan, Bunching and antibunching in the fluorescence of semiconductor nanocrystals  Opt. Lett. 26, 23, (2001).
  * (40) G.H. Weiss and P.P. Calabrese, Occupation times of a CTRW on a lattice with anomalous sites, Physica A, 234, 443, (1996).
  * (41) M. D. Donsker, S.R.S. Varadhan, Assymptotics of Wiener Sausage, Comm. Pure Appl. Math. 28, 1 (1975).
  * (42) M. D. Donsker, S.R.S. Varadhan, On the principal eigenvalue of second-order elliptic differential operators, Comm. Pure Appl. Math. 29, 389 (1976).
  * (43) S. N. Majumdar and A. J. Bray, Large-deviation functions for nonlinear functionals of a Gaussian stationary Markov process, Phys. Rev. E 65, 051112, 2002.
  * (44) S. N. Majumdar, Brownian Functionals in Physics and Computer Science, Current Science 89, 2076 (2005).
  * (45) N. R. Smith, Large deviations in chaotic systems: Exact results and dynamical phase transition, Phys. Rev. E 106, L042202, (2022).
  * (46) R. Graham, T. Tél, Existence of a Potential for Dissipative Dynamical Systems, Phys. Rev. Lett. 52(2), 9–12 (1984).
  * (47) P. -G. de Gennes, Soluble Model for Fibrous Structures with Steric Constraints, J. Chem. Phys. 48, 2257 (1968) .
  * (48) M.E. Fisher,  Walks, walls, wetting, and melting J. Stat. Phys. 34, 667–729 (1984).
  * (49) J. W. Essam, A. J. Guttmann, Vicious walkers and directed polymer networks in general dimensions, Phys. Rev. E 52, 5849, (1995).
  * (50) A.J. Bray, K. Winkler, Vicious walkers in a potential J. Phys. A Gen. Phys. 37, 2 (2004).
  * (51) G. Schehr, S. N. Majumdar, A. Comtet, J. Randon-Furling, Exact Distribution of the Maximal Height of p Vicious Walkers Phys. Rev. Lett. 101, 150601, (2008).
  * (52) G. Schehr, Extremes of N Vicious Walkers for Large N: Application to the Directed Polymer and KPZ Interfaces, J Stat Phys 149, 385–410 (2012).
  * (53) G. Schehr, S. N. Majumdar, A. Comtet, P. J. Forrester, Reunion Probability of N Vicious Walkers: Typical and Large Fluctuations for Large N, J Stat Phys 150, 491–530 (2013).
  * (54) G. B. Nguyen, D. Remenik, Extreme statistics of non-intersecting Brownian paths, Electron. J. Probab. 22, 1 - 40, (2017).
  * (55) T. Gautié, P. Le Doussal, S. N. Majumdar, G. Schehr, Non-crossing Brownian paths and Dyson Brownian motion under a moving boundaryJ. Stat. Phys. 177,752 (2019).
  * (56) T. Gautié, N. R. Smith, Constrained non-crossing Brownian motions, fermions and the Ferrari–Spohn distribution, J. Stat. Mech. 2021, 033212, (2021).
  * (57) F. J. Dyson, A Brownian‐motion model for the eigenvalues of a random matrix, J. Math. Phys. 3, 1191 (1962).
  * (58) K. Johansson, From Gumbel to Tracy-Widom, Probab. Theory Relat. Fields 138, 75–112 (2007).
  * (59) M. Katori, Determinantal Process Starting from an Orthogonal Symmetry is a Pfaffian Process, J Stat Phys 146, 249–263 (2012).
  * (60) J. Gärtner, On Large Deviations from the Invariant Measure, Th. Prob. Appl. 22, 24 (1977).
  * (61) R. S. Ellis, Large Deviations for a General Class of Random Vectors,  Ann. Prob. 12, 1 (1984).
  * (62) H. Touchette, Legendre-Fenchel transforms in a nutshell, https://www.ise.ncsu.edu/fuzzy-neural/wp-content/uploads/sites/9/2019/01/or706-LF-transform-1.pdf (2005).
  * (63) D. A. Huse and M. E. Fisher, Commensurate melting, domain walls, and dislocations Phys. Rev. B 29, 239.
  * (64) C. Krattenthaler, A.J. Guttmann, X.G. Viennot, Vicious walkers, friendly walkers and Young tableaux: II. With a wall J. Phys. A Math. Gen. 33(48), 8835 (2000).
  * (65) S. N. Majumdar, Brownian Functionals in Physics and Computer Science, Current Science 89, 2076 (2005).
  * (66) H. Touchette, Introduction to dynamical large deviations of Markov processes, Physica A 504, 5 (2018).
  * (67) D. J. Griffiths, Introduction to Quantum Mechanics, 2nd Ed., 2005.
  * (68) B. C. Hall, Quantum theory for mathematicians, Springer, (2013).
  * (69) The absence of the usual normalization factor $1/\sqrt{N!}$ before the Slater determinant in Eq. (27) is because it is constrained to $X_{1}<\dots<X_{N}$ in contrast to the usual case with fermions, thus leading to a different normalization factor.
  * (70) Y. Castin, in Quantum gases in low dimensions, J. Phys. IV France 116, 89 (2004); arXiv:0407118.
  * (71) P. Calabrese, M. Mintchev, E. Vicari, Exact relations between particle fluctuations and entanglement in Fermi gases, EPL 98, 20003 (2012).
  * (72) V. Eisler, Universality in the full counting statistics of trapped fermions, Phys. Rev. Lett. 111, 080402 (2013).
  * (73) R. Marino, S. N. Majumdar, G. Schehr, P. Vivo, Phase transitions and edge scaling of number variance in Gaussian random matrices, Phys. Rev. Lett. 112, 254101 (2014).
  * (74) D. S. Dean, P. Le Doussal, S. N. Majumdar, G. Schehr, Noninteracting fermions at finite temperature in a d-dimensional trap: universal correlations, Phys. Rev. A 94, 063622 (2016).
  * (75) I. Bloch, J. Dalibard, and W. Zwerger, Many-body physics with ultracold gases, Rev. Mod. Phys. 80, 885 (2008).
  * (76) L. W. Cheuk, M. A. Nichols, M. Okan, T. Gersdorf, R.Vinay, W. Bakr, T. Lompe, and M. Zwierlein, Quantum-gas microscope for fermionic atoms, Phys. Rev. Lett. 114, 193001 (2015).
  * (77) E. Haller, J. Hudson, A. Kelly, D. A. Cotta, B. Peaudecerf, G. D. Bruce, and S. Kuhr, Single-atom imaging of fermions in a quantum-gas microscope, Nat. Phys. 11, 738 (2015).
  * (78) M. F. Parsons, F. Huber, A. Mazurenko, C. S. Chiu, W. Setiawan, K. Wooley-Brown, S. Blatt, and M. Greiner, Site-resolved imaging of fermionic 6Li in an optical lattice, Phys. Rev. Lett. 114, 213002 (2015).
  * (79) B. Mukherjee, Z. Yan, P. B. Patel, Z. Hadzibabic, T.Yefsah, J. Struck, and M. W. Zwierlein, Homogeneous atomic Fermi gases, Phys. Rev. Lett. 118, 123401 (2017).
  * (80) M. Holten, L. Bayha, K. Subramanian, C. Heintze, P. M. Preiss, and S. Jochim, Observation of Pauli crystals, Phys. Rev. Lett. 126, 020401 (2021).
  * (81) P. Calabrese, P. Le Doussal, S. N. Majumdar, Random matrices and entanglement entropy of trapped Fermi gases, Phys. Rev. A 91, 012303 (2015).
  * (82) D. S. Dean, P. Le Doussal, S. N. Majumdar, G. Schehr, Noninteracting fermions in a trap and random matrix theory, J. Phys. A: Math. Theor. 52, 144006 (2019).
  * (83) D. S. Dean, P. Le Doussal, S. N. Majumdar, G. Schehr, and N. R. Smith, Kernels for noninteracting fermions via a Green’s function approach with applications to step potentials, J. Phys. A: Math. Theor. 54, 084001 (2021).
  * (84) N. R. Smith, P. Le Doussal, S. N. Majumdar, and G. Schehr, Full counting statistics for interacting trapped fermions, SciPost Phys. 11, 110 (2021).
  * (85) The $\sum$ in Eq. (46) is to be understood as a sum over the discrete part of the spectrum of the tilted operator plus an integral over the continuous part of the spectrum.