# Triamese-ViT: A 3D-Aware Method for Robust Brain Age Estimation from MRIs

Zhaonian Zhang
Lancaster University
Lancaster LA1 4YZ, England
<EMAIL_ADDRESS>Richard Jiang
Lancaster University
Lancaster LA1 4YZ, England
<EMAIL_ADDRESS>
###### Abstract

The integration of machine learning in medicine has significantly improved
diagnostic precision, particularly in the interpretation of complex structures
like the human brain. Diagnosing challenging conditions such as Alzheimer’s
disease has prompted the development of brain age estimation techniques. These
methods often leverage three-dimensional Magnetic Resonance Imaging (MRI)
scans, with recent studies emphasizing the efficacy of 3D convolutional neural
networks (CNNs) like 3D ResNet. However, the untapped potential of Vision
Transformers (ViTs), known for their accuracy and interpretability, persists
in this domain due to limitations in their 3D versions. This paper introduces
Triamese-ViT, an innovative adaptation of the ViT model for brain age
estimation. Our model uniquely combines ViTs from three different orientations
to capture 3D information, significantly enhancing accuracy and
interpretability. Tested on a dataset of 1351 MRI scans, Triamese-ViT achieves
a Mean Absolute Error (MAE) of 3.84, a 0.9 Spearman correlation coefficient
with chronological age, and a -0.29 Spearman correlation coefficient between
the brain age gap (BAG) and chronological age, significantly better than
previous methods for brian age estimation. A key innovation of Triamese-ViT is
its capacity to generate a comprehensive 3D-like attention map, synthesized
from 2D attention maps of each orientation-specific ViT. This feature is
particularly beneficial for in-depth brain age analysis and disease diagnosis,
offering deeper insights into brain health and the mechanisms of age-related
neural changes.

##  1 Introduction

Figure 1: The structure of Triamese-ViT. We reshape MRI scans into three
distinct viewpoints, dividing each into fixed-size patches. These patches are
then linearly embedded, enhanced with position embeddings, and subsequently
inputted into a standard Transformer encoder. The encoder’s output is directed
through MLP Heads to generate three separate predictions. These predictions
are then integrated using the Triamese MLP, culminating in the final result.

Aging naturally impacts all body parts, including the brain, which is
particularly sensitive to age-related changes, heightening the risk of
diseases like Alzheimer’s. Traditional diagnostic methods for brain diseases
are often slow and heavily reliant on subjective clinical judgment [5, 35],
leading to potential delays or inaccuracies in diagnosis. Such shortcomings
can be critical, risking the loss of crucial treatment time and possibly
exacerbating the patient’s condition.

Recent advances in machine learning, particularly deep learning, have
revolutionized brain diagnostics [14, 21, 4, 54]. Deep learning’s ability to
estimate brain age from images is pivotal in detecting age-related diseases
[14, 7, 22]. The brain age gap (BAG), the disparity between estimated and
actual brain age, is a crucial indicator [14]. A younger-appearing brain
indicates health, while an older-looking brain may signal conditions like
Alzheimer’s [7], psychosis [13], mild cognitive impairment [22], or depression
[23]. Enhancing brain age estimation algorithms is essential, facilitating
early disease detection and offering patients hope through improved treatment
prospects.

Current brain age estimation largely depends on convolutional neural networks
(CNNs) trained on either 3D MRI scans [14] or 2D slices from these scans [26,
8]. While CNNs excel in image processing and can utilize full 3D MRI data for
comprehensive predictions, they struggle with global feature representation
due to their focus on small, local pixel groups [27]. This results in a loss
of critical details, particularly in the analysis of the brain’s complex
structures. Moreover, CNNs’ lack of transparency poses challenges in
Explainable AI, making their predictions difficult to interpret in medical
diagnostics, especially when pinpointing specific brain abnormalities [3].

The Vision Transformer (ViT) presents a notable advantage over traditional
CNNs in image detail analysis [20]. By segmenting images into patches and
transforming each via a convolutional layer into a high-dimensional space, ViT
effectively extracts intricate details and understands inter-part interactions
[31]. Its attention map feature provides insight into its focus and prediction
rationale. However, ViT has limitations: it may overlook broader, global
contexts [47] and is primarily suited for 2D images [1]. For brain age
prediction, since researchers often use flat 2D slices of MRI scans as input
for ViTs [24], although this approach can greatly reduce the time of training
and prediction, it doesn’t take into account the full 3D structure of the
brain, leading to a loss of important depth information.

To address the need for both detailed analysis and understanding of
interactions across different areas within an image, we introduce a new model,
Triamese-ViT, whose design is a first in the realm of brain age estimation.
This model is inspired by the Siamese Networks, which has been successfully
applied to image classification, detection and comparison [16, 52, 51]. The
novelty of our model is that it adeptly transforms 3D imaging challenges into
a 2D analytical framework, effectively preserving the rich information
contained in three-dimensional images, capturing the relationship between
different areas in the brain very well, and has the same fast process speed as
2D images at the same time.

A standout feature of our model is its capability to generate 3D-like
attention maps, offering profound interpretability. We meticulously compared
these maps with conventional Explainable AI (XAI) methods and corroborated
with findings from relevant medical research [44, 45, 43, 50, 48], affirming
their scientific validity, logical coherence, and precision. This synergy of
advanced imaging analysis and interpretability positions our model as a
significant tool in medical diagnostics and research, especially in areas
demanding high fidelity in three-dimensional medical data understanding.

We trained and tested our model on MRI scans from 1351 healthy individuals.
When compared with the state-of-the-art algorithms, Triamese-ViT not only
showed superior predictive accuracy but also demonstrated improvements in
fairness and interoperability. Our model achieves impressive results that mean
absolute error (MAE) of 3.87, a 0.93 Spearman correlation coefficient when
comparing predicted brain age to chronological age, and a -0.29 Spearman
correlation coefficient between the chronological age and the brain age gap
(the absolute value of the gap between the predicted and actual age),
significantly better than previous methods for brain age estimation.

In addition, we subjected our model to Explainable AI (XAI) techniques, such
as 3D occlusion analysis, to validate the high-correlation areas identified by
the model. These areas matched those highlighted by the attention maps
produced by Triamese-ViT, reinforcing the validity of our model’s interpretive
output.

##  2 Related work

This section will introduce the background which is related to our project.
While rapid advancements in brain age estimation have been made possible
through deep learning, the field still faces significant challenges. One key
issue is the lack of comprehensive research into which specific brain regions
most significantly impact age estimation [47]. Limited studies in this area
have indicated the presence of ageism, as accuracy varies across different age
groups [47, 53]. Additionally, deep learning models, often described as ’black
boxes,’ lack transparent processes for mapping input components to output
values [6]. This limitation is particularly evident when compared to the more
interpretable, yet less effective, algorithmic models like decision trees. The
emergence of post-hoc explainability methods, such as SHAP and LIME[36], marks
a crucial development in addressing this issue, though their application in
brain age estimation remains limited. Future research should prioritize the
inclusion of these explainability methods, focusing not just on performance,
but also on reliability and trustworthiness for clinical application, to
ensure that models are both effective and comprehensible in a medical context
[25].

Convolutional neural networks (CNNs) are very popular to be applied for the
field of brain age estimation. Since 2017, CNNs have been at the forefront,
favored for their ability to autonomously extract features and deliver high-
accuracy results [8, 40, 30, 26]. Cole et al. [14] demonstrated this by
training a 3D CNN on MRI samples, finding that models trained on gray matter
outperformed those trained on white matter. To enhance efficiency, researchers
have turned to 2D MRI slices, with [32] using recurrent neural networks (RNN)
to understand the connections between slices, offering a balance between
performance and computational demand. Then, the advanced and special CNN
structure such as SFCN, ResNet, and DenseNet also highlighted their strength
in this area [39, 38, 49]. When the transformer generally become a hot topic,
since it can help pay more attention to details of the images, it is very
suitable for use in the medical field. For instance, the Global–Local
Transformer [24] melds the strengths of CNNs with transformers, utilizing 2D
MRI slices to capture both local and global information in the images. Cai et
al. [9] furthered this approach by leveraging a graph transformer network in a
multimodal method, efficiently harnessing both global and local features for
brain age estimation.

Since we will use XAI techniques to prove our model’s ability to
interpretability, we will introduce the background of XAI. Despite its wide-
ranging utility, AI often operates as a ”black box,” with complex models and a
multitude of parameters obscuring the decision-making process [34, 33, 19, 2].
This lack of transparency is particularly critical in sensitive fields such as
finance [10] and healthcare [11], where understanding AI’s reasoning is
crucial, so this need boosts the development of XAI.

XAI methods generally fall into two categories [37, 17]: backpropagation-based
and perturbation-based. Backpropagation-based XAI involves algorithms that
provide insight during the backpropagation stage of neural network processing,
typically by using derivatives to produce attribution maps like class
activation or saliency maps. On the other hand, perturbation-based approaches
modify the input features in various ways—through occlusion, substitution, or
generative techniques—to observe changes in the output. Our project employs
occlusion analysis, a perturbation-based technique, to validate our model’s
interpretability.

Our Triamese-ViT model is an adaptation of Siamese Networks, a class of neural
networks comprising twin subnetworks that share weights while processing
distinct inputs. Renowned for their performance, Siamese Networks have been
notably utilized in various fields. Zeng et al.[51] employed this network
structure to develop an anchor-free tracking method, leading in scale variance
for video attribute analysis. Zhang et al.[52] applied these networks for
regression tasks on physicochemical datasets, achieving high accuracy.
Additionally, SiamCAR, created by Cui et al. [16] using the Siamese
architecture, excelled in real-time visual tracking of generic objects.

##  3 Method

###  3.1 Proposed Triamese-ViT

In this section, we will introduce the structure of Triamese-ViT. As we show
in Fig.1, the idea of Triamese-ViT is inspired by ViT [20] and Siamese
Networks. The input of Triamese-ViT are 3D MRIs, we called the images
$I$$\in$$R^{H\times W\times C}$, here H, W, C are image height, width and
number of channels. Then we reshape the image $I$ into 3 different viewpoints,
$I$$\to$($I_{x}$,$I_{y}$,$I_{z}$), $I_{x}$$\in$$R^{H\times W\times
C}$,$I_{y}$$\in$$R^{H\times C\times W}$,$I_{z}$$\in$$R^{W\times H\times C}$.
Firstly, we focus on $I_{x}$, we change it to a sequence of flattened 2D
patches called $I_{x,p}$$\in$$R^{N\times(p^{2}\cdot C)}$, the patches are
squares of length $P$. So the number of patches is $N=H\times W/P^{2}$. In all
of the layers in the transformer encoder, vectors with dimension D will be
processed as the object, so we need to map $I_{x}$ to D dimensions with a
trainable linear projection. The formulate is shown below:

$z_{x,0}=Concat(I_{x,class};I_{x,p}^{1}E;I_{x,p}^{2}E;...;I_{x,p}^{N}E)+E_{pos}$
(1)

In Eq. 1, $I_{x,class}$ means a learnable token, in other words, the class
token, which is added into ViT, this method is similar to [18]. Note that
$I_{x,class}$ will finally output from the Transformer Encoder as
$z_{x,L}^{0}$, which represents the image representation $P$ (Eq. 7). $E\in
R^{(p^{2}\cdot C)\times D}$ means Linear Projection, $Concat$ means token
concatenation and $E_{pos}\in R^{(N+1)\times D}$ represents positional
information which is added to each token embedding. We use the same pre-
processing methods on $I_{y}$ and $I_{z}$, and we can get $z_{y,0}$ and
$z_{z,0}$.

Now, the original data is processed into three suitable matrices
$z_{x,0},z_{y,0},z_{z,0}\in R^{(N+1)\times D}$, then we feed them to the
transformer encoder. Each Transformer encoder is a multi-layered system where
the input sequentially passes through several key components in each layer: it
first encounters a Layer Normalization (LN), followed by a Multi-Head
Attention mechanism, another Layer Normalization, and then a Multi-Layer
Perceptron (MLP). The Multi-Head Attention (MSA) operates by conducting
parallel attention calculations across multiple ’heads’, diversifying the
focus and allowing for a richer understanding of the input data.

$[Q,K,V]=FC(z_{x,0})$ (2)

Here, $Q\in R^{(N+1)\times d}$,$K\in R^{(N+1)\times d}$,$V\in R^{(N+1)\times
d}$ represent Query, Keyword and Value. We assume the number of heads in MSA
is n, and $D=n\times d$. Each ’head’ independently processes the input,
allowing the model to simultaneously attend to information from different
representation subspaces at different positions.

$head_{i}=softmax(\frac{Q_{i}K_{i}^{\mathsf{T}}}{\sqrt{d}}V_{i})$ (3)
$MSA(z_{x,0})=Concat(head_{1},head_{2},...,head_{n})$ (4)

Let $z_{x,0}$ be the input of the first layer Transformer Encoder, so the
feedforward calculation in the Encoder is written as:

$z_{x,l}^{{}^{\prime}}=MSA(LN(z_{x,l-1}))+z_{x,l-1}$ (5)
$z_{x,l}=MLP(LN(z_{x,l}^{{}^{\prime}}))+z_{x,l}^{{}^{\prime}}$ (6)

The $l\in[1,2,...L]$. The outputs from each Transformer Encoder are channeled
into an MLP (Multi-Layer Perceptron) head, comprising a hidden layer followed
by an output layer, to generate the final prediction for each view. We denote
the prediction from the $I_{x}$ (the first view) as $P_{x}$. Applying the same
procedure to the $I_{y}$ and $I_{z}$ views, we obtain two additional
predictions, $P_{y}$ and $P_{z}$, corresponding to these orientations.

In the final stage, these three view-based predictions ($P_{x}$, $P_{y}$, and
$P_{z}$) are input into the Triamese MLP. This step integrates the insights
from all three views to produce the model’s comprehensive final prediction.

$P=MLP(P_{x},P_{y},P_{z})$ (7)

We experimented with various strategies to integrate the outputs from the
three viewpoint-specific ViTs. These strategies included averaging the
results, selecting the best-performing output, and combining the outputs
through an MLP. Our experiments demonstrated that processing the results
through an MLP yielded the most effective performance. The loss function used
in the training is the Mean Squared Error (MSE) between the predicted age
$P_{i}$ and chronological age $C_{i}$.

$Loss=\frac{\sum_{i=1}^{n}(P_{i}-C_{i})^{2}}{n}$ (8)

###  3.2 Occlusion Sensitivity Analysis of Triamese-ViT

To demonstrate the interpretability of our Triamese-ViT model, we employ
Occlusion Sensitivity Analysis to generate saliency maps, which we then
compare to the model’s inherent attention maps. In this section, we will delve
into the specifics of how Occlusion Sensitivity Analysis is conducted.

This analysis method systematically obscures different parts of the input data
— in this case, regions of the brain within MRI images — to assess their
impact on the model’s output. By applying a cube-shaped mask (7x7x7 in our
case) that sets the covered voxels to zero, and moving this mask throughout
the entire volume of the brain without overlap, we can monitor changes in the
model’s predictions. The difference in prediction accuracy, quantified by the
Mean Absolute Error (MAE) with and without occlusion, indicates the relative
importance of each region.

By mapping these changes, we create a saliency map that highlights critical
areas the model focuses on for its predictions. Comparing these saliency maps
with the attention maps produced by Triamese-ViT provides a dual perspective
on the model’s decision-making process, offering a clear view of how it
interprets the brain images to estimate age, and proving the ability of
interpretability of the Triamese-ViT. The illustration is shown in Fig.2

Figure 2: Illustration of the framework for occlusion analysis.

##  4 Experiments

###  4.1 Dataset

Our research utilized the IXI https://brain-development.org/ixi-dataset/ and
ABIDE datasets https://fcon_1000.projects.nitrc.org/indi/abide/, encompassing
1351 brain MRI scans from healthy individuals aged 6 to 90 years. The data
allocation was as follows: 70$\%$ for training the Triamese-ViT model, 15$\%$
for validation, and the remaining 15$\%$ for testing to evaluate model
performance.

The MRI scans in our datasets underwent a standardized preprocessing routine
using FSL 5.10 [29]. This process involved several key steps: nonlinear
registration to the standard MNI space, brain extraction [46], and
normalization of voxel values within the brain area (achieved by subtracting
the mean and dividing by the standard deviation of these values). Post-
preprocessing, all MRI scans were unified to a voxel size of $91\times
109\times 91$ and an isotropic spatial resolution of 2mm.

###  4.2 Experimental Settings

In our experiments, all models used the ADAM optimizer on PyTorch, with a
0.001 initial learning rate, $10^{-6}$ decay rate, and $\beta_{1}=0.9$,
$\beta_{2}=0.999$. He initialization and L2 norm weight regularization (weight
$5\times 10^{-4}$) were employed. Batch size was 100, with regularization
weights $\lambda_{1}$ and $\lambda_{2}$ both at 10, a decision informed by the
models’ performance on the validation set, and a manual seed of 3407 due to
research [41].

Data augmentation on MRIs, in order to further mitigate the risk of
overfitting, involved 50$\%$ probability of spatial transformations like 3D
translation (up to 10 voxels), rotation (-20° to 20°), and random flips. For
Triamese-ViT, each ViT had a $7\times 7$ patch size, 768 embedding dimension,
12 attention heads, 10 layers, and 0.1 dropout rate. The MLP inside each ViT
had one hidden layer with 3072 dimensionality. The Triamese MLP’s 9 layers
formed a pyramid sequence (3, 128, 256, 512, 1024, downscaling back to 3),
synthesizing data from the ViTs into a single output.

###  4.3 Performance Evaluation of Age Estimation

In evaluating the performance of our age estimation model, we employ three key
metrics. First is the mean absolute error (MAE) between the predicted age and
chronological age, it is a direct measure of the model’s accuracy; a lower MAE
indicates higher accuracy. The MAE is calculated as follows:

$MAE=\frac{\sum_{i=1}^{n}\left|(P_{i}-C_{i})\right|}{n}$ (9)

The second metric is the correlation coefficient (r), which is computed as the
Spearman correlation between the predicted ages and the chronological ages. A
higher value of ’rp signifies better model performance. The Spearman
correlation is computed using the following formula:

$r=\frac{\sum_{i=1}^{n}(P_{i}-\bar{P})(C_{i}-\bar{C})}{\sqrt{\sum_{i}^{n}(P_{i}-\bar{P})^{2}\sum_{i}^{n}(C_{i}-\bar{C})^{2}}}$
(10)

The last metric is the Spearman correlation between the chronological age and
BAG (the absolute value of the gap between the predicted and actual age) (rp).
This metric evaluates the fairness of the model, particularly checking for age
bias in predictions. A higher correlation in this context suggests more
pronounced ageism. The formula for this metric is:

$G_{i}=\left|P_{i}-C_{i}\right|$ (11)
$rp=\frac{\sum_{i=1}^{n}(G_{i}-\bar{G})(C_{i}-\bar{C})}{\sqrt{\sum_{i}^{n}(G_{i}-\bar{G})^{2}\sum_{i}^{n}(C_{i}-\bar{C})^{2}}}$
(12)

###  4.4 Comparison With State-of-the-Art Algorithms for Brain Age Estimation

In this section, we’ll present a comparative analysis of our Triamese-ViT
model against various other state-of-the-art brain age estimation algorithms.
This comparison, based on their performance on our dataset, aims to highlight
the advancements and superior capabilities of our proposed model.

Algorithm | MAE | r | rp
---|---|---|---
ScaleDense [12] | 3.92 | 0.92 | 0.39
5-layer CNN [15] | 4.51 | 0.78 | 0.45
ResNet [14] | 4.01 | 0.83 | -0.31
VGG19 [28] | 4.07 | 0.7 | 0.5
VGG16 [28] | 5.24 | 0.6 | 0.42
Global-Local Transformer [24] | 4.66 | 0.78 | -0.32
Efficient Net [42] | 4.56 | 0.89 | -0.4
Our Triamese-ViT | $\mathbf{3.87}$ | $\mathbf{0.93}$ | $\mathbf{-0.29}$
Table 1:  The details of tested algorithms’ performance. Since the input of
Global-Local Transformer should be 2D image, we extract 2D slices around the
center of the 3D brain volumes in the axial as input, which is the same
process method as [24]. Other algorithms’ input are 3D MRIs with dimensions
(91,109,91). Our u-DemAI has consistently achieved the best among all
measures.

Table 1 in our study provides a detailed comparison of our Triamese-ViT model
with seven other models, encompassing both classic and contemporary approaches
to brain age estimation. The comparison includes five well-known 3D CNN-based
models: ScaleDense, a 5-layer CNN, ResNet, VGG16, and VGG19. Additionally, we
evaluated against two other high-performing methods: the Global-Local
Transformer, which is trained on 2D slices of the brain, and Effificent Net,
known for its ensemble structure.

According to Table 1, our Triamese-ViT model leads in Mean Absolute Error
(MAE) performance with a score of 3.87. ScaleDense also performs impressively,
achieving an MAE of 3.92, followed by ResNet with 4.01. The highest MAE,
indicating the least accuracy, was recorded for VGG16 at 5.24.

In terms of the Spearman Correlation between predicted and chronological ages,
our Triamese-ViT tops the list with a correlation of 0.93. ScaleDense is close
behind with a correlation of 0.92, and ResNet follows with 0.83. VGG16 trails
in this metric as well, showing a correlation of only 0.6.

Regarding the Spearman correlation between the Brain Age Gap (BAG) and
chronological age, which reflects the fairness of the model, Triamese-ViT
demonstrates the best outcome with a -0.29 correlation, indicating lower age
bias. ResNet is next with a -0.31 correlation. On the other hand, VGG19 shows
the most pronounced age bias with a 0.5 positive correlation.

Overall, this comparison underscores the effectiveness of Triamese-ViT in
brain age estimation, both in terms of accuracy and fairness, when benchmarked
against other leading methods in the field.

###  4.5 Ablation Experiments

In this part of our study, we conduct ablation experiments to explore and
justify the design choices in the structure of Triamese-ViT. First of all, we
specifically focus on the number of layers in the Triamese MLP. While keeping
all other variables constant, we varied the number of MLP layers and observed
their impact on the model’s performance.

The findings depicted in Fig. 3 show a distinctive trend in the Mean Absolute
Error (MAE) relative to the MLP layers in Triamese MLP. The MAE initially
rises when increasing layers from 4 to 6, then decreases after 6 layers,
reaching a minimum at 9 layers before rising again at 10 layers. This
indicates an optimal layer count for balancing model complexity and accuracy.
The observed MAE variation with different layer counts underscores the
intricate relationship between model depth and performance, emphasizing the
need for precise architectural tuning in the model.

Figure 3: The impact of the number of MLP layers in Triamese-Encoder

Continuing our ablation studies, we turned our focus to the backbone of
Triamese-ViT. To assess the impact of different backbone architectures, we
substituted the original ViT with alternative models like ResNet, a 5-layer
CNN, and VGG19. These were then integrated with the Triamese MLP to evaluate
how they influenced overall performance. The results of this experiment are
detailed in Fig. 4.

For ease of interpretation, we chose to display the absolute value of the
Spearman correlation coefficient between BAG and chronological age (denoted as
$\left\|rp\right\|$) in the figure. A larger value of $\left\|rp\right\|$
indicates a stronger age bias in the model’s predictions. According to our
findings, the original ViT backbone proves to be the most effective for the
Triamese structure. The 5-layer CNN also shows commendable adaptability,
registering an MAE of 6, a Spearman correlation (r) of 0.85, and
$\left\|rp\right\|$ of 0.45.

In stark contrast, ResNet and VGG19 appear significantly less suited for the
Triamese framework. Both these architectures yielded MAEs exceeding 10, which
are highly unfavorable outcomes for brain age estimation. This experiment
underscores the importance of selecting an appropriate backbone model for the
Triamese structure to ensure optimal performance.

Figure 4: The impact of the backbone architectures

Next, we explore the unique structures within our Triamese-ViT model,
particularly focusing on the individual contributions of the three Vision
Transformers (ViTs) oriented along different axes of the MRIs. These are the
$ViT_{x}$ with dimensions (91,109,91), $ViT_{y}$ with dimensions (91,91,109),
and $ViT_{z}$ with dimensions (109,91,91). The performance of each of these
orientation-specific ViTs is crucial in understanding the efficacy of the
combined Triamese-MLP structure.

Additionally, we tested a variant model, $Triamese_{map}$, which also utilizes
three ViTs on different viewpoints. However, unlike the standard Triamese-ViT,
each ViT in $Triamese_{map}$ outputs a feature map from the Transformer
Encoder, rather than a direct prediction from the MLP Head. The Triamese MLP
in this variant then takes as input the concatenated feature maps from the
three ViTs to make the final prediction.

The comparative performance of these models, including each individual
orientation ViT and the $Triamese_{map}$ variant, is presented in Table 2.
This comparison is key to demonstrating the value added by the Triamese MLP in
synthesizing the perspectives from the three distinct ViT orientations,
highlighting the importance of integrating these viewpoints for more accurate
brain age estimation.

Algorithm | MAE | r | rp
---|---|---|---
Triamese-ViT | $\mathbf{3.87}$ | $\mathbf{0.93}$ | -$\mathbf{0.29}$
$ViT_{x}$ | 4.42 | 0.78 | 0.33
$ViT_{y}$ | 4.99 | 0.92 | -0.29
$ViT_{z}$ | 5.29 | 0.73 | -0.37
$ViT_{map}$ | 5.04 | 0.61 | -0.55
Table 2: Unique structures performance

Table 2 says Triamese MLP supports a great improvement of performance, for
MAE, $ViT_{x}$ is the second best with 4.42, $ViT_{z}$ is the worst with 5.29.
As for r, $ViT_{y}$ has the highest value with 0.92, This is closely followed
by the combined Triamese-ViT model. Notably, $ViT_{map}$, which uses
concatenated feature maps for prediction, shows the lowest correlation value
at 0.61. Regarding the aspect of fairness, only $ViT_{map}$ displays a strong
negative correlation. This suggests a significant reduction in age bias.
Conversely, the other models, including the individual orientation-specific
ViTs, exhibit minimal ageism in their predictions.

Overall, the data in Table 2 strongly supports the efficacy of the Triamese
MLP in enhancing both the accuracy and fairness of brain age estimation,
validating the design of our Triamese-ViT model.

###  4.6 Explainable Results for Triamese-ViT

Figure 5: Comparison between the attention map and occlusion analysis from
Triamese-ViT. The upper half showcases the results from the occlusion
analysis, while the lower half displays the results from the attention map.
Both halves collectively highlight the specific regions of the brain that the
Triamese-ViT model prioritizes and considers most informative for determining
age.

In this section, we delve into the explainable results generated by the
Triamese-ViT model. Fig. 5 displays outcomes obtained through two different
methods: occlusion analysis and the attention map feature of the Triamese-ViT.

The upper half of the figure, representing the occlusion analysis, indicates
that the removal of the Basal Ganglia and Thalamus significantly impacts the
model’s predictions. This finding underscores the importance of these brain
areas in age estimation. Additionally, the analysis suggests that the left
side of the brain has a more pronounced effect on the results.

The lower half of the figure, showing the attention map results, offers more
nuanced insights. It corroborates the significance of the Basal Ganglia and
Thalamus but also highlights the Midbrain as a key influencer. This detailed
analysis reaffirms the greater importance of the left brain in the prediction
process.

Comparing these two methods, it is evident that the Triamese-ViT’s attention
maps not only align with the findings from traditional explainable AI (XAI)
methods but also provide additional, detailed insights.

Supporting this, various studies have linked critical brain functions and
diseases to these regions. [48, 50] have associated conditions like
Parkinson’s disease, Tourette’s syndrome, Huntington’s disease, Alzheimer’s
Disease, and addiction with the Basal Ganglia. [43, 45] discuss the role of
the substantia nigra in the midbrain, a dopamine-producing area crucial for
movement regulation and significantly affected in Parkinson’s disease.
Furthermore, the Thalamus, as noted in [44], acts as a central hub for sensory
information processing and plays a vital role in attention coordination.

These scientific findings validate the attention map’s emphasis on the Basal
Ganglia, Midbrain, and Thalamus. Their crucial roles in various brain diseases
directly relate to brain age estimation. Thus, the results from our Triamese-
ViT not only demonstrate its strong capability for explainability but also
have significant implications for understanding and studying brain diseases.

##  5 Discussion

In our research, we introduce a groundbreaking deep-learning algorithm,
Triamese-ViT, used for brain age estimation. This model is benchmarked against
other leading models in the field, demonstrating its remarkable superiority in
performance. The pivotal contribution of Triamese-ViT is its innovative
Triamese structure. This design is a first in the realm of brain age
estimation, merging the benefits of global context understanding with detailed
image analysis. It delves into the intricate relationships between image
patches, leading to predictions that are not only more comprehensive and
accurate but also highly interpretable. This blend of detailed analysis and
contextual understanding sets Triamese-ViT apart, marking a significant
advancement in the field of brain age estimation.

In this project, while traditional CNN-based models use complete 3D MRI scans
for detailed predictions, their focus on small, localized pixels can miss
crucial global features, affecting prediction accuracy, especially in complex
brain structure analysis. Conversely, Vision Transformers (ViTs), adapted from
natural language processing, enhance brain age estimation by dissecting images
into patches and analyzing their interrelations, offering detailed insights.
However, ViTs often overlook the overall image context and typically process
3D MRIs as 2D slices, which may result in the loss of important depth
information.

Our innovation, Triamese-ViT, draws inspiration from the Siamese Network,
which shows great performance in various fileds. Triamese-ViT is constructed
with three ViTs, each analyzing 3D images from different viewpoints. This
setup, combined with a Triamese MLP for feature extraction and prediction,
effectively harnesses the strengths of both CNNs and ViTs while mitigating
their respective weaknesses.

The ViT backbone enables detailed image analysis and understanding of inter-
patch relationships. The Triamese structure, on the other hand, ensures a
comprehensive assessment of the whole image from multiple perspectives,
preserving the depth aspect in the estimation process.

Triamese-ViT, when tested on a public dataset, demonstrated excellent
performance: a Mean Absolute Error (MAE) of 3.87, a 0.93 Spearman correlation
with chronological age, and a -0.29 Spearman correlation between Brain Age Gap
(BAG) and chronological age. These results signify not only high predictive
accuracy but also a reduction in age bias, marking a notable advancement in
brain age estimation.

Moreover, Triamese-ViT excels in interpretability, crucial in medicine. Its
attention maps provide more detailed insights compared to occlusion analysis
and are integrated into the prediction process, offering a faster, user-
friendly interpretation. This feature is especially valuable in medical
settings where swift, accurate decision-making is essential.

The attention maps generated by Triamese-ViT pinpoint the Basal Ganglia,
Thalamus, and Midbrain as crucial areas for brain age estimation. These
regions play a significant role in determining the health of a patient’s brain
according to the model. Supporting this, several medical studies [44, 45, 43,
50, 48] have established a strong correlation between these brain areas and
various severe neurological diseases, underlining the model’s robust
interpretability.

Furthermore, both the attention maps and occlusion analysis consistently
indicate a greater influence of the left brain in age estimation. This finding
could be attributed to the larger proportion of right-handed individuals in
the dataset, as right-handedness is often associated with more developed left-
brain regions. This aspect of our findings opens avenues for further research
and underscores the depth and reliability of the insights provided by
Triamese-ViT, making it an invaluable tool for advancing our understanding of
brain health and aging.

##  6 Conclusion

In this paper, we introduced Triamese-ViT, a novel deep-learning architecture
applied to brain age estimation. Triamese-ViT exhibits exceptional accuracy,
fairness, and interpretability, surpassing existing advanced algorithms in the
field. Its high performance, low bias, and robust interpretability make it
well-suited for medical research. The model’s user-friendly nature enhances
its applicability in clinical settings where efficiency and clarity are
crucial. Triamese-ViT represents a meaningful contribution to the integration
of AI in medicine, offering potential advancements in research and
applications. We envision its utility not only in advancing brain age
estimation but also as a valuable tool for broader medical AI research and
development.

## References

  * Al-Hammuri et al. [2023] Khalid Al-Hammuri, Fayez Gebali, Awos Kanan, and Ilamparithi Thirumarai Chelvan.  Vision transformer architecture and applications in digital health: a tutorial and survey.  _Visual Computing for Industry, Biomedicine, and Art_ , 6(1):14, 2023.
  * Ali et al. [2023] Sajid Ali, Tamer Abuhmed, Shaker El-Sappagh, Khan Muhammad, Jose M Alonso-Moral, Roberto Confalonieri, Riccardo Guidotti, Javier Del Ser, Natalia Díaz-Rodríguez, and Francisco Herrera.  Explainable artificial intelligence (xai): What we know and what is left to attain trustworthy artificial intelligence.  _Information Fusion_ , 99:101805, 2023.
  * Alzubaidi et al. [2021] Laith Alzubaidi, Jinglan Zhang, Amjad J Humaidi, Ayad Al-Dujaili, Ye Duan, Omran Al-Shamma, José Santamaría, Mohammed A Fadhel, Muthana Al-Amidie, and Laith Farhan.  Review of deep learning: Concepts, cnn architectures, challenges, applications, future directions.  _Journal of big Data_ , 8:1–74, 2021.
  * Armanious et al. [2021] Karim Armanious, Sherif Abdulatif, Wenbin Shi, Shashank Salian, Thomas Küstner, Daniel Weiskopf, Tobias Hepp, Sergios Gatidis, and Bin Yang.  Age-net: An mri-based iterative framework for brain biological age estimation.  _IEEE Transactions on Medical Imaging_ , 40(7):1778–1791, 2021.
  * Avants et al. [2008] Brian B Avants, Charles L Epstein, Murray Grossman, and James C Gee.  Symmetric diffeomorphic image registration with cross-correlation: evaluating automated labeling of elderly and neurodegenerative brain.  _Medical image analysis_ , 12(1):26–41, 2008.
  * Barredo Arrieta et al. [2020] Alejandro Barredo Arrieta, Natalia Díaz-Rodríguez, Javier Del Ser, Adrien Bennetot, Siham Tabik, Alberto Barbado, Salvador Garcia, Sergio Gil-Lopez, Daniel Molina, Richard Benjamins, Raja Chatila, and Francisco Herrera.  Explainable artificial intelligence (xai): Concepts, taxonomies, opportunities and challenges toward responsible ai.  _Information Fusion_ , 58:82–115, 2020.
  * Beheshti et al. [2018] Iman Beheshti, Norihide Maikusa, and Hiroshi Matsuda.  The association between “brain-age score”(bas) and traditional neuropsychological screening tools in alzheimer’s disease.  _Brain and Behavior_ , 8(8):e01020, 2018.
  * Bellantuono et al. [2021] Loredana Bellantuono, Luca Marzano, Marianna La Rocca, Dominique Duncan, Angela Lombardi, Tommaso Maggipinto, Alfonso Monaco, Sabina Tangaro, Nicola Amoroso, and Roberto Bellotti.  Predicting brain age with complex networks: From adolescence to adulthood.  _NeuroImage_ , 225:117458, 2021.
  * Cai et al. [2022] Hongjie Cai, Yue Gao, and Manhua Liu.  Graph transformer geometric learning of brain networks using multimodal mr images for brain age estimation.  _IEEE Transactions on Medical Imaging_ , 42(2):456–466, 2022.
  * Cao [2022] Longbing Cao.  Ai in finance: challenges, techniques, and opportunities.  _ACM Computing Surveys (CSUR)_ , 55(3):1–38, 2022.
  * Caruana et al. [2015] Rich Caruana, Yin Lou, Johannes Gehrke, Paul Koch, Marc Sturm, and Noemie Elhadad.  Intelligible models for healthcare: Predicting pneumonia risk and hospital 30-day readmission.  In _Proceedings of the 21th ACM SIGKDD international conference on knowledge discovery and data mining_ , pages 1721–1730, 2015.
  * Cheng et al. [2021] Jian Cheng, Ziyang Liu, Hao Guan, Zhenzhou Wu, Haogang Zhu, Jiyang Jiang, Wei Wen, Dacheng Tao, and Tao Liu.  Brain age estimation from mri using cascade networks with ranking loss.  _IEEE Transactions on Medical Imaging_ , 40(12):3400–3412, 2021.
  * Chung et al. [2018] Yoonho Chung, Jean Addington, Carrie E Bearden, Kristin Cadenhead, Barbara Cornblatt, Daniel H Mathalon, Thomas McGlashan, Diana Perkins, Larry J Seidman, Ming Tsuang, et al.  Use of machine learning to determine deviance in neuroanatomical maturity associated with future psychosis in youths at clinically high risk.  _JAMA psychiatry_ , 75(9):960–968, 2018.
  * Cole and Franke [2017] James H Cole and Katja Franke.  Predicting age using neuroimaging: innovative brain ageing biomarkers.  _Trends in neurosciences_ , 40(12):681–690, 2017.
  * Couvy-Duchesne et al. [2020] Baptiste Couvy-Duchesne, Johann Faouzi, Benoît Martin, Elina Thibeau-Sutre, Adam Wild, Manon Ansart, Stanley Durrleman, Didier Dormont, Ninon Burgos, and Olivier Colliot.  Ensemble learning of convolutional neural network, support vector machine, and best linear unbiased predictor for brain age prediction: Aramis contribution to the predictive analytics competition 2019 challenge.  _Frontiers in Psychiatry_ , 11:593336, 2020.
  * Cui et al. [2022] Ying Cui, Dongyan Guo, Yanyan Shao, Zhenhua Wang, Chunhua Shen, Liyan Zhang, and Shengyong Chen.  Joint classification and regression for visual tracking with fully convolutional siamese networks.  _International Journal of Computer Vision_ , pages 1–17, 2022.
  * Das and Rad [2020] Arun Das and Paul Rad.  Opportunities and challenges in explainable artificial intelligence (xai): A survey.  _arXiv preprint arXiv:2006.11371_ , 2020.
  * Devlin et al. [2018] Jacob Devlin, Ming-Wei Chang, Kenton Lee, and Kristina Toutanova.  Bert: Pre-training of deep bidirectional transformers for language understanding.  _arXiv preprint arXiv:1810.04805_ , 2018.
  * Doshi-Velez and Kim [2017] Finale Doshi-Velez and Been Kim.  Towards a rigorous science of interpretable machine learning.  _arXiv preprint arXiv:1702.08608_ , 2017.
  * Dosovitskiy et al. [2020] Alexey Dosovitskiy, Lucas Beyer, Alexander Kolesnikov, Dirk Weissenborn, Xiaohua Zhai, Thomas Unterthiner, Mostafa Dehghani, Matthias Minderer, Georg Heigold, Sylvain Gelly, et al.  An image is worth 16x16 words: Transformers for image recognition at scale.  _arXiv preprint arXiv:2010.11929_ , 2020.
  * Feng et al. [2020] Xinyang Feng, Zachary C Lipton, Jie Yang, Scott A Small, Frank A Provenzano, Alzheimer’s Disease Neuroimaging Initiative, Frontotemporal Lobar Degeneration Neuroimaging Initiative, et al.  Estimating brain age based on a uniform healthy population with deep learning and structural magnetic resonance imaging.  _Neurobiology of aging_ , 91:15–25, 2020.
  * Gaser et al. [2013] Christian Gaser, Katja Franke, Stefan Klöppel, Nikolaos Koutsouleris, Heinrich Sauer, and Alzheimer’s Disease Neuroimaging Initiative.  Brainage in mild cognitive impaired patients: predicting the conversion to alzheimer’s disease.  _PloS one_ , 8(6):e67346, 2013.
  * Han et al. [2021] Laura KM Han, Richard Dinga, Tim Hahn, Christopher RK Ching, Lisa T Eyler, Lyubomir Aftanas, Moji Aghajani, André Aleman, Bernhard T Baune, Klaus Berger, et al.  Brain aging in major depressive disorder: results from the enigma major depressive disorder working group.  _Molecular psychiatry_ , 26(9):5124–5139, 2021.
  * He et al. [2021] Sheng He, P Ellen Grant, and Yangming Ou.  Global-local transformer for brain age estimation.  _IEEE transactions on medical imaging_ , 41(1):213–224, 2021.
  * Holzinger et al. [2022] Andreas Holzinger, Matthias Dehmer, Frank Emmert-Streib, Rita Cucchiara, Isabelle Augenstein, Javier Del Ser, Wojciech Samek, Igor Jurisica, and Natalia Díaz-Rodríguez.  Information fusion as an integrative cross-cutting enabler to achieve robust, explainable, and trustworthy medical artificial intelligence.  _Information Fusion_ , 79:263–278, 2022.
  * Hong et al. [2020] Jin Hong, Zhangzhi Feng, Shui-Hua Wang, Andrew Peet, Yu-Dong Zhang, Yu Sun, and Ming Yang.  Brain age prediction of children using routine brain mr images via deep learning.  _Frontiers in Neurology_ , 11:584682, 2020.
  * Hu et al. [2022] Yixiao Hu, Haolin Wang, and Baobin Li.  Sqet: Squeeze and excitation transformer for high-accuracy brain age estimation.  In _2022 IEEE International Conference on Bioinformatics and Biomedicine (BIBM)_ , pages 1554–1557. IEEE, 2022.
  * Huang et al. [2017] Tzu-Wei Huang, Hwann-Tzong Chen, Ryuichi Fujimoto, Koichi Ito, Kai Wu, Kazunori Sato, Yasuyuki Taki, Hiroshi Fukuda, and Takafumi Aoki.  Age estimation from brain mri images using deep learning.  In _2017 IEEE 14th International Symposium on Biomedical Imaging (ISBI 2017)_ , pages 849–852. IEEE, 2017.
  * Jenkinson et al. [2012] Mark Jenkinson, Christian F Beckmann, Timothy EJ Behrens, Mark W Woolrich, and Stephen M Smith.  Fsl.  _Neuroimage_ , 62(2):782–790, 2012.
  * Jiang et al. [2020] Huiting Jiang, Na Lu, Kewei Chen, Li Yao, Ke Li, Jiacai Zhang, and Xiaojuan Guo.  Predicting brain age of healthy adults based on structural mri parcellation using convolutional neural networks.  _Frontiers in neurology_ , 10:1346, 2020.
  * Khan et al. [2022] Salman Khan, Muzammal Naseer, Munawar Hayat, Syed Waqas Zamir, Fahad Shahbaz Khan, and Mubarak Shah.  Transformers in vision: A survey.  _ACM computing surveys (CSUR)_ , 54(10s):1–41, 2022.
  * Lam et al. [2020] Pradeep K Lam, Vigneshwaran Santhalingam, Parth Suresh, Rahul Baboota, Alyssa H Zhu, Sophia I Thomopoulos, Neda Jahanshad, and Paul M Thompson.  Accurate brain age prediction using recurrent slice-based networks.  In _16th International Symposium on Medical Information Processing and Analysis_ , pages 11–20. SPIE, 2020.
  * Langer et al. [2021a] Markus Langer, Kevin Baum, Kathrin Hartmann, Stefan Hessel, Timo Speith, and Jonas Wahl.  Explainability auditing for intelligent systems: a rationale for multi-disciplinary perspectives.  In _2021 IEEE 29th international requirements engineering conference workshops (REW)_ , pages 164–168. IEEE, 2021a.
  * Langer et al. [2021b] Markus Langer, Daniel Oster, Timo Speith, Holger Hermanns, Lena Kästner, Eva Schmidt, Andreas Sesing, and Kevin Baum.  What do we want from explainable artificial intelligence (xai)?–a stakeholder perspective on xai and a conceptual model guiding interdisciplinary xai research.  _Artificial Intelligence_ , 296:103473, 2021b.
  * Liu et al. [2015] Xiaoxiao Liu, Marc Niethammer, Roland Kwitt, Nikhil Singh, Matt McCormick, and Stephen Aylward.  Low-rank atlas image analyses in the presence of pathologies.  _IEEE transactions on medical imaging_ , 34(12):2583–2591, 2015.
  * Lombardi et al. [2021] Angela Lombardi, Domenico Diacono, Nicola Amoroso, Alfonso Monaco, João Manuel RS Tavares, Roberto Bellotti, and Sabina Tangaro.  Explainable deep learning for personalized age prediction with brain morphology.  _Frontiers in neuroscience_ , 15:578, 2021.
  * Longo et al. [2023] Luca Longo, Mario Brcic, Federico Cabitza, Jaesik Choi, Roberto Confalonieri, Javier Del Ser, Riccardo Guidotti, Yoichi Hayashi, Francisco Herrera, Andreas Holzinger, et al.  Explainable artificial intelligence (xai) 2.0: A manifesto of open challenges and interdisciplinary research directions.  _arXiv preprint arXiv:2310.19775_ , 2023.
  * Mouches et al. [2022] Pauline Mouches, Matthias Wilms, Deepthi Rajashekar, Sönke Langner, and Nils D Forkert.  Multimodal biological brain age prediction using magnetic resonance imaging and angiography with the identification of predictive regions.  _Human brain mapping_ , 43(8):2554–2566, 2022.
  * Ning et al. [2021] Kaida Ning, Ben A Duffy, Meredith Franklin, Will Matloff, Lu Zhao, Nibal Arzouni, Fengzhu Sun, and Arthur W Toga.  Improving brain age estimates with deep learning leads to identification of novel genetic factors associated with brain aging.  _Neurobiology of Aging_ , 105:199–204, 2021.
  * Pardakhti and Sajedi [2020] Nastaran Pardakhti and Hedieh Sajedi.  Brain age estimation based on 3d mri images using 3d convolutional neural network.  _Multimedia Tools and Applications_ , 79:25051–25065, 2020.
  * Picard [2021] David Picard.  Torch. manual_seed (3407) is all you need: On the influence of random seeds in deep learning architectures for computer vision.  _arXiv preprint arXiv:2109.08203_ , 2021.
  * Poloni et al. [2022] Katia Maria Poloni, Ricardo José Ferrari, Alzheimer’s Disease Neuroimaging Initiative, et al.  A deep ensemble hippocampal cnn model for brain age estimation applied to alzheimer’s diagnosis.  _Expert Systems with Applications_ , 195:116622, 2022.
  * Russo and Riessland [2022] Taylor Russo and Markus Riessland.  Age-related midbrain inflammation and senescence in parkinson’s disease.  _Frontiers in Aging Neuroscience_ , 14:917797, 2022.
  * Shine et al. [2023] James M Shine, Laura D Lewis, Douglas D Garrett, and Kai Hwang.  The impact of the human thalamus on brain-wide information processing.  _Nature Reviews Neuroscience_ , pages 1–15, 2023.
  * Smajić et al. [2022] Semra Smajić, Cesar A Prada-Medina, Zied Landoulsi, Jenny Ghelfi, Sylvie Delcambre, Carola Dietrich, Javier Jarazo, Jana Henck, Saranya Balachandran, Sinthuja Pachchek, et al.  Single-cell sequencing of human midbrain reveals glial activation and a parkinson-specific neuronal state.  _Brain_ , 145(3):964–978, 2022.
  * Smith [2002] Stephen M Smith.  Fast robust automated brain extraction.  _Human brain mapping_ , 17(3):143–155, 2002.
  * Tanveer et al. [2023] M Tanveer, MA Ganaie, Iman Beheshti, Tripti Goel, Nehal Ahmad, Kuan-Ting Lai, Kaizhu Huang, Yu-Dong Zhang, Javier Del Ser, and Chin-Teng Lin.  Deep learning for brain age estimation: A systematic review.  _Information Fusion_ , 2023.
  * Wang et al. [2022] Mengya Wang, Huayuan Liu, and Zegang Ma.  Roles of the cannabinoid system in the basal ganglia in parkinson’s disease.  _Frontiers in cellular neuroscience_ , 16:832854, 2022.
  * Wood et al. [2022] David A Wood, Sina Kafiabadi, Ayisha Al Busaidi, Emily Guilhem, Antanas Montvila, Jeremy Lynch, Matthew Townend, Siddharth Agarwal, Asif Mazumder, Gareth J Barker, et al.  Accurate brain-age models for routine clinical mri examinations.  _Neuroimage_ , 249:118871, 2022.
  * Xiong et al. [2022] Yu Xiong, Chenghui Ye, Ying Chen, Xiaochun Zhong, Hongda Chen, Ruxin Sun, Jiaqi Zhang, Zhanhua Zhong, and Min Huang.  Altered functional connectivity of basal ganglia in mild cognitive impairment and alzheimer’s disease.  _Brain Sciences_ , 12(11):1555, 2022.
  * Zeng et al. [2022] Yulin Zeng, Bi Zeng, Xiuwen Yin, and Guangke Chen.  Siampcf: siamese point regression with coarse-fine classification network for visual tracking.  _Applied Intelligence_ , 52(5):4973–4986, 2022.
  * Zhang et al. [2023] Yumeng Zhang, Janosch Menke, Jiazhen He, Eva Nittinger, Christian Tyrchan, Oliver Koch, and Hongtao Zhao.  Similarity-based pairing improves efficiency of siamese neural networks for regression tasks and uncertainty quantification.  _Journal of Cheminformatics_ , 15(1):75, 2023.
  * [53] Zhaonian Zhang and Richard Jiang.  User-centric democratization towards social value aligned medical ai services.
  * Zhang et al. [2022] Zhaonian Zhang, Richard Jiang, Ce Zhang, Bryan Williams, Ziping Jiang, Chang-Tsun Li, Paul Chazot, Nicola Pavese, Ahmed Bouridane, and Azeddine Beghdadi.  Robust brain age estimation based on smri via nonlinear age-adaptive ensemble learning.  _IEEE Transactions on Neural Systems and Rehabilitation Engineering_ , 30:2146–2156, 2022.


# Finite-size effects in the diffusion dynamics of a glassforming binary
mixture with large size ratio

Vinay Vaibhav  The Institute of Mathematical Sciences, IV Cross Road, CIT
Campus, Taramani, Chennai 600 113, Tamil Nadu, India  Jürgen Horbach  Institut
für Theoretische Physik II: Weiche Materie, Heinrich-Heine-Universität
Düsseldorf, Universitätsstraße 1, 40225 Düsseldorf, Germany  Pinaki Chaudhuri
The Institute of Mathematical Sciences, IV Cross Road, CIT Campus, Taramani,
Chennai 600 113, Tamil Nadu, India

###### Abstract

Extensive molecular dynamics (MD) computer simulations of an equimolar
glassforming AB mixture with large size ratio are presented. While the large A
particles show a glass transition around the critical density of mode coupling
theory $\rho_{c}$, the small B particles remain mobile with a relatively weak
decrease of their self-diffusion coefficient $D_{\rm B}$ with increasing
density. Surprisingly, around $\rho_{c}$, the self-diffusion coefficient of
the A particles, $D_{\rm A}$, also starts to follow a rather weak dependence
on density. We show that this is due to finite-size effects that can be
understood from the analysis of the collective interdiffusion dynamics.

##  I Introduction

Many soft matter systems as well as many biological systems consist of
particles of very different sizes bechinger2013 ; weiss2014 ; hoefling2013 .
These systems may show a glassy dynamics with a time-scale separation of
relaxation processes among the different constituents. Examples of such
systems are glassforming mixtures of small and large particles that have been
studied experimentally via various colloidal and organic systems imhof1995_1 ;
imhof1995_2 ; kurita2010 ; blochowicz2012 ; bierwirth2018 ; sentjabrskaja2016
; laurati2019  and numerically via hard or soft sphere systems in computer
simulations moreno2006 ; horbach2009 ; xu2012 ; xu2015 ; lazaro2019  as well
as in the framework of mode-coupling theory (MCT) bosse1987 ; bosse1995 ;
voigtmann2011 . A common feature in these studies is a freezing of the large
particles into a glass state while the small particles remain mobile. Here,
the dynamics of the small particles is typically associated with anomalous
diffusion on long transient time scales, as reflected, e.g., by a sublinear
growth of the mean-squared displacement $\delta r^{2}(t)$ as a function of
time $t$, i.e. $\delta r^{2}(t)\propto t^{\alpha}$ with $\alpha<1$. In
computer simulations as well as experiments of disparate-sized mixtures
kurita2010 ; blochowicz2012 ; horbach2009 ; schnyder2018 ; kurzidim2011 , one
finds values for the exponent $\alpha$ that in general depend on the
temperature, the total density of the system, the concentration of small
mobile particles, and the interactions between the particle, especially those
between the large and the small particles. Thus, the values of $\alpha$ are
non-universal and there is typically the lack of a sharp critical point at
which one observes an asymptotic subdiffusive behavior in the longtime limit
with a universal value of the exponent $\alpha$. This non-universal behavior
can be due to the thermal motion of the particles or soft interactions between
small and large particles.

In a binary mixture of small and large particles, there are the two
corresponding selfdiffusion coefficients $D_{\rm s}$ and $D_{\rm l}$,
respectively, that characterize on one hand the glassy dynamics of the large
particles ($D_{\rm l}$) and on the other hand the transport of the mobile
small particles ($D_{\rm s}$). However, in addition to these single-particle
transport coefficients, there is also a collective diffusion coefficient,
namely the interdiffusion coefficient $D_{\rm AB}$, that characterizes the
mass transport in the binary mixture fitts1962 ; akcasu1997 ; horbach2007 . In
good approximation, $D_{\rm AB}$ can be often expressed as a simple linear
combination of the selfdiffusion coefficients,

$D_{\rm AB}=\Phi\left(x_{\rm l}D_{\rm s}+x_{\rm s}D_{\rm l}\right),$ (1)

with $x_{\rm l}$ and $x_{\rm s}$ the concentration of the large and the small
particles, respectively, and $\Phi$ the thermodynamic factor (see below).
Equation (1) is often called the Darken equation darken1949  or the Hartley-
Crank equation hartley1949 . Computer simulations of glassforming metallic
systems Al-Ni and Zr-Ni with different compositions horbach2007 ; kuhn2014
have shown that Eq. (1) qualitatively reproduces the temperature dependence of
the interdiffusion coefficient, especially at low temperatures.

The question of whether the interdiffusion coefficient can be expressed in
terms of the selfdiffusion coefficients has been extensively discussed in the
literature, especially in the context of (binary) polymer mixtures akcasu1997
; bearman1960 ; brochard1986 ; sillescu1987 ; akcasu1991 ; hess1990 . In this
context, Eq. (1) is often referred to as the result of a “fast mode theory”
akcasu1997 , because according to Eq. (1) for a disparate-sized binary mixture
the interdiffusion coefficient would be essentially given by the selfdiffusion
coefficient, $D_{\rm s}$, of the fast mobile species. In a “slow mode theory”,
however, the opposite behavior is predicted. Here, the relation between the
interdiffusion and the selfdiffusion coefficients is given by hess1990

$D_{\rm AB}=\frac{\Phi}{\frac{x_{\rm l}}{D_{\rm s}}+\frac{x_{\rm s}}{D_{\rm
l}}},$ (2)

This result can be obtained in the framework of a random phase approximation
akcasu1997 . It implies that in a disparate-sized mixture $D_{\rm AB}$ is
dominated by the selfdiffusion coefficient of the slow species, $D_{\rm l}$.
Note that in the framework of MCT one also finds that $D_{\rm AB}$ tends to
“follow” the slow species such that it always vanishes in a glass state
latz1990 .

For our study, we consider an equimolar binary AB mixture of soft spheres for
which the size ratio of the two species is $\approx 2.85$ and, in addition,
the strength of the interaction between AB pairs is weaker than that between
AA and BB pairs. In an earlier molecular dynamics (MD) simulation study of
this system horbach2009 , it has been demonstrated that on the typical time
scale accessible in the MD simulation the A species falls out-of-equilibrium
around the MCT critical number density which is at $\rho_{c}=2.23$,
corresponding to a number density of A particles $\rho_{c}^{\rm A}=1.115$.
While the A species is in a frozen-in state above $\rho_{c}$, the B species
remains mobile and there is a relatively weak decrease of the corresponding
selfdiffusion coefficient $D_{\rm B}$ with increasing density above
$\rho_{c}$. We demonstrate that Eq. (1) very well describes the density
dependence of the interdiffusion coefficients and thus at high density $D_{\rm
AB}$ is proportional to $D_{\rm B}$.

We show that the approximative proportionality $D_{\rm AB}\propto D_{\rm B}$
is associated with strong finite-size effects of the selfdiffusion coefficient
of the slow large species, $D_{\rm A}$. These finite-size effects are due to
the relation of $D_{\rm AB}$ to the diffusion coefficient of the centre of
mass of species $\alpha$ (with $\alpha={\rm A,B}$), $D_{\rm cm}^{(\alpha)}$.
Note that $D_{\rm cm}^{\rm(A)}=D_{\rm cm}^{\rm(B)}$ holds because the total
system’s centre of mass is fixed. As we shall see below, $D_{\rm
cm}^{\rm(A)}\propto D_{\rm AB}/N$, with $N$ the total number of particles in
the system. Thus, the self-diffusion coefficient of the A species, $D_{\rm
A}$, has a finite-size correction $\propto D_{\rm AB}/N\propto D_{\rm B}/N$
that may be the dominant contribution to $D_{\rm A}$ for $\rho^{\rm
A}\gtrsim\rho_{c}^{\rm A}$ and small system sizes. Only if one corrects the
selfdiffusion coefficient $D_{\rm A}$ by computing it relative to the center
of mass of the A species, one can extract the true value of $D_{\rm A}$
without the $1/N$ correction. As we argue below, similar features could be
observed in any glassforming system with strong dynamic heterogeneities. In
such systems, there can be clusters of slow particles with a relatively fast
center-of-mass motion on time scales where no particle rearrangements inside
the cluster occur. Therefore, our study reveals a common feature in the
dynamics of glassforming liquids.

The rest of the paper is organized as follows: In Sec. II, we present the
model of the AB mixture, the details of the simulation and quantities used to
analyze the structure and dynamics of the system. The results of the analysis
of structure and dynamics are given in Sec. III, followed by a summary and
conclusions in Sec. IV.

##  II Model and methods

###  II.1 Interaction potential and details of the simulation

The system that we study horbach2009  is a binary $50-50$ mixture of repulsive
particles, where the diameter of the bigger particles (species A) is sampled
from a uniform distribution, i.e. $d_{\rm A}\in[0.85,1.15]$, while the
diameter of the smaller particles (species B) is $d_{\rm B}=0.35$ (see Fig.
1). The average size ratio of A and B particles is $\langle d_{\rm
A}\rangle/d_{\rm B}\approx 2.85$ where $\langle d_{\rm A}\rangle\approx 1$. A
pair of particles $\\{\alpha,\beta\\}$ (with $\alpha={\rm A,B}$ and
$\beta={\rm A,B}$), separated by a distance $r$, interacts via a Weeks-
Chandler-Andersen (WCA) potential weeks1971 , i.e. a Lennard-Jones potential
that is cut off at its minimum and shifted to zero. To further smoothen the
WCA potential, we also add a term that provides the continuity of its
derivative. Thus, the potential is defined by

$\displaystyle\textrm{V}_{\alpha\beta}(r)$ $\displaystyle=$ $\displaystyle
u_{\alpha\beta}(r)-u_{\alpha\beta}(R_{c})-\left(r-R_{c}\right)\left.\frac{du_{\alpha\beta}}{dr}\right|_{r=R_{c}},$
$\displaystyle u_{\alpha\beta}(r)$ $\displaystyle=$ $\displaystyle
4\epsilon_{\alpha\beta}\left[\left(\sigma_{\alpha\beta}/r\right)^{12}-\left(\sigma_{\alpha\beta}/r\right)^{6}\right]\>,$
(3)

for $r<R_{c}=2^{1/6}\sigma_{\alpha\beta}$, else
$\textrm{V}_{\alpha\beta}(r)=0$, with $\alpha,\beta$ as particle index. Here,
$\sigma_{\alpha\beta}=(d_{\alpha}+d_{\beta})/2$ and
$\epsilon_{\alpha\beta}=\epsilon=1.0$ if both interacting particles are of the
same type, else $\sigma_{\alpha\beta}=(1.00+0.35)/2$ and
$\epsilon_{\alpha\beta}=0.1$. All particles have the same mass $m_{\rm
A}=m_{\rm B}=m=1$. In the following, length, energy, and time are measured in
units of $\langle d_{\rm A}\rangle$, $\epsilon$, and $\tau_{\rm WCA}=[m\langle
d_{\rm A}\rangle^{2}/\epsilon]^{1/2}$, respectively.

Figure 1: Snapshot of the system at the density $\rho=2.5$. A and B particles
are represented by green and black spheres, respectively.

Using LAMMPS plimpton1995 , extensive molecular dynamics (MD) simulations are
performed for systems with $N=1000$, 2000, and 4000 particles, placed in a
three-dimensional cubic box with periodic boundary conditions. The equations
of motion are integrated using the velocity form of the Verlet algorithm
allen2017  with a time step $dt=0.00075\,\tau_{\rm WCA}$. The simulations are
done for different number densities in the range $2.1\leq\rho\leq 3.5$ at the
temperature $T=2/3$. For each density, 30 independent samples are simulated.
During the equilibration of the samples, the temperature is kept constant
using a dissipative particle dynamics (DPD) thermostat soddemann2003  where
the damping coefficient is set to 1.0. For the high densities, $\rho\geq
2.25$, the MD simulations are combined with the swap Monte Carlo (SMC)
algorithm grigera2001 ; berthier2019  to obtain well-annealed samples. In a
trial SMC move, one randomly selects a pair of particles, exchanges their
diameters, and accepts or rejects this move according to a Metropolis
criterion. In our scheme, only the diameters of the A particles are swapped.
Every 100 MD steps, $N_{\rm A}$ trial SMC moves are done. The longest
equilibration runs with the hybrid MD-SMC method were over $6\times 10^{8}$
time steps. After the equilibration, the thermostat as well as the SMC are
switched off, performing the production runs in the microcanonical ensemble.

A snapshot of the system at the density $\rho=2.5$ is shown in Fig. 1. From
this snapshot, one can infer that the large A particles form a close-packed
structure while the small B particles can explore the free volume between the
A particles.

###  II.2 Structural and dynamic properties

In this section, we define the correlation functions and transport
coefficients that we use to analyse the simulation results for our AB mixture.
A central static correlation function for our analysis is the concentration-
concentration structure factor $S_{cc}(q)$. In the limit $q\to 0$, this
function is related to the thermodynamic factor $\Phi$ in Eq. (1). After
having introduced $S_{cc}(q)$ and its relation to $\Phi$, we show how the
selfdiffusion as well as the interdiffusion coefficients can be computed via
Einstein relations, i.e. via long-time limits of mean-squared displacements.

We consider an AB mixture that contains a total number of $N=N_{\rm A}+N_{\rm
B}$ particles. Thus, the concentration of A and B particles is given by
$x_{\rm A}=N_{\rm A}/N$ and $x_{\rm B}=N_{\rm B}/N$, respectively. The local
number density in reciprocal space for particles of type $\alpha$ can be
defined as follows hansen1986 :

$\rho_{\alpha}(\vec{q})=\sum_{j=1}^{N_{\alpha}}\exp\left(i\vec{q}\cdot\vec{r}_{j}\right)$
(4)

with $\vec{q}$ the wavevector and $\vec{r}_{j}$ the position of the $j$’th
particle of type $\alpha$. The autocorrelation functions of the density
variables, as defined by Eq. (4), are the partial structure factors hansen1986
,

$\displaystyle S_{\alpha\beta}(q)$ $\displaystyle=$
$\displaystyle\frac{1}{N}\left\langle\rho_{\alpha}(\vec{q})\rho_{\beta}(-\vec{q})\right\rangle$
(5) $\displaystyle=$
$\displaystyle\frac{1}{N}\sum_{j=1}^{N_{\alpha}}\sum_{k=1}^{N_{\beta}}\left\langle\exp\left[-i\vec{q}\cdot\left(\vec{r}_{j}-\vec{r}_{k}\right)\right]\right\rangle,$

where $\langle.\rangle$ indicates an ensemble average. Note that we assume in
Eq. (5) that the system is isotropic and thus the partial structure factors
only depend on the magnitude of the wavevector, $q$.

From the densities (4), we can introduce local concentration fluctuation
variables hansen1986  as

$c_{\alpha}(\vec{q})=\rho_{\alpha}(\vec{q})-x_{\alpha}\left(\rho_{\rm
A}(\vec{q})+\rho_{\rm B}(\vec{q})\right)\,,$ (6)

describing the local deviation from a homogeneous distribution of particles of
type $\alpha$. Since $c_{\rm A}+c_{\rm B}=0$, the concentration variables are
not independent of each other and it suffices to define one concentration-
concentration structure factor for the binary mixture as

$S_{cc}(q)=\frac{1}{N}\left\langle c_{\rm A}(\vec{q})c_{\rm
A}(-\vec{q})\right\rangle\,.$ (7)

This function can be also expressed as a linear combination of the partial
structure factors,

$S_{cc}(q)=x_{B}^{2}S_{\rm AA}(q)+x_{\rm A}^{2}S_{\rm BB}(q)-2x_{\rm A}x_{\rm
B}S_{\rm AB}(q).$ (8)

In the limit $q\to\infty$, $S_{cc}(q)$ approaches $x_{\rm A}x_{\rm B}$,
corresponding to the concentration-concentration structure factor of an ideal
binary mixture. In the limit $q\to 0$, $S_{cc}(q)$ is related to the
concentration susceptibility and thus to the second derivative of the Gibbs
free energy $G$ with respect to $x_{\rm A}$ and $x_{\rm B}$ via

$\Phi=\frac{x_{\rm A}x_{\rm B}}{k_{\rm B}T}\frac{\partial^{2}G}{\partial
x_{\rm A}x_{\rm B}}=\frac{x_{\rm A}x_{\rm B}}{S(q=0)}\,,$ (9)

with $k_{\rm B}$ the Boltzmann constant. Below, we use Eq. (9) to compute the
thermodynamic factor from an extrapolation of $S_{cc}(q)$ to $q=0$.

Now we introduce quantities that characterize dynamic properties of the AB
mixture. At a single-particle level, we consider the incoherent intermediate
scattering function $F_{\rm s}^{\alpha}(q,t)$ of a tagged particle of type
$\alpha$. This is the correlation function of the time-displaced one-particle
density and defined by hansen1986

$F_{\rm
s}^{\alpha}(q,t)=\frac{1}{N_{\alpha}}\sum_{j=1}^{N_{\alpha}}\left\langle\exp\left[-i\vec{q}\cdot\left(\vec{r}_{j}(t)-\vec{r}_{j}(0)\right)\right]\right\rangle.$
(10)

Of special interest for the analysis of the dynamics of glassforming liquids
is the decay of $F_{\rm s}^{\alpha}(q,t)$ as a function of time around values
of $q$ corresponding to the location of the first peak of the static structure
factor, $q_{\rm max}$. This is due to the fact that slowing down of the
structural relaxation of the glassforming liquid is associated with the cage
effect binder2011  and the typical size of a cage is of the order of
$2\pi/q_{\rm max}$. In our case, the A particles exhibit a typical glassy
dynamics and the first peak of $S_{\rm AA}(q)$ is at $q_{\rm max}\approx 6$.
Thus, below we consider $F_{\rm s}^{\alpha}(q,t)$ at this value of $q$.

We determine the selfdiffusion coefficient of a tagged particle of type
$\alpha$ from the mean-squared displacement (MSD), defined by

$\left\langle\delta
r^{2}_{\alpha}(t)\right\rangle=\frac{1}{N_{\alpha}}\sum_{j=1}^{N_{\alpha}}\left\langle\left(\vec{r}_{j}(t)-\vec{r}_{j}(0)\right)^{2}\right\rangle,$
(11)

The corresponding selfdiffusion coefficient $D_{\alpha}$ is obtained from the
long-time limit of $\left\langle\delta r^{2}_{\alpha}(t)\right\rangle$ via the
Einstein relation hansen1986 ; binder2011

$D_{\alpha}=\lim_{t\to\infty}\frac{\left\langle\delta\tilde{r}^{2}_{\alpha}(t)\right\rangle}{6t}.$
(12)

We also compute modified versions of the incoherent intermediate scattering
function and the MSD where we replace the coordinates of the particles
$\vec{r}_{j}(t)$ in Eqs. (10) and (11) by

$\vec{r}^{\;\prime}_{j}(t)=\vec{r}_{j}(t)-\vec{R}_{\rm A}(t)\,.$ (13)

In this equation, $\vec{R}_{\rm A}(t)$ is the center-of-mass coordinate of the
A particles at time $t$,

$\vec{R}_{\rm A}(t)=\frac{1}{N_{\rm A}}\sum_{j=1}^{N_{\rm
A}}\vec{r}_{j}(t)\,,$ (14)

with $\vec{r}_{j}$ the position of the $j$’th particle of type $A$. The
purpose of calculating $F_{\rm s}^{\alpha}(q,t)$ and $\left\langle\delta
r^{2}_{\alpha}(t)\right\rangle$ with the center-of-mass-corrected coordinates
(13) will become clear below.

The interdiffusion coefficient $D_{\rm AB}$ can be also calculated from an
Einstein relation, i.e. from the long-time limit of a MSD. In this case, the
MSD of the variable $\vec{R}_{\rm A}(t)$ has to be considered,

$\left\langle\delta R^{2}(t)\right\rangle_{\rm
cm}=\left\langle\left[\vec{R}_{\rm A}(t)-\vec{R}_{\rm
A}(0)\right]^{2}\right\rangle\,.$ (15)

Then, the interdiffusion coefficient is given by horbach2007

$D_{\rm AB}=\Phi L\,,$ (16)

where $\Phi$ is the thermodynamic factor, as defined by Eq. (9), and the
Onsager coefficient

$L=\left[1+\frac{x_{\rm A}}{x_{\rm B}}\right]^{2}Nx_{\rm A}x_{\rm
B}\lim_{t\to\infty}\frac{\left\langle\delta R^{2}(t)\right\rangle_{\rm
cm}}{6t}\,,$ (17)

describes the kinetic part of $D_{\rm AB}$.

One can also define a center-of-mass diffusion coefficient for each species
which, for species A, is given by

$D_{\rm A}^{\rm
cm}=\lim_{t\to\infty}\frac{\left\langle\delta{R^{2}(t)}\right\rangle_{\rm
cm}}{6t}\,.$ (18)

For the mixture that is being studied, where $x_{\rm A}=x_{\rm B}=1/2$, we
therefore obtain $D_{\rm AB}=\Phi ND_{\rm A}^{\rm cm}$ and $L=ND_{\rm A}^{\rm
cm}$.

Formally, the interdiffusion coefficient $D_{\rm AB}$ can be written as a
linear combination of the selfdiffusion coefficients,

$D_{\rm AB}=\Phi S\left(x_{\rm A}D_{\rm B}+x_{\rm B}D_{\rm A}\right)\,,$ (19)

where the “Manning factor” manning1961  $S$ contains all the cross
correlations horbach2007  that contribute to $L$. The value $S=1$ implies
vanishing cross correlations akcasu1997 ; horbach2007 . In this case, the
Darken equation (1) holds.

##  III Results

###  III.1 Structural and thermodynamic properties

Figure 2: Pressure $P$ as a function of density. The inset zooms in the region
around $\rho=2.296$ (this density is indicated by an arrow).

With the aid of the SMC technique, we are able to obtain well-annealed samples
at very high densities that are above the critical MCT density
$\rho_{c}\approx 2.23$ horbach2009 . However, at such high densities, one may
expect that at least for the A particles, the thermodynamic equilibrium is an
ordered crystalline phase. Although the large polydispersity of the A
particles in our model suppresses crystallization to some extent, the use of
the SMC technique tends to also accelerate the formation of crystalline
clusters and therefore we check especially the high-density samples whether
they are purely amorphous structures and thus free of any crystalline
clusters. To this end, we measure static structure factors and analyse the
samples in terms of local bond order parameters. Furthermore, we determine the
thermodynamic factor $\Phi$ from the concentration-concentration structure
factor $S_{cc}(q)$.

The pressure $P$ increases monotonously with increasing density $\rho$ (Fig.
2). However, the lines connecting the points between $\rho=2.25$ and
$\rho=2.296$ in $P(\rho)$ indicate a slight change of the slope (cf. the inset
of Fig. 2). This could be due to a liquid-solid coexistence, occurring around
these densities. And indeed our analysis for the density $\rho=2.296$ (see
below) indicates the occurrence of crystalline clusters in some of the samples
at that density.

Figure 3: Top panels: Partial structure factors $S_{\rm AA}(q)$, $S_{\rm
BB}(q)$, and $S_{\rm AB}(q)$ at different densities. Emerging Bragg peaks are
indicated by circles. The inset of the plot of $S_{\rm AA}(q)$ shows the total
structure factor $S_{nn}(q)$ at the density $\rho=2.296$. Bottom panel:
Concentration-concentration structure factor $S_{cc}(q)$ at different
densities. The dashed lines indicate the extrapolation to $q=0$ via a fit
function (see text). The thermodynamic factor $\Phi$, as obtained from the
fitted $S(q=0)$, is shown in the inset as a function of density.

To quantify the structural changes in our AB mixture with increasing density,
we now consider the partial structure factors $S_{\alpha\beta}(q)$. The top
panels of Fig. 3 show these functions for different values of $\rho$. With
increasing $\rho$, the first peak in both $S_{\rm AA}(q)$ and $S_{\rm BB}(q)$
shifts to larger $q$ as the inter-particle separation decreases. In between,
at the density $\rho=2.296$, $S_{\rm AA}(q)$ shows signatures of possible
formation of local crystallites, with discrete spikes being clearly visible.
Although less prominent, this is also reflected in the cross correlation,
$S_{\rm AB}(q)$, as well as in the total structure factor $S_{nn}(q)=S_{\rm
AA}(q)+S_{\rm BB}(q)+2S_{\rm AB}(q)$ that is shown in the inset of the plot of
$S_{\rm AA}(q)$. However, for densities higher than $\rho=2.296$ there is no
sign of any Bragg peaks, suggesting that the samples at these high densities
is purely amorphous. In fact, this is confirmed by our analysis of the samples
in terms of local bond order parameters (see below).

In the bottom panel of Fig. 3, we show how $S_{cc}(q)$ varies with increasing
density. We extract $S_{\rm cc}(0)$ by extrapolating $S_{cc}(q)$ to $q=0$,
using the fit function $f(q)=S_{cc}(0)[1-A\,q^{2}+B\,q^{4}]$ (with
$S_{cc}(0)$,$A$, and $B$ being fit parameters). For the fits (dashed lines),
we have only taken into account the data for $q\leq 2.0$. Then, from
$S_{cc}(0)$ we compute the thermodynamic factor $\Phi$ via Eq. (9). The
variation of $\Phi$ is shown in the inset. We observe that $\Phi$ decreases
from a value of about 1.0 at $\rho=2.0$ to a value of about 0.5 at the highest
considered density, $\rho=3.5$. This rather weak variation of $\Phi$ over a
broad range of densities also implies that the thermodynamic factor does not
strongly affect the density dependence of the interdiffusion coefficient,
which we discuss further below.

Figure 4: $\bar{q}_{4}-\bar{q}_{6}$ plot for system at different marked
densities (calculation with cutoff). The closed line in the plot for
$\rho=2.296$ indicates the values of $\bar{q}_{4}-\bar{q}_{6}$ pairs that
correspond to crystalline BCC clusters (see text).

Our analysis to check for local crystalline order is based on the local bond
order parameters that have been proposed by Steinhardt et al. steinhardt1983 .
For such analysis, at first we identify the nearest neighbours for each
particle within the cutoff $1.5$ ($N_{b}(i)$ is the number of such neighbours
for each particle $i$). After this, for each particle, a complex local
orientational order vector $q_{lm}(i)$ with $(2l+1)$ components is constructed
using the following definition:

$q_{lm}(i)=\frac{1}{N_{b}(i)}\sum_{j=1}^{N_{b}(i)}Y_{lm}(\vec{r}_{ij}).$ (20)

Here, $Y_{lm}(\vec{r}_{ij})$ is the spherical harmonic of degree $l$ and order
$m$. As per definition, $l$ is always a non-negative integer and $m$ can take
the integral values from $m=-l$ to $m=l$ for a given value of $l$. Also,
$\vec{r}_{ij}$ is the vector from the particle $i$ to particle $j$. Now, as
suggested in Ref. dellago2008 , the locally averaged quantity
$\bar{q}_{lm}(i)$ is calculated for each particle $i$, using

$\bar{q}_{lm}(i)=\frac{1}{\tilde{N}_{b}(i)}\sum_{k=0}^{\tilde{N}_{b}(i)}q_{lm}(k),$
(21)

where averaging (summation from $k=0$ to $\tilde{N}_{b}(i)$) has been done
over the neighboring particles of $i$ (as above defined, based on the cutoff)
and particle $i$ itself. Such averaging procedure takes into account even the
information of the structure beyond the cutoff and has shown better
identification of crystal structures in simulations. Using these averaged form
of local bond order vector components a norm $\bar{q}_{l}(i)$ is defined for
each particle:

$\bar{q}_{l}(i)=\sqrt{\frac{4\pi}{2l+1}\sum_{m=-l}^{l}|\bar{q}_{lm}(i)|^{2}}.$
(22)

Depending on the definition of the spherical harmonics used, the factor
$\frac{4\pi}{2l+1}$ before summation in the above equation, is used or not
used. This quantity $\bar{q}_{l}(i)$ is sensitive to different crystal
structure systems depending on the choice of $l$. Specially, bond
orientational order parameter for $l=4$ and $l=6$ have been used to identify
structures similar to cubic and hexagonal systems. We have calculated the
correlation between $\bar{q}_{4}$ and $\bar{q}_{6}$ for the A particles in our
system.

The measured values of $\bar{q}_{4}$ and $\bar{q}_{6}$, for each particle
across several representative configurations, can be visualized in the form of
a scatter plot in the $\bar{q}_{4}$-$\bar{q}_{6}$ plane, following Ref.
dellago2008 , to check if local environments exhibit any formation of ordered
structures. This is shown in Fig. 4, for a few densities across the range that
we have studied, using configurations sampled from $30$ independent
trajectories. First thing to notice is that, across all the densities, the
values of $\bar{q}_{4}$ are small ($<0.06$) and in fact shrink with increasing
density. Similarly for $\bar{q}_{6}$, the numbers are in the liquid-like
regime, except for $\rho=2.296$ where some pockets of BCC-like structures seem
to be visible for some clusters within some trajectories, if we compare with
observations reported in Ref. dellago2008 .

To summarize, the structural changes with increasing density are on expected
lines. There is some hint of occurrence of a small number of local
crystallites, in the vicinity of the mode coupling density. However, when
supercooled to higher densities, the system remains disordered, even when
equilibriated via SMC.

###  III.2 Diffusion dynamics

We will now discuss the dynamic properties of the mixture. As mentioned
earlier, for the model system that we are studying, the larger A species are
known to exhibit a mode coupling transition around $\rho_{c}\approx 2.23$
horbach2009 . For $\rho\gtrsim\rho_{c}$, the A particles are essentially in a
frozen-in configuration on the diffusive time scale of the B particles.
However, the diffusivity of the B particles also continues to decrease with
increasing density such that they are expected to be in an arrested state at
very large densities.

Figure 5: (a) Center-of-mass MSD, $\left\langle\delta
R^{2}(t)\right\rangle_{\rm cm}$, for different densities. (b) Single-particle
MSDs of A species and (c) of B species for the same densities. In b) and c),
solid and dotted lines represent calculations with and without the center-of-
mass correction, respectively (see text). In all sub-plots, the solid lines
with dots are straight lines $\propto t$ to indicate the diffusional regime of
the MSDs at long times.

Figure 6: Top panel: Selfdiffusion coefficients $D_{\rm A}$ and $D_{\rm B}$,
the center-of-mass diffusion coefficient $D_{\rm A}^{\rm cm}$ and the
interdiffusion coeffient $D_{\rm AB}$ as a function of density $\rho$. Also
shown is the rectified self-diffusion coefficient $\tilde{D}_{\rm A}$ (see
text). The vertical dashed line marks the previously estimated critical MCT
density for $D_{\rm A}$ at $\rho_{c}=2.23$. Bottom panel: Manning factor $S$
and the product $S\Phi$ as a function of density.

In the following, our main objective is to study the density dependence of the
interdiffusion coefficient. To this end, one has to monitor the trajectory of
the center of mass of the A species and the corresponding time evolution of
$\langle\delta R^{2}(t)\rangle_{\rm cm}$, as defined by Eq. (15), needs to be
computed. This MSD is shown in Fig. 5a for different densities. For the same
densities, Figs. 5b and 5c display the single-particle MSDs, $\langle\delta
r^{2}_{\alpha}(t)\rangle$ (see Eq. (11)), of the A and B particles,
respectively. For the A species, $\langle\delta r^{2}_{\rm A}(t)\rangle$ shows
the behaviour that is expected for a typical glassforming liquid. There is a
short-time ballistic regime ($\propto t^{2}$) and a long-time diffusive regime
($\propto t$). In between these two regimes, there is a plateau-like region
that becomes more pronounced and broader with increasing density. The latter
regime is due to the intermediate caging of the particles. Importantly, we
note that even for $\rho>\rho_{c}$, the MSD of A particles still displays
diffusive dynamics at long times. Below, we show that this feature is a
finite-size effect, i.e. the observed diffusive regime for $\rho>\rho_{c}$
shifts to longer time scales with increasing system size. For the case of B
species, the shape of the MSD, i.e. of $\langle\delta r^{2}_{\rm
B}(t)\rangle$, is very different from that of the A particles, with no
intermediate plateau at all. We will also discuss this, below. Now, if we look
at the MSD curves for the center-of-mass motion, it is evident that they bear
close resemblance with that of the B species, the reason for which will be
evident once we analyse the diffusion coefficients.

The long-time dynamical behaviour is well characterized by measuring the
respective diffusion coefficients, measured from the corresponding MSD data,
as defined in Eqs. (12), (16), and (18). The density variation of the measured
single-particle diffusion coefficients, $D_{\rm A}$ and $D_{\rm B}$, the
center-of-mass diffusion coefficient $D_{\rm A}^{\rm cm}$ as well as the
interdiffusion coeffient $D_{\rm AB}$ are shown in Fig. 6. First, we note that
the diffusion coefficient of the smaller B particles remains finite beyond
$\rho_{c}$ (indicated via the dashed vertical line), consistent with previous
study horbach2009 . Now, if we look at the data for $D_{\rm A}$, we observe
that it decreases as it approaches $\rho_{c}$, as reported earlier, and beyond
that there seems to be a separate branch of weaker decrease with density.
Interestingly, for $\rho>\rho_{c}$, the measured center-of-mass diffusion
coefficient, $D_{\rm A}^{\rm cm}$, shows exactly the same density dependence.
As already noted in the previous paragraph, the measured interdiffusion
coefficient has a density dependence which resembles that of the diffusivity
of B species, over the entire density range.

We now analyse the motion of A particles for $\rho>\rho_{c}$. Since the
density dependence of $D_{\rm A}$, in this density regime, matches that of
$D_{\rm A}^{\rm cm}$, it implies that the observed motion of the individual A
particles is actually due to the motion of the center of mass of the A
species. To disentangle that, we compute the MSD of the A particles by
shifting to the frame of reference of the population’s center of mass, see Eq.
(13). The redefined MSDs are plotted in the top panel of Fig. 5, using dashed
lines. We observe that following this rectification, there is a significant
change in the MSD for $\rho>\rho_{c}$. The rectified MSDs, in that regime,
exhibit a prolonged plateau over the time window of our observation implying
that in the frame of reference of the center of mass, the particles are
essentially caged and there is hardly any cage-breaking, especially for
$\rho\geq 2.5$. Therefore, we can conclude that only due to the motion of the
center of mass, deviations from the plateau are exhibited in the bare MSDs at
these densities. At small enough density ($\rho=2.1$), this rectification is
not observed, and there is a mild rectification for $\rho=2.205$. The
rectified selfdiffusion coefficient of the A particles is also shown in Fig.
6, from which it is evident that the rectified $D_{\rm A}$, labelled as
$\tilde{D}_{\rm A}$, sharply decreases around $\rho_{c}$.

Figure 7: Top panel: Center-of-mass MSD for A species, scaled with system size
$N$, at density $\rho=2.5$. Bottom panel: Single particle MSD of the A species
at the same density, for the different system sizes. Solid and dotted lines
represent calculations without and with center-of-mass correction,
respectively.

Having clarified the actual density dependence of $D_{\rm A}$, we can now use
Eq. (19) to analyse the behaviour of the interdiffusion coefficient. If the
Manning factor is $S=1.0$, Eq. (19) reduces to the Darken equation (1). Then,
when $D_{\rm A}=0$, $D_{\rm AB}\sim D_{\rm B}$, which explains the observed
behaviour of the interdiffusion coefficient in the regime $\rho\gg\rho_{c}$.
Physically, even though there is a dynamical arrest of the A species, because
of momentum conservation, the mobility of B species leads to the center-of-
mass motion of the A species and the consequent finite interdiffusion
coefficient. We also note that even at lower densities, $D_{\rm B}>D_{\rm A}$,
and thus there too, the behaviour of $D_{\rm AB}$ is also dependent on the
diffusive motion of the smaller but faster B particles.

Furthermore, we test the validity of the Darken approximation for the mixture
that we are studying. To do that, we consider two different quantities, viz
$S=D_{\rm AB}/(x_{\rm A}D_{\rm B}+x_{\rm B}D_{\rm A})$ and $S\Phi$, which are
plotted in the bottom panel of Fig. 5. If the Darken approximation works, then
$S$ should be around 1. As we can infer from Fig. 5, over the whole range of
densities, both $S$ and $S\Phi$ deviate from unity by less than a factor of 2,
while the interdiffusion coefficient decreases by about two orders of
magnitude. So neither the cross correlations nor the thermodynamic factor
strongly affect the density dependence of $D_{\rm AB}$ which is therefore
essentially given by the linear combination of selfdiffusion coefficients,
i.e. $D_{\rm AB}\approx 0.5(D_{\rm A}+D_{\rm B}$).

Figure 8: Self-intermediate scattering function, $F_{s}^{\alpha}(q,t)$, for
bigger A (top) and smaller B (bottom) particles, measured at $q=6$, for
different densities. In the top panel, solid and dashed lines represent
calculations without and with center-of-mass correction, respectively. The
dotted lines are fits with stretched exponentials (see text). The values of
the stretching exponent are $\beta=0.73$, 0.68, 0.89, 0.86, and 0.87 for the
densities $\rho=2.1$, 2.205, 2.296, 2.5, and 3.0, respectively.

In the above discussions, we have noted that the MSD curves for B species do
not show any plateau-like feature, i.e. the absence of any signatures of
caging by neighbouring particles. Rather, these curves show anomalous
intermediate sub-diffusive dynamics prior to diffusion. This is similar to
what is observed in the case of interacting particles moving in a quenched
environment of soft obstacles schnyder2018 ; schnyder2015 . For the binary
mixture that we are studying, this is a reasonable scenario considering the
large size ratio (see Fig. 1 for the visualisation). At the large density
beyond $\rho_{c}$ where we are probing the dynamics, the A population is
nearly frozen, and the B particles are diffusing through this quenched
environment with which they interact via some soft interaction. For the case
of interacting particles moving in a quenched environment of soft obstacles,
an avoided localization transition is observed schnyder2018 , and one would
expect a similar situation for the eventual dynamical arrest of the B species
at large enough densities horbach2009 .

Next, we analyse the finite size effects in the measured diffusivities. Since
the interdiffusion coefficient can be written as $D_{\rm AB}=\Phi L$, where
$L=ND_{\rm A}^{\rm cm}$, as discussed in the previous section, $D_{\rm AB}$
has no finite size effects. This is illustrated in the top panel of Fig. 7,
where we plot the center-of-mass MSD of A species scaled with system size, at
density of $\rho=2.5$, for a density larger than $\rho_{c}$. We observe that
the curves for the different system sizes collapse for the scaled MSD,
implying that ${N}D_{\rm A}^{\rm cm}$ is the same and therefore $D_{\rm AB}$
remains unaffected. On the other hand, the bare single particle MSDs show
system size dependence at the same density, as shown in the bottom panel of
Fig. 7. However, if we now measure the MSD of the particles in the frame of
reference of the center of mass of the A particles, see Eq. (13), we observe
that all collapse on the same plateau; see top panel of Fig. 7. This implies,
that the finite size dependence in the single particle dynamics is coming from
the change in centre-of-mass motion. However, as discussed above, this effect
gets scaled out when computing the interdiffusion coefficient.

So far, we have been discussing the anomalous effects in the mean squared
displacement for $\rho>\rho_{c}$ and the rectification that is needed. This is
also the case in other dynamical observables. In the top panel of Fig. 8, we
illustrate this for the incoherent intermediate scattering function
$F_{s}^{\rm A}(q,t)$ at a wavenumber of $q=6$, corresponding to correlations
on nearest-neighbour distances. The bare measurements are shown with solid
lines and the rectified curves as dashed lines, i.e. measurements done in the
frame of reference of the center of mass. Consistent with the MSD data, at low
densities (here $\rho=2.1$), the bare and the rectified $F_{s}^{\rm A}(q,t)$
coincide, while at higher densities the rectified function shows a slower and
a more stretched decay (cf. the curves for $\rho=2.205$). For
$\rho\geq\rho_{c}$, the center-of-mass-corrected functions show extended
plateaus in time, reflecting that in the frame of reference of the center of
mass, the A particles are stuck in their locations.

The bare curves show a relaxation process, which can be fitted with stretched
exponential functions, $f(q,t)=A_{\rm sef}\exp[-(t/\tau_{q})^{\beta}]$ (with
$A_{\rm sef}$, $\tau_{q}$, and $\beta$ being fit parameters). The values of
the exponent $\beta$, as obtained from the fits, are in the range between 0.7
and 0.9 (see caption of Fig. 8). With increasing density, the stretching
exponent $\beta$ tends to increase towards 1.0, i.e. towards a simple
exponential decay. This decay of $F_{s}^{\rm A}(q,t)$ at $q=6.0$ reflects the
center-of-mass motion of the A species. At very long times, however, there is
the crossover to a much slower decay, which is due to the fact that there are
no rearrangements in the configuration of A particles. The bottom panel of
Fig. 8 shows $F_{s}^{\rm B}(q,t)$. Here, for the low densities, the decay
happens very quickly. Only for $\rho>\rho_{c}$, the relaxation timescales
start increasing significantly.

##  IV Summary and conclusions

Our work elucidates several aspects of the diffusion dynamics in an equimolar
glassforming binary soft-sphere mixture with large size ratio. We observe a
very pronounced time scale separation between the motion of the slow A species
and that of the fast B species at high densities. As a consequence, the A
particles show the typical glassy dynamics of a densely packed system, while
the B particles remain mobile at very high densities, $\rho>\rho_{c}$,
exploring the void space in between the A particles. We have seen that the
interdiffusion coefficient $D_{\rm AB}$ can be well approximated by the Darken
equation (1). Since at high density $D_{\rm A}\ll D_{\rm B}$, one therefore
obtains $D_{\rm AB}\approx\Phi x_{\rm A}D_{\rm B}$, i.e. the interdiffusion
coefficient is dominated by the fast species. For our equimolar mixture, we
can also write $D_{\rm AB}=\Phi ND_{\rm A}^{\rm cm}$ and thus at high density
we have $D_{\rm A}^{\rm cm}\approx D_{\rm B}/(2N)$. This implies that even if
the A particles form a frozen-in structure with essentially no rearrangements
of the relative positions of particles, they may follow collectively the
diffusive motion of their center of mass and one finds for the selfdiffusion
coefficient $D_{\rm A}\approx D_{\rm A}^{\rm cm}\propto N^{-1}$. One can, of
course, correct for this finite-size effect by computing one-particle
quantities such as $\langle\delta r^{2}_{\rm A}(t)\rangle$ and $F_{s}^{\rm
A}(q,t)$ from the center-of-mass-corrected coordinates, as defined by Eq.
(13).

The finite-size effects, that we have reported in this work with respect to
the selfdiffusion coefficient of the slow species, are expected to be a
typical feature in glassforming binary mixtures with large size ratio.
Moreover, one may expect similar effects in any glassforming system with
strong dynamic heterogeneities.

###### Acknowledgements.

We thank IMSc HPC facility for providing computational resources for our work.

## References

  * (1) C. Bechinger, F. Sciortino, and P. Ziherl (eds.), Physics of Complex Colloids (IOS Press, Amsterdam, 2013).
  * (2) M. Weiss, in New Models of the Cell Nucleus: Crowding, Entropic Forces, Phase Separation, and Fractals, edited by R. Hancock and K. W. Jeon, International review of cell and molecular biology Vol. 307 (Academic Press, San Diego, CA, 2014), Chap. 11, pp. 383–417.
  * (3) F. Höfling and T. Franosch, Rep. Prog. Phys. 76, 046602 (2013).
  * (4) A. Imhof and J. K. G. Dhont, Phys. Rev. Lett. 75, 1662 (1995).
  * (5) A. Imhof and J. K. G. Dhont, Phys. Rev. E 52, 6344 (1995).
  * (6) R. Kurita and E. R. Weeks, Phys. Rev. E 82, 041402 (2010).
  * (7) T. Blochowicz, S. Schramm, S. Lusceac, M. Vogel, B. Stühn, P. Gutfreund, and B. Frick, Phys. Rev. Lett. 109, 035702 (2012).
  * (8) S. P. Bierwirth, C. Gainaru, and R. Böhmer, J. Chem. Phys. 149, 044509 (2018).
  * (9) T. Sentjabrskaja, E. Zaccarelli, C. De Michele, F. Sciortino, P. Tartaglia, T. Voigtmann, S. U. Egelhaaf, and M. Laurati, Nat. Commun. 7, 11133 (2016).
  * (10) E. Martinez-Sotelo, M. A. Escobedo-Sánchez, and M. Laurati, J. Chem. Phys. 151, 164504 (2019).
  * (11) A. J. Moreno and J. Colmenero, J. Chem. Phys. 125, 164507 (2006).
  * (12) J. Horbach, and T. Voigtmann, Phys. Rev. Lett. 109, 205901 (2009).
  * (13) W.-S. Xu, Z.-Y. Sun, and L.-J. An, J. Chem. Phys. 137, 104509 (2012).
  * (14) W.-S. Xu, Z.-Y. Sun, and L.-J. An, Soft Matter 11, 627 (2015).
  * (15) E. Lázaro-Lázaro, J. A. Perera-Burgos, P. Laermann, T. Sentjabrskaja, G. Pérez-Ángel, M. Laurati, S. U. Egelhaaf, M. Medina-Noyola, T. Voigtmann, R. Castaneda-Priego, and L. F. Elizondo-Aguilera, Phys. Rev. E 99, 042603 (2019).
  * (16) J. Bosse and J. S. Thakur, Phys. Rev. Lett. 59, 998 (1987).
  * (17) J. Bosse and Y. Kaneko, Phys. Rev. Lett. 74, 4023 (1995).
  * (18) T. Voigtmann, Europhys. Lett. 96, 36006 (2011).
  * (19) S. K. Schnyder and J. Horbach, Phys. Rev. Lett. 120, 078001 (2018).
  * (20) J. Kurzidim, D. Coslovich, and G. Kahl, J. Phys.: Condens. Matter 23, 234122 (2011).
  * (21) D. D. Fitts, Non-equilibrium Thermodynamics (MacGraw-Hill, New York, 1962).
  * (22) A. Z. Akcasu, Macromol. Theory Simul. 6, 679 (1997).
  * (23) J. Horbach, S. K. Das, A. Griesche, M.-P. Macht, G. Frohberg, and A. Meyer, Phys. Rev. B 75, 174304 (2007).
  * (24) L. S. Darken, Trans. AIME 180, 430 (1949).
  * (25) G. S. Hartley and J. Crank, Trans. Faraday Soc. 45, 801 (1949).
  * (26) P. Kuhn, J. Horbach, F. Kargl, A. Meyer, and Th. Voigtmann, Phys. Rev. B 90(2), 024309 (2014).
  * (27) R. J. Bearman, J. Chem. Phys. 32, 1308 (1960).
  * (28) F. Brochard and P. G. de Gennes, Europhys. Lett. 1, 221 (1986).
  * (29) H. Sillescu, Makromol. Chem., Rapid Commun. 8, 393 (1987).
  * (30) W. Hess, G. Nägele, and A. Z. Akcasu, J. Polym. Sci., Part B: Polym. Phys. 28, 2233 (1990).
  * (31) A. Z. Akcasu, G. Nägele, and R. Klein, Macromolecules 24, 4408 (1991).
  * (32) A. Latz, Verallgemeinerte konstituierende Gleichungen und Formfaktoren für einfache Glasbildner, Ph.D. thesis, TU München, Germany, 1990.
  * (33) J. D. Weeks, D. Chandler, and H. C. Andersen, J. Chem. Phys. 54, 5237 (1971).
  * (34) S. Plimpton, J. Comp. Phys. 117, 1 (1995).
  * (35) M. P. Allen and D. J. Tildesley, Computer Simulation of Liquids, 2nd ed. (Oxford University Press, Oxford, 2017).
  * (36) T. Soddemann, B. Dünweg, and K. Kremer, Phys. Rev. E 68, 046702 (2003).
  * (37) T. S. Grigera and G. Parisi, Phys. Rev. E 63, 045102 (2001).
  * (38) L. Berthier, E. Flenner, C. J. Fullerton, C. Scalliet, and M. Singh, J. Stat. Mech. 064004 (2019).
  * (39) J.-P. Hansen and I. R. McDonald, Theory of Simple Liquids (Academic Press, London, 1986).
  * (40) K. Binder and W. Kob, Glassy Materials and Disordered Solids: An Introduction to Their Statistical Mechanics, Rev. Ed. (World Scientific, Singapore, 2011).
  * (41) J. R. Manning, Phys. Rev. 124, 470 (1961).
  * (42) P. J. Steinhardt, D. R. Nelson, and M. Ronchetti, Phys. Rev. B 28 (2), 784 (1983).
  * (43) W. Lechner and C. Dellago, J. Chem. Phys. 129, 114707 (2008).
  * (44) S. K. Schnyder, M. Spanner, F. Höfling, T. Franosch, and J. Horbach, Soft Matter 11, 701 (2015).