# CLSA: Contrastive Learning-based Survival Analysis for Popularity Prediction
in
MEC Networks

Zohreh Hajiakhondi-Meybodi, , Arash Mohammadi, , Jamshid Abouei, , and
Konstantinos N. Plataniotis Z. Hajiakhondi-Meybodi is with Electrical and
Computer Engineering (ECE), Concordia University, Montreal, Canada. (E-mail:
z_hajiak@encs.concordia.ca). A. Mohammadi (corresponding author) is with
Concordia Institute of Information Systems Engineering (CIISE), Concordia
University, Montreal, Canada. (P: +1 (514) 848-2424 ext. 2712 F: +1 (514)
848-3171, E-mail: arash.mohammadi@concordia.ca). J. Abouei was with the
Department of Electrical and Computer Engineering, University of Toronto,
Toronto, Canada. He is now with the Department of Electrical Engineering, Yazd
University, Yazd 89195-741, Iran (E-mail: abouei@yazd.ac.ir). K. N.
Plataniotis is with Electrical and Computer Engineering (ECE), University of
Toronto, Toronto, Canada. (E-mail: kostas@ece.utoronto.ca).This Project was
partially supported by Department of National Defence’s Innovation for Defence
Excellence & Security (IDEaS) program, Canada.

###### Abstract

Mobile Edge Caching (MEC) integrated with Deep Neural Networks (DNNs) is an
innovative technology with significant potential for the future generation of
wireless networks, resulting in a considerable reduction in users’ latency.
The MEC network’s effectiveness, however, heavily relies on its capacity to
predict and dynamically update the storage of caching nodes with the most
popular contents. To be effective, a DNN-based popularity prediction model
needs to have the ability to understand the historical request patterns of
content, including their temporal and spatial correlations. Existing state-of-
the-art time-series DNN models capture the latter by simultaneously inputting
the sequential request patterns of multiple contents to the network,
considerably increasing the size of the input sample. This motivates us to
address this challenge by proposing a DNN-based popularity prediction
framework based on the idea of contrasting input samples against each other,
designed for the Unmanned Aerial Vehicle (UAV)-aided MEC networks. Referred to
as the Contrastive Learning-based Survival Analysis (CLSA), the proposed
architecture consists of a self-supervised Contrastive Learning (CL) model,
where the temporal information of sequential requests is learned using a Long
Short Term Memory (LSTM) network as the encoder of the CL architecture.
Followed by a Survival Analysis (SA) network, the output of the proposed CLSA
architecture is probabilities for each content’s future popularity, which are
then sorted in descending order to identify the Top-$K$ popular contents.
Based on the simulation results, the proposed CLSA architecture outperforms
its counterparts across the classification accuracy and cache-hit ratio.

###### Index Terms:

Mobile Edge Caching (MEC), Popularity Prediction, Deep Neural Network (DNN),
Contrastive Learning, and Survival Analysis.

##  I Introduction

Mobile Edge Caching (MEC) is a promising technology empowering the Beyond
Fifth-Generation (5G) wireless networks to cope with the unprecedented
expansion of the world’s mobile data traffic [1]. MEC networks bring
multimedia content closer to the Internet of Things (IoT) devices by offering
storage capacities at the edge of the network, resulting in lower
communication latency. The edge network’s storage capacity, however, is
constricted, making it impractical to keep all contents there. Existing
reactive caching schemes [2, 3] identify the most popular contents based on
the observed request, where the content will be cached after being requested.
Users’ preferences, however, dynamically change over time. As a result,
identifying the required contents before being requested and storing them on
edge devices is a crucial part of the MEC networks.

As a prominent approach, Deep Neural Networks (DNN)-based proactive caching
schemes [4, 5, 6, 7, 8, 9, 10, 11] have been introduced. These models predict
the Top-$K$ popular contents in the upcoming time using the historical request
patterns of contents. To further enhance the performance of these schemes,
several time-series learning models have recently been introduced [12, 13, 14]
to capture the temporal information of historical requests. Additionally,
there is a meaningful correlation between the historical patterns of
popular/unpopular content requests, making it necessary to learn the spatial
correlation, along with capturing the temporal dependencies. The paper aims to
further advance this emerging field.

Literature Review: Generally speaking, popularity prediction models can be
classified into three groups, i.e., $(i)$ Statistical models [15, 16], e.g.,
content-based filtering, item-to-item correlation systems, and collaborative
filtering; $(ii)$ Machine Learning (ML)-based schemes [17], such as Random
Forest (RF) [18], Generalized Linear Model (GLM) [19], and Decision Tree (DT)
[20], and; $(iii)$ DNN-based architectures [12, 13, 14, 4, 5, 6, 7, 8, 9, 10,
11, 21, 22], such as Vision Transformers [12, 14], Transformers [13],
Autoencoders [7], Long Short Term Memory (LSTM) [6, 11], and Convolutional
Neural Network (CNN) [23]. Despite all the advantages of existing statistical
and ML-based techniques, there are several limitations that render them
ineffective for real-time caching solutions. On the one hand, statistical
models would not be time-efficient when there are a huge number of
users/multimedia contents. On the other hand, the sparsity and cold-start
issues that occur when insufficient data is provided about a new mobile
user/multimedia content affect the performance of the statistical and ML
models. These issues, however, can be resolved by using existing DNN-based
models. Last but not least, DNN-based frameworks predict the popularity of
contents utilizing the raw historical request patterns of contents. Therefore,
there is no requirement for feature engineering or pre-processing. These
encourage the researchers to concentrate on DNN-based frameworks for
popularity prediction in MEC networks.

Leveraging the CNN for feature extraction, a content-aware popularity
prediction framework was introduced in [4]. Using statistical clustering
methods, contents were categorized into different groups based on users’
preferences. Similarly, Ndikumana et al. [24] used CNN and Multi-Layer
Perceptron (MLP) to predict users’ features and the probabilities that
contents will be requested. Tang et al. [25] implemented a Deep Reinforcement
Learning (DRL) framework to model users’ requesting behaviors. Sadeghi et al.
[26] introduced an adaptive caching framework in hierarchical networks, where
a two-way interactive influence between the cloud and edge devices is
represented by a DRL model. These models, however, are inapplicable in highly
dynamic practical networks since they make the assumption that the content
popularity will not alter over time.

To tackle the aforementioned issue, the main focus of recent works has been
shifted to developing time-series prediction models [5, 27] to process the
sequential and time-variant historical request patterns of contents. For
instance, LSTM has been widely deployed in recent works, including but not
limited to [6, 11, 28] for predicting the number of users’ requests in the
upcoming time. Although LSTM is capable of learning long-term dependencies,
the correlation between the current sample and the earlier ones degrades over
time. Moreover, they are unable to capture the spatial correlation between
multiple contents. Another time-series learning model is Transformer
architectures [29], where there is no need to analyze sequential data in the
same order. In our previous works [12, 14], we deployed Vision Transformer
(ViT) to predict the Top-$K$ popular contents with high accuracy. To
simultaneously predict the popularity of multiple contents and capture the
spatial correlation between different contents, $2$D images were created to be
used as the input of the ViT architecture, where each column of the image was
related to the requests pattern of one content. In another work [13], we
deployed a multi-channel Transformer architecture, where the historical
request pattern of each content is given to a channel. Consequently, the
spatial correlation of multiple contents is captured. Although providing
higher accuracies, the complexity of the learning models in [12, 14, 13]
exponentially increases by increasing the number of contents (i.e., the number
of columns of the input image in the ViT model [12, 14] and the number of
channels in the Transformer architecture [13]). Another issue we encountered
with these models [12, 14, 13] is that the input samples lack the contextual
information of users such as age, and gender, making it ineffective to
precisely understand distinct users’ interests. To tackle these issues, we
focus on the incorporation of Contrastive Learning (CL) [30], which is a self-
supervised learning paradigm that has been widely used (recently) in the
computer vision domain. Intuitively speaking, by capitalizing on a unique
characteristic of the CL paradigm, i.e., learning common attributes and
differences in input samples, there is no need to simultaneously input the
request patterns of all contents to capture the spatial correlations.

Survival Analysis (SA) [31] is another line of research that has been designed
to learn the relationship between the contextual information of different
users requesting the same content, and the distribution of the first time that
the content will be requested again in the future. While SA models are widely
used in medical studies, they are also applied in recommender systems [32] to
predict users’ preferences in the future. For example, Desirena et al. [33]
introduced a survival neural network to better interact with customers by
modeling non-linear user preferences. To the best of our knowledge, there has
been no study conducted on the joint use of CL and SA models for popularity
prediction in MEC networks. This necessitates an urgent quest to develop and
design a new and innovative SA-based popularity prediction architecture, which
is the focus of this paper.

Contributions: Motivated by the above discussion, we introduce the Contrastive
Learning-based Survival Analysis (CLSA) caching framework with the application
to the MEC networks. The proposed CLSA architecture is a time-series
popularity prediction framework learning the real-time content placement from
the historical contextual information of users requesting a content. The
proposed CLSA architecture consists of Reconstruction Network (RN), CL, and SA
blocks, with the following characteristics:

  * •

To provide a better understanding of historical patterns of content requests,
longitudinal measurements are incorporated into a SA model. The output of the
SA model is the probability that one content gets popular within a time
window, which is then sorted in descending order to identify the Top-$K$
popular contents. Moreover, a self-supervised CL network is used, focusing on
the differences between the request patterns of distinct contents.

  * •

Within the CL network, a shared encoder (developed based on the LSTM model) is
used to learn a meaningful representation of historical patterns of distinct
contents to capture the temporal information as well. Consequently, it is
unnecessary to input multiple historical patterns of contents simultaneously
into the CLSA architecture to capture the spatial content correlations. This
results in a significant reduction in the input size of the learning model.
Since a CL network needs positive and negative samples, two augmentation
methods are employed to create positive samples. Finally, a decoder is
employed in the RN block to recreate the original input sample from the
augmented one.

The effectiveness of the proposed CLSA framework is evaluated through
comprehensive studies on the real-trace multimedia request patterns, focusing
on classification accuracy and cache-hit ratio. Simulation results corroborate
that the CLSA framework outperforms its counterparts in all the above-
mentioned aspects.

The remainder of the paper is organized as follows: In Section II, the system
model and the dataset are described. Section III presents the proposed CLSA
architecture. Simulation results are presented in Section IV. Finally, Section
V concludes the paper.

Figure 1: A typical structure of the UAV-aided cellular network.

##  II System Model and Problem Description

We consider an Unmanned Aerial Vehicle (UAV)-aided MEC network, as shown in
Fig. 1, which includes several UAVs as aerial caching nodes and Femto Access
Points (FAPs) as terrestrial infrastructure. We define the set of UAVs as
$u_{s}$, where $s\in\\{1,\ldots,N_{u}\\}$, and the set of FAPs as $f_{i}$,
with $i\in\\{1,\ldots,N_{f}\\}$, serving requests of users, denoted by
$u_{l}$, where $l\in\\{1,\ldots,U\\}$. It is assumed that users request
different multimedia contents, denoted by $c_{m}$, where
$m\in\\{1,\ldots,M\\}$, belonging to a content library with the size of $M$.
Moreover, UAVs and FAPs have a limited storage capacity. Following the common
approach in the literature [34], it is assumed that all contents are the same
size, and users only request one content at a time. In addition, FAPs and
users are randomly distributed through the network with a Poisson Point
Process (PPP) and Gaussian mixture distributions [35], respectively, and UAVs’
locations are obtained as the result of the $K$-Means clustering model [36].
It is assumed that UAVs continue to hover over their locations as they fulfill
a request for data delivery [34]. The link quality and topology of UAVs are
controlled by a Software Defined Network (SDN) controller, which also manages
the aerial and terrestrial links [34].

###  II-A Dataset

MovieLens dataset [37], as one of the most well-known movie recommendation
services that provide users’ contextual information, is used in this study to
evaluate the proposed CLSA architecture. MovieLens was generated on October
$17$, $2016$, including the movie rates that $943$ users gave to $1682$ movies
between September $19$, $1997$, and April $22$, $1998$. As shown in Fig. 2(a),
u.user document provided in the MovieLens dataset, contains users’ contextual
information, including gender, age, occupation, and ZIP code. The ZIP codes
are converted to latitude and longitude coordinates to extract users’
locations during their requests. Another document provided in this dataset is
u.data, including user ID, item ID (content ID), the rate the user gave to the
corresponding content, and the timestamp that the user watched and rated the
content.

As shown in Fig. 2, the following steps are performed to adopt the Movielens
dataset with the CLSA architecture: $(i)$ The u.data and u.user documents are
concatenated on “user ID” (see Fig. 2(b)) and the ZIP code is dropped from the
concatenated dataset; $(ii)$ The concatenated dataset in Fig. 2(b) is sorted
by “item ID”, and “timestamp”; $(iii)$ The “user ID” and “item ID” are dropped
from the concatenated dataset, since these columns are not informative; $(iv)$
The categorical features including gender, age, and occupation are encoded
using a one-hot encoder; $(v)$ The timestamp is discretized with a resolution
of one day since the storage of edge devices should be updated during the off-
peak times (i.e., midnight) [38]. Then, the discretized timestamp, named
“day”, is replaced with the timestamp column; and, $(vi)$ Finally, a column
named “label” is added to the concatenated dataset, indicating the content
popularity (popular or unpopular), which will be described shortly.

We define an observational window for each content $c_{m}$, for ($1\leq m\leq
M$) having the length of $T_{\tau}^{m}$ at time $\tau$, where the request
pattern of contents within this range are studied to predict their popularity
in the future, i.e., the study window, denoted by $T_{s}$. Following the
Reference [39], we consider the same number of requests $N_{o}$ in the
observational window for all the contents, where zero padding is used for
contents with less number of requests. Given the request patterns of contents
over the observational window at time $\tau$, the popularity of contents is
predicted within the time window $[\tau,\tau+T_{s}]$. Without loss of
generality, since the storage capacity of edge devices should be updated each
day, we assume that the length of the study window is $T_{s}=1$. Given $M$
number of multimedia contents and $U$ number of users, the dataset
$\mathcal{D}=\\{d_{\tau}^{m}\\}^{M}_{m=1}$ is created, where
$d_{\tau}^{m}=\\{(x_{k}^{m},t_{k}^{m},y_{k}^{m})\\}^{N_{o}}_{k=1}$ is the
time-series observational data for content $c_{m}$. Term
$\\{x_{k}^{m}\\}^{N_{o}}_{k=1}$ includes the contextual information of users
requesting content $c_{m}$ and the rating that these users gave to content
$c_{m}$ (see Fig. 2(c)). Moreover, rating time is represented by
$\\{t_{k}^{m}\\}^{N_{o}}_{k=1}$. Finally, term $y_{k}^{m}\in\\{0,1\\}$ is
associated with the content popularity, which is $1$ if content $c_{m}$ gets
popular during the study window $T_{s}$, otherwise, $y_{k}^{m}=0$.
Consequently, given $d_{\tau}^{m}$, the output of the proposed CLSA model
predicts whether or not content $c_{m}$ will be popular during the study
window $T_{s}$.

Figure 2: (a) A typical sample of the Movielens dataset, (b) the concatenated
dataset, and (c) the adopted version of the Movielens dataset used for the
CLSA architecture. Figure 3: (a) An illustrative example of request patterns
of contents to construct input samples and their labels. (b) Zero padding
technique is used for contents with the number of requests less than $N_{o}$.

To clarify the benefits of the proposed strategy, an illustrative example of
constructing input samples and their labels is depicted in Fig. 3. In this
example, it is assumed that there are $5$ contents $\\{c_{m}\\}^{5}_{m=1}$
through the network, and the goal is to update the storage of edge devices at
time $\tau+T_{s}=9$ using the historical requests of contents up to time
$\tau=8$. Moreover, $N_{o}=4$ number of requests for each content are
investigated during the observational window. As shown in Fig. 3(a), the
following cases can occur:

  * 1.

In the first case, content $c_{1}$ is requested at time $\tau=8$. Ending at
time $\tau=8$, the observational window $T_{8}^{1}$ is determined in such a
way that the number of requests during this observational window is $N_{o}=4$,
where $T_{8}^{1}=3$. Therefore, the requests within $T_{8}^{1}=3$ ending at
$\tau=8$ are considered to evaluate the popularity of content $c_{1}$ at time
$\tau+T_{s}=9$. Consequently, $\\{x_{k}^{1}\\}^{4}_{k=1}$ includes the
contextual information of users requesting content $c_{1}$ and their ratings.
As shown in Fig. 3(a), since content $c_{1}$ is requested at time $9$, its
label is $\\{y_{k}^{1}\\}^{4}_{k=1}=1$.

  * 2.

In the second case (i.e., content $c_{2}$), there is a request at time
$\tau=8$ and the observational window is $T_{8}^{2}=6$ for having $N_{o}=4$
number of requests within this range. However, since this content is not
requested at time $9$, its label is $\\{y_{k}^{2}\\}^{4}_{k=1}=0$.

  * 3.

In the third case (i.e., content $c_{3}$), the content is requested at
$\tau=8$, but the number of existing requests from the beginning to $\tau=8$
is less than $N_{o}=4$. Therefore, we use zero padding to create the same
length of input samples, shown in Fig 3(b), where zero pad will be added at
the beginning of the observational window.

  * 4.

In the fourth case, content $c_{4}$ is not requested at time $\tau=8$, but
there is at least $N_{o}=4$ number of requests before $\tau=8$. Therefore, the
last $N_{o}=4$ of requests within the observational window $T_{o}^{4}$ are
used as the input sample.

  * 5.

In the fifth case, content $c_{5}$ is not requested at time $\tau=8$, and the
number of requests from the beginning to time $\tau=8$ is less than $N_{o}=4$.
Similar to the case $3$, the zero padding technique is used to create the
input sample.

  * 6.

Finally, if there is a content with zero requests up to time $\tau=8$, this
content will be removed from the study.

This completes the problem description and system modeling, next, we present
the proposed CLSA framework.

##  III Proposed CLSA Framework

Figure 4: The overall perspective of the CLSA architecture.

In this section, we present the constituent components of the proposed CLSA
framework. As can be seen from Fig. 4, the CLSA architecture consists of three
modules, i.e., CL, RN, and SA models. The CA network is utilized to capture
the spatial and temporal correlations of input samples by converting the
longitudinal input data $x_{k}^{m}\in\\{d^{m}_{\tau}\\}^{M}_{m=1}$ to the
latent representations. The RN block is then used to decode the latent
representation. Given a meaningful latent representation, the SA model is used
to predict the probability that content $c_{m}$ will become popular. In the
following, first, we briefly introduce the CL block, then, we explain RN and
the SA models.

###  III-A Contrastive Learning (CL) Block

CL model, as one of the widely used self-supervised learning paradigms in
computer vision, has recently been applied to tabular and longitudinal
datasets. The fundamental goal of the CL network is to learn a latent
representation in which similar samples stay nearby and pairs that are
dissimilar to one another move farther apart. As shown in Fig. 4, among all
samples in each batch, one input sample is considered the anchor sample.
Relying on the data augmentation, a positive pair is created using the anchor
sample, and other samples in the batch are considered negative samples. In
this work, we use the following two data augmentation techniques to generate
positive samples (in what follows, for simplicity, we drop subscript $\tau$
wherever there is no ambiguity):

  * •

Masking: Given the longitudinal data $d_{m}$ of content $c_{m}$, several
users’ information $\\{x_{k}^{m}\\}^{N_{o}}_{k=1}$ is randomly masked to
generate a positive sample for content $c_{m}$, denoted by $d_{m}^{(MA)}$.
Other historical request patterns of contents in the batch are represented by
negative samples.

  * •

Shuffling: In this type of data augmentation, the positive sample is generated
by randomly shuffling the time order of users’ information requesting content
$c_{m}$. The shuffled sample is denoted by $d_{m}^{(SH)}$.

Given the augmented and negative samples, a shared encoder is utilized to
convert the longitudinal input data to a meaningful latent representation.
This representation places contents with similar request patterns closer
together and dissimilar contents farther apart. To preserve the temporal
correlation of request patterns, an LSTM architecture is used as the encoder,
where the latent representation of $x_{k}^{m}$ at time $t_{k}^{m}$, denoted by
$h_{k}^{m}$, is given by

$h_{k}^{m}=LSTM(x_{k}^{m},h_{k-1}^{m}),~{}~{}~{}~{}~{}~{}k=1,\ldots,N_{o}.$
(1)

In addition to the anchor sample $d_{m}$, the shuffled and masked version of
$d_{m}$ are given as the input of the encoder to efficiently learn the latent
representation of $d_{m}$. The shuffled and masked encoded samples are
represented by $h_{k}^{m,(SH)}$, and $h_{k}^{m,(MA)}$, respectively.

Finally, we utilize the masked and shuffled CL loss functions, denoted by
$L_{cl}^{(MA)}$, and $L_{cl}^{(SH)}$, respectively, where the shuffled/masked
learned representation is considered the positive sample, and other
representations in the batch are considered the negative ones, as follows

$L_{cl}^{(MA)}=-\sum\limits_{m=1}^{M}\log\dfrac{\exp\big{(}h_{N_{o}}^{m}(h_{N_{o}}^{m,(MA)})^{T}\big{)}}{\sum\limits_{j=1,j\neq
m}^{M}\exp\big{(}h_{N_{o}}^{m}(h_{N_{o}}^{j,(MA)})^{T}\big{)}},$ (2)
$L_{cl}^{(SH)}=-\sum\limits_{m=1}^{M}\log\dfrac{\exp\big{(}h_{N_{o}}^{m}(h_{N_{o}}^{m,(SH)})^{T}\big{)}}{\sum\limits_{j=1,j\neq
m}^{M}\exp\big{(}h_{N_{o}}^{m}(h_{N_{o}}^{j,(SH)})^{T}\big{)}},$ (3)

where $(.)^{T}$ is the transpose function, and the total CL loss is
$L_{cl}=L_{cl}^{(MA)}+L_{cl}^{(SH)}$. This completes the presentation of the
CL component of the CLSA architecture. Next, we present the reconstruction
network.

###  III-B Reconstruction Network (RN) Block

To learn a better latent representation, the masked version of $x_{k}^{m}$ is
regenerated from the encoded masked input sample. Following Reference [40] and
to preserve the temporal information of $\\{x_{k}^{m}\\}^{N_{o}}_{k=1}$, the
Time-LSTM2 is utilized as the decoder, where the regenerated sample is given
by

$\overline{x}_{k}^{m}=Time-LSTM2(H_{k}^{m}),$ (4)

where $\overline{x}_{k}^{m}$ is the decoded version of $x_{k}^{m}$, and
$H_{k}^{m}=[(h_{1}^{m,(MA)},t_{2}^{m}-t_{1}^{m}),(h_{2}^{m,(MA)},t_{3}^{m}-t_{2}^{m}),\ldots,(h_{k}^{m,(MA)},t_{k+1}^{m}-t_{k}^{m})]$
with $k\in\\{1,\ldots,N_{o}\\}$. Term $h_{k}^{m,(MA)}$ is the masked encoded
sample at time $t_{k}^{m}$, and $t_{k}^{m}-t_{k-1}^{m}$ is the time difference
between two consecutive requests of content $m$. By minimizing the difference
between the original input sample $x_{k}^{m}$ and the decoded one
$\overline{x}_{k}^{m}$, our goal is to provide a better latent representation.
Therefore, the RN loss, denoted by $L_{re}$, is given by

$L_{rn}=\sum_{k=1}^{N_{o}}||\overline{x}_{k}^{m}-x_{k}^{m}||^{2}.$ (5)

###  III-C Survival Analysis (SA) Block

We utilize a SA model to capture a meaningful relationship between the
longitudinal and contextual information of users requesting content $c_{m}$
and the probability of getting popular in the future. To this end, we use an
MLP network, where the input of this model is the latent representation
$\\{h_{N_{o}}^{m}\\}^{M}_{m=1}$, and the output is the probability that
content $c_{m}$ gets popular at time $t$, denoted by
$p_{t}^{m}=P(t|h_{N_{o}}^{m})$, during the total time window, denoted by
$T_{total}$, where $T_{total}=T_{s}+\max\limits_{1\leq m\leq M}T_{o}^{m}$. We
use the softmax function as the activation function of the output layer of the
MLP network to compute $p_{t}^{m}$. Consequently, the estimated Cumulative
Incidence Function (CIF) for content $c_{m}$ at time $\tau+T_{s}$ is
calculated as follows

$\overline{F}^{m}(\tau+T_{s}|h_{N_{o}}^{m})=\dfrac{\sum\limits_{\tau\leq
t\leq\tau+T_{s}}p_{t}^{m}}{1-\sum\limits_{t\leq t_{N_{o}}^{m}}p_{t}^{m}},$ (6)

where $\overline{F}^{m}(\tau+T_{s}|h_{N_{o}}^{m})$ indicates the probability
that content $m$ has gained popularity up to time $\tau+T_{s}$. Finally, the
negative log-likelihood is used as the loss function for the SA network, given
by

$\displaystyle
L_{sa}=-\sum\limits_{m=1}^{M}\log\big{(}1-\overline{F}^{m}(\tau^{m}|h_{T_{m}}^{m})\big{)}.$
(7)

Finally, $\\{p_{t}^{m}\\}^{M}_{m=1}$ are sorted in the descending order to
identify the Top-$K$ popular contents.

##  IV Simulation Results

###  IV-A Experimental Configurations

To investigate the effectiveness of the proposed CLSA architecture, a UAV-
aided MEC network, consisting of $4$ terrestrial and $2$ aerial caching nodes
with $943$ users and $1682$ multimedia contents. Following the common
assumption [2], the size of the storage capacity of caching nodes is $10\%$ of
the total contents, where the size of all multimedia contents is the same.
Using the five-fold cross-validation strategy, $80\%$ and $20\%$ samples are
used as the training dataset and the test dataset, respectively. Adam
optimizer was employed to train the model, where betas are ($0.9,0.999$) and
the weight decay is set to $1e-7$. Moreover, the $l_{2}$ regularization was
set to $1e-3$ to avoid over-fitting. We trained the proposed CLSA architecture
by minimizing the total loss function, denoted by $L_{total}$, obtained as
follows

$L_{total}=\omega_{cl}L_{cl}+\omega_{rn}L_{rn}+\omega_{sa}L_{sa},$ (8)

where $\omega_{cl}$, $\omega_{rn}$, and $\omega_{sa}$ represent the weight of
CL, RN, and SA blocks, respectively, where the summation of them is one. In
the following, the details of each learning block are described:

  * •

Encoder: An LSTM network is used as the encoder, with Rectified Linear Unit
(ReLU) activation and sigmoid as the recurrent activation function, where the
output size of this block is denoted by $D_{\text{E}}$. Then, the output is
given to an MLP with three layers, with the size of $\alpha_{l}D_{\text{FI}}$,
where $l\in\\{1,2,3\\}$. Term $D_{\text{FI}}$ represents the feature dimension
of the input sample, which is $28$ in this work, and the hyperparameter
$\alpha_{l}$, with $l\in\\{1,2,3\\}$ is set to $\alpha_{1}=1,\alpha_{2}=3$,
and $\alpha_{3}=5$.

  * •

Decoder: There is a decoder in the RN block, performing based on the Time-
LSTM2, where the recurrent activation function is sigmoid and the general
activation is Tanh. The output size of the Time-LSTM2 network is denoted by
$D_{\text{D}}$.

  * •

MLP Network: Using batch-normalization technique, this network consists of
three dense layers with the same size of $D_{\text{M}}$ and the Exponential
Linear Unit (ELU) activation function. There are more $5$ dense layers after
that, where the size of each layer is $(5,3,2,1,1)\times T_{total}$. The
activation function of the last layer is softmax, while the rest is ReLU.

###  IV-B Effectiveness of the CLSA Architecture

To evaluate the performance of the proposed CLSA architecture, we first
consider different variants of the CLSA model by changing different
hyperparameters, such as the batch size, $D_{\text{E}}$, $D_{\text{D}}$,
$D_{\text{M}}$, and the learning rate, denoted by $lr$. According to the
information provided in Table I, five different models are defined. Moreover,
we investigate the effect of the number of requests studied in the
observational window, denoted by $N_{o}$, on the classification accuracy. As
it can be seen from Table II, decreasing the batch size from $512$ to $256$
(Model $1$ and $2$) improves the classification accuracy over all folds with
different $N_{o}$. Comparing Models $2$ and $3$, it can be seen that reducing
the encoder and decoder dimensions from $512$ to $256$ decreases the
classification accuracy. Similarly, reducing $D_{\text{M}}$ from $128$ to $32$
results in degrading the accuracy (Model $2$ and Model $4$). Finally, by
comparing Model $1$ and Model $5$, it is evident that $lr=1e-3$ outperforms
$lr=1e-4$.

Moreover, we investigate the effect of $N_{o}$ on the classification accuracy
over $5$ models. According to the results provided in Table II, increasing the
number of requests studied over an observational window provides more
information about the behavior of users’ interests in the past, improving the
classification accuracy. Accordingly, from the aspect of classification
accuracy, it can be seen that Model $2$ with $N_{o}=20$ outperforms other
variants. For this reason, we have selected this model to conduct further
research. Table III demonstrates the precision, recall, and F1-score for Model
$2$ with $N_{o}=20$ over $5$ folds and their average values. According to the
information provided in Table III, the high value of the aforementioned
parameters illustrates the effectiveness of the CLSA architecture.

Fig. 5 represents the confusion matrix of the Model $2$ with $N_{o}=20$. It
should be noted that a challenging issue encountered in the Movielens dataset
pertains to imbalanced data, which arises due to a vast proportion of
multimedia contents being unpopular. To address this problem, we implemented
the random oversampling technique to increase the number of popular samples.
As depicted in Fig. 5, there is a misclassification rate of $2.84\%$ for
unpopular contents (labeled as 0) being incorrectly classified as popular
contents, which leads to the wastage of storage capacity on edge devices.
Similarly, a misclassification rate of $2.06\%$ for popular contents being
classified as unpopular contents can result in failure to place highly
requested contents on edge devices.

Moreover, we utilize the T-distributed Stochastic Neighbor Embedding (TSNE)
method [41] to evaluate the efficiency of the CL block in generating latent
representations for discriminating between popular and unpopular contents. For
instance, in Fig. 6, the latent representation of a test set from one of the
$5$ fold cross-validation experiments is employed to illustrate the embedded
space of popular and unpopular contents. As shown in Fig.6, the embedded space
of popular and unpopular contents are distinguishable, indicating that the CL
network has been effectively trained.

TABLE I: Variants of the CLSA architecture. Model ID | Batch Size | $D_{\text{E}}$ | $D_{\text{D}}$ | $D_{\text{M}}$ | $lr$
---|---|---|---|---|---
1 | 512 | 512 | 512 | 128 | $1e-3$
2 | 256 | 512 | 512 | 128 | $1e-3$
3 | 256 | 256 | 256 | 128 | $1e-3$
4 | 256 | 512 | 512 | 32 | $1e-3$
5 | 256 | 512 | 512 | 128 | $1e-4$
TABLE II: 5 fold cross-validation accuracy $\pm$ standard deviation for different variants of the proposed CLSA architecture using different window sizes (10, 15, and 20 days). $N_{o}$ | Model ID | Fold 1 | Fold 2 | Fold 3 | Fold 4 | Fold 5 | Average
---|---|---|---|---|---|---|---
10 | 1 | $0.850\pm 0.026$ | $0.836\pm 0.037$ | $0.837\pm 0.024$ | $0.875\pm 0.022$ | $0.831\pm 0.050$ | $0.846\pm 0.035$
2 | $0.892\pm 0.031$ | $0.909\pm 0.011$ | $0.885\pm 0.031$ | $0.905\pm 0.016$ | $0.892\pm 0.021$ | $0.896\pm 0.024$
3 | $0.853\pm 0.023$ | $0.829\pm 0.030$ | $0.826\pm 0.025$ | $0.848\pm 0.036$ | $0.815\pm 0.029$ | $0.834\pm 0.031$
4 | $0.864\pm 0.018$ | $0.826\pm 0.012$ | $0.835\pm 0.015$ | $0.849\pm 0.018$ | $0.844\pm 0.031$ | $0.843\pm 0.023$
5 | $0.866\pm 0.012$ | $0.831\pm 0.015$ | $0.834\pm 0.010$ | $0.857\pm 0.012$ | $0.859\pm 0.030$ | $0.849\pm 0.022$
15 | 1 | $0.890\pm 0.031$ | $0.868\pm 0.032$ | $0.843\pm 0.052$ | $0.861\pm 0.027$ | $0.861\pm 0.050$ | $0.865\pm 0.041$
2 | $0.949\pm 0.019$ | $0.940\pm 0.007$ | $0.934\pm 0.010$ | $0.937\pm 0.005$ | $0.932\pm 0.012$ | $0.938\pm 0.012$
3 | $0.861\pm 0.022$ | $0.836\pm 0.024$ | $0.842\pm 0.031$ | $0.861\pm 0.018$ | $0.843\pm 0.033$ | $0.849\pm 0.027$
4 | $0.867\pm 0.026$ | $0.832\pm 0.025$ | $0.835\pm 0.024$ | $0.868\pm 0.023$ | $0.850\pm 0.029$ | $0.850\pm 0.029$
5 | $0.864\pm 0.009$ | $0.826\pm 0.010$ | $0.823\pm 0.012$ | $0.850\pm 0.021$ | $0.870\pm 0.019$ | $0.847\pm 0.024$
20 | 1 | $0.878\pm 0.037$ | $0.881\pm 0.041$ | $0.859\pm 0.050$ | $0.868\pm 0.058$ | $0.855\pm 0.023$ | $0.868\pm 0.043$
2 | $0.959\pm 0.003$ | $0.942\pm 0.006$ | $0.941\pm 0.006$ | $0.946\pm 0.002$ | $0.965\pm 0.003$ | $\textbf{0.951}\pm\textbf{0.010}$
3 | $0.866\pm 0.033$ | $0.875\pm 0.032$ | $0.876\pm 0.019$ | $0.854\pm 0.024$ | $0.845\pm 0.032$ | $0.863\pm 0.029$
4 | $0.858\pm 0.016$ | $0.845\pm 0.039$ | $0.834\pm 0.013$ | $0.869\pm 0.034$ | $0.829\pm 0.065$ | $0.847\pm 0.039$
5 | $0.871\pm 0.009$ | $0.828\pm 0.013$ | $0.830\pm 0.022$ | $0.859\pm 0.009$ | $0.878\pm 0.006$ | $0.853\pm 0.024$
TABLE III: Precision, recall, and F1-score for two classes (i.e., popular (class 1) and unpopular (class 0)) for Model 2 with $N_{o}=20$ using $5$ fold cross-validation. | Class | Fold 1 | Fold 2 | Fold 3 | Fold 4 | Fold 5 | Average
---|---|---|---|---|---|---|---
Precision | 0 | $0.973\pm 0.012$ | $0.934\pm 0.019$ | $0.911\pm 0.005$ | $0.928\pm 0.005$ | $0.956\pm 0.010$ | $0.941\pm 0.024$
1 | $0.946\pm 0.009$ | $0.938\pm 0.019$ | $0.976\pm 0.011$ | $0.951\pm 0.006$ | $0.973\pm 0.005$ | $0.957\pm 0.018$
Recall | 0 | $0.945\pm 0.010$ | $0.937\pm 0.022$ | $0.978\pm 0.010$ | $0.952\pm 0.006$ | $0.974\pm 0.006$ | $0.957\pm 0.019$
1 | $0.974\pm 0.012$ | $0.934\pm 0.021$ | $0.904\pm 0.006$ | $0.926\pm 0.006$ | $0.956\pm 0.010$ | $0.939\pm 0.027$
F1-score | 0 | $0.959\pm 0.004$ | $0.935\pm 0.004$ | $0.943\pm 0.006$ | $0.940\pm 0.005$ | $0.965\pm 0.002$ | $0.948\pm 0.012$
1 | $0.960\pm 0.003$ | $0.935\pm 0.003$ | $0.939\pm 0.003$ | $0.938\pm 0.005$ | $0.964\pm 0.003$ | $0.947\pm 0.013$
Figure 5: Confusion matrix of the proposed CLSA architecture (Model 2,
$N_{o}=20$).  Figure 6: The embedded space of the latent representation of
popular/unpopular contents using the TSNE technique.

###  IV-C Ablation Study

To further investigate the effectiveness of the proposed CLSA architecture, we
conduct an ablation study, where different variants of the CLSA model are
introduced in Table IV. We evaluate the significance of each block in the CLSA
architecture on classification accuracy by analyzing Models $L1$ to $L7$. For
instance, the results presented in Table IV demonstrate that the CL block
alone (Model $L1$) achieves an accuracy of approximately $85\%$, underscoring
the importance of the other two blocks in attaining higher performance (Model
$L7$ with about $95\%$ accuracy). Consequently, it is evident that the
presence of all three blocks in the CLSA architecture is of paramount
importance, where the optimal loss weights of each block are represented in
Model $L7$.

###  IV-D Effectiveness of the CLSA Architecture

In this subsection, we compare the performance of the proposed CLSA
architecture with the following state-of-the-art:

TABLE IV: The accuracy of the proposed CLSA architecture using different loss weights. Model ID | $\omega_{cl}$ | $\omega_{rn}$ | $\omega_{sa}$ | Average Accuracy $\pm$ STD
---|---|---|---|---
$L1$ | 1.0 | 0.0 | 0.0 | $0.851\pm 0.026$
$L2$ | 0.0 | 1.0 | 0.0 | $0.850\pm 0.022$
$L3$ | 0.0 | 0.0 | 1.0 | $0.903\pm 0.029$
$L4$ | 0.5 | 0.5 | 0.0 | $0.859\pm 0.034$
$L5$ | 0.0 | 0.5 | 0.5 | $0.888\pm 0.025$
$L6$ | 0.5 | 0.0 | 0.5 | $0.936\pm 0.013$
$L7$ | 0.3 | 0.2 | 0.5 | $\textbf{0.951}\pm\textbf{0.010}$

  * •

Transformer-based Edge Caching (TEDGE) Scheme [12], which is based on a simple
ViT architecture acting as a multi-label classification model with the aim of
predicting the Top-K popular contents in the upcoming time. To capture the
spatial correlation of contents, $2$D images of historical requests pattern of
contents were created, where the number of columns and rows of this image were
corresponding to the number of contents and the number of historical requests
for each content, respectively. The size of the input sample, therefore,
significantly increases to capture as much spatial correlation as possible.

  * •

Multiple-model Transformer-based Edge Caching (MTEC) [13] consists of two
parallel multi-channel Transformer networks with a dense layer as the fusion
layer. Similarly, the output of the model is Top-K popular contents, while it
first predicted the request patterns of contents in the future. The input
sample is $1$D historical requests patterns of contents. To capture the
spatial correlation of contents, the multi-channel Transformer networks were
employed, where the sequential request pattern of each content is given to a
channel of the Transformer model.

  * •

Vision Transformers with Cross Attention (ViT-CAT) [14], consisting of two
parallel ViT networks with different patching techniques, with a cross
attention mechanism as the fusion layer. The input and output of the network
are similar to the TEDGE caching scheme.

  * •

Self-Supervised Contrastive Learning using Random Feature Corruption (SCARF)
[42], which was applied on tabular datasets. We modified it for the popularity
prediction task, where the input sample is the $1$D historical requests of
content, without any need for labeling. SCARF attempts to learn the latent
representation of contents to classify them as popular and unpopular contents
using random feature corruption.

  * •

Deep Learning-based Content Caching (DLCC) [43] used CNN to predict the
popularity of contents and used the RL model for the content placement phase.

We compare the proposed CLSA framework with the aforementioned baselines from
the classification accuracy perspective in Table V. As shown in Table V, the
proposed CLSA architecture outperforms other baselines, while there is no need
to create large input samples to capture the spatial correlation of contents.

TABLE V: Comparison with state-of-the-art based on the classification accuracy. Model | Accuracy
---|---
TEDGE (ViT) [12] | $93.72\%$
MTEC [13] | $94.13\%$
ViT-CAT [14] | $94.84\%$
SCARF [42] | $87.17\%$
DLCC [43] | $92.81\%$
Proposed CLSA | 95.10 %
Figure 7: Comparison with state-of-the-art based on the cache-hit ratio.

Finally, we compare the performance of the proposed CLSA architecture with
other baselines in terms of the cache-hit ratio. Note that the cache-hit ratio
is commonly used in MEC networks to evaluate the effectiveness of the
popularity prediction framework. This metric shows the number of requests that
are being handled by caching nodes versus the overall number of requests made
throughout the network. As shown in Fig. 7, there are other baselines in
addition to the aforementioned schemes, including Least Recently Used (LRU)
[3], Least Frequently Used (LFU) [3], PopCaching [44], and LSTM-C [6]. As
depicted in Fig. 7, the optimal approach [6] is a caching scheme where caching
nodes handle all requests throughout the network, which is not feasible in
real-world scenarios. Based on the results presented in Fig. 7, the proposed
CLSA architecture achieves the highest cache-hit ratio when compared to other
baselines.

##  V Conclusion

In this paper, we developed the Contrastive Learning-based Survival Analysis
(CLSA) popularity prediction framework with the application to the Mobile Edge
Caching (MEC) networks. To learn the temporal information of sequential
requests, the proposed architecture utilized a self-supervised Contrastive
Learning (CL) model that employed a Long Short-Term Memory (LSTM) network as
the encoder. Unlike existing research works that used multiple contents’
historical request patterns simultaneously to capture spatial dependency, the
proposed CLSA architecture used the CL network, eliminating the need for large
input samples. The output of the architecture was the probabilities of each
content’s future popularity, which were sorted in descending order to
determine the Top-K popular contents. Employing the CL network not only
reduced the input size significantly but also enhanced the scalability of the
learning model. Adding new content to the network does not require recreating
input samples, making it easier to scale up the popularity prediction
framework. Simulation results illustrated that the proposed CLSA architecture
improved the cache-hit ratio and classification accuracy when compared to its
state-of-the-art. Going forward, several directions deserve further
investigation. One area of interest is the implementation of a Contrastive
Learning-based Graph Neural Network (GNN) to establish a meaningful connection
between different users and their interests. Another important aspect is the
role of the augmentation scheme in CL performance. Therefore, developing an
effective augmentation method to create more relevant positive samples is
essential for improving the effectiveness of the popularity prediction
framework.

## References

  * [1] M. Sheraz, M. Ahmed, X. Hou, Y. Li, D. Jin, Z. Han, T. Jiang,  “Artificial Intelligence for Wireless Caching: Schemes, Performance, and Challenges,”  IEEE Communications Surveys & Tutorials, vol. 23, no. 1, pp. 631-661, First quarter 2021.
  * [2] Z. HajiAkhondi-Meybodi, J. Abouei, M. Jaseemuddin and A. Mohammadi,  “Mobility-Aware Femtocaching Algorithm in D2D Networks Based on Handover,”  IEEE Transactions on Vehicular Technology, vol. 69, no. 9, pp. 10188-10201, June 2020.
  * [3] A. Giovanidis, and A. Avranas,  “Spatial multi-LRU caching for wireless networks with coverage overlaps,”  ACM SIGMETRICS Performance Evaluation Review, vol. 44, no. 1, pp. 403-405, June 2016.
  * [4] K. N. Doan, T. Van Nguyen, T. Q. S. Quek, and H. Shin,  “Content-Aware Proactive Caching for Backhaul Offloading in Cellular Network,”  IEEE Transactions on Wireless Communications, vol. 17, no. 5, pp. 3128-3140, May 2018.
  * [5] L. Ale, N. Zhang, H. Wu, D. Chen and T. Han,  “Online Proactive Caching in Mobile Edge Computing Using Bidirectional Deep Recurrent Neural Network,”  IEEE Internet of Things Journal, vol. 6, no. 3, pp. 5520-5530, June 2019.
  * [6] C. Zhang et al.,  “Toward Edge-Assisted Video Content Intelligent Caching With Long Short-Term Memory Learning,”  IEEE Access, vol. 7, pp. 152832-152846, 2019.
  * [7] Y. Lin, C. Yen and J. Wang,  “Video Popularity Prediction: An Autoencoder Approach With Clustering,”  IEEE Access, vol. 8, pp. 129285-129299, 2020.
  * [8] C. Zhong, M. C. Gursoy and S. Velipasalar,  “Deep Reinforcement Learning-Based Edge Caching in Wireless Networks,”  IEEE Transactions on Cognitive Communications and Networking, vol. 6, no. 1, pp. 48-61, March 2020.
  * [9] P. Wu, J. Li, L. Shi, M. Ding, K. Cai and F. Yang,  “Dynamic Content Update for Wireless Edge Caching via Deep Reinforcement Learning,”  IEEE Communications Letters, vol. 23, no. 10, pp. 1773-1777, Oct. 2019.
  * [10] Y. Wang, Y. Li, T. Lan and V. Aggarwal,  “DeepChunk: Deep Q-Learning for Chunk-Based Caching in Wireless Data Processing Networks,”  IEEE Transactions on Cognitive Communications and Networking, vol. 5, no. 4, pp. 1034-1045, Dec. 2019.
  * [11] H. Mou, Y. Liu and L. Wang,  “LSTM for Mobility Based Content Popularity Prediction in Wireless Caching Networks,”  IEEE Globecom Workshops, 2019, pp. 1-6.
  * [12] Z. Hajiakhondi-Meybodi, A. Mohammadi, E. Rahimian, S. Heidarian, J. Abouei, and K. N. Plataniotis,  “TEDGE-Caching: Transformer-based Edge Caching Towards 6G Networks,”  IEEE International Conference on Communications (ICC), Apr. 2022.
  * [13] Z. HajiAkhondi-Meybodi, et al.,  “Multi-Content Time-Series Popularity Prediction with Multiple-Model Transformers in MEC Networks,”  arXiv preprint arXiv:2210.05874, Oct. 2022.
  * [14] Z. HajiAkhondi-Meybodi, A. Mohammadi, M. Hou, J. Abouei, and K. N. Plataniotis.  “ViT-CAT: Parallel Vision Transformers with Cross Attention Fusion for Popularity Prediction in MEC Networks,”  arXiv preprint arXiv:2210.15125, 2022.
  * [15] B. M. Marlin, R. S. Zemel, S. T. Roweis, and M. Slaney,  “Recommender systems: missing data and statistical model estimation,”  Twenty-Second International Joint Conference on Artificial Intelligence, June 2011.
  * [16] A. Odic, M. Tkalcic, J. F. Tasic, and A. Kosir,  “Predicting and detecting the relevant contextual information in a movie-recommender system,”  Interacting with Computers vol. 25, no. 1, pp. 74-90, 2013.
  * [17] S. M. R. Abidi, X. Yonglin, N. Jianyue, W. Xiangmeng, and W. Zhang,  “Popularity prediction of movies: from statistical modeling to machine learning techniques,”  Multimedia Tools and Applications, vol. 79, no. 47, pp. 35583-35617, 2020.
  * [18] G. Mendez, T. D. Buskirk, S. Lohr, and S. Haag,  “Factors associated with persistence in science and engineering majors: An exploratory study using classification trees and random forests,”  Journal of Engineering Education, vol. 97, no. 1, pp. 57-70, 2008.
  * [19] V. K. Y. Ng, and R. A. Cribbie,  “The gamma generalized linear model, log transformation, and the robust Yuen-Welch test for analyzing group means with skewed and heteroscedastic data,”  Comm. in Statistics-Simulation and Computation, vol. 48, no. 8, pp. 2269-2286, 2019.
  * [20] R. R. Kabra, and R. S. Bichkar,  “Performance prediction of engineering students using decision trees,”  International Journal of computer applications, vol. 36, no. 11, pp. 8-12, 2011.
  * [21] Q. Fan, X. Li, J. Li, Q. He, K. Wang and J. Wen,  “PA-Cache: Evolving Learning-Based Popularity- Aware Content Caching in Edge Networks,”  IEEE Transactions on Network and Service Management, vol. 18, no. 2, pp. 1746-1757, June 2021.
  * [22] S. Rathore, J. H. Ryu, P. K. Sharma and J. H. Park,  “DeepCachNet: A Proactive Caching Framework Based on Deep Learning in Cellular Networks,”  IEEE Network, vol. 33, no. 3, pp. 130-138, May/June 2019.
  * [23] K. C. Tsai, L. Wang, and Z. Han,  “Mobile Social Media Networks Caching with Convolutional Neural Network,”  IEEE Wireless Communications and Networking Conference Workshops, 2018, pp. 83-88.
  * [24] A. Ndikumana, N. H. Tran, D. H. Kim, K. T. Kim and C. S. Hong,  “Deep Learning Based Caching for Self-Driving Cars in Multi-Access Edge Computing,”  IEEE Trans. Intell. Transp. Syst., vol. 22, no. 5, pp. 2862-2877, May 2021.
  * [25] J. Tang, H. Tang, X. Zhang, K. Cumanan, G. Chen, K.-K. Wong, and J. A. Chambers,  “Energy Minimization in D2D-Assisted Cache-Enabled Internet of Things: A Deep Reinforcement Learning Approach,”  IEEE Transactions on Industrial Informatics, vol. 16, no. 8, pp. 5412-5423, Aug. 2020.
  * [26] A. Sadeghi, G. Wang and G. B. Giannakis,  “Deep Reinforcement Learning for Adaptive Caching in Hierarchical Content Delivery Networks,”  IEEE Transactions on Cognitive Communications and Networking, vol. 5, no. 4, pp. 1024-1033, Dec. 2019.
  * [27] Y. Jiang, H. Feng, F. -C. Zheng, D. Niyato and X. You,  “Deep Learning-Based Edge Caching in Fog Radio Access Networks,”  IEEE Transactions on Wireless Communications, vol. 19, no. 12, pp. 8442-8454, Dec. 2020.
  * [28] Z. Zhang and M. Tao,  “Deep Learning for Wireless Coded Caching With Unknown and Time-Variant Content Popularity,”  IEEE Transactions on Wireless Communications, vol. 20, no. 2, pp. 1152-1163, Feb. 2021.
  * [29] A. Vaswani, N.Shazeer, N. Parmar, J. Uszkoreit, L. Jones, A. N. Gomez, L. Kaiser, and I. Polosukhin,  “Attention is all you need,”  Advances in neural information processing systems, pp. 5998-6008. 2017.
  * [30] P. H. Le-Khac, G. Healy and A. F. Smeaton,  “Contrastive Representation Learning: A Framework and Review,”  IEEE Access, vol. 8, pp. 193907-193934, 2020.
  * [31] J. L. Katzman, U. Shaham, A. Cloninger, J. Bates, T. Jiang, and Y. Kluger, Y,  “DeepSurv: personalized treatment recommender system using a Cox proportional hazards deep neural network,”  BMC medical research methodology, 18(1), 1-12.
  * [32] H. Jing, and A. J. Smola,  “Neural survival recommender,”  Tenth ACM International Conference on Web Search and Data Mining, 2017.
  * [33] G. Desirena, A. Diaz, J. Desirena, I. Moreno, and D. Garcia,  “Maximizing Customer Lifetime Value using Stacked Neural Networks: An Insurance Industry Application,”  IEEE International Conference On Machine Learning And Applications (ICMLA), Boca Raton, FL, USA, 2019, pp. 541-544.
  * [34] Z. M. Fadlullah and N. Kato,  “HCP: Heterogeneous Computing Platform for Federated Learning Based Collaborative Content Caching Towards 6G Networks,”  IEEE Trans. Emerg. Topics Comput., Apr. 2020.
  * [35] Z. Chen, J. Lee, T. Q. S. Quek and M. Kountouris,  “Cooperative Caching and Transmission Design in Cluster-Centric Small Cell Networks,”  IEEE Trans. Wireless Commun., vol. 16, no. 5, pp. 3401-3415, May 2017.
  * [36] Z. Hajiakhondi-Meybodi, A. Mohammadi and J. Abouei,  “Deep Reinforcement Learning for Trustworthy and Time-Varying Connection Scheduling in a Coupled UAV-Based Femtocaching Architecture,”  IEEE Access, vol. 9, pp. 32263-32281, Feb. 2021.
  * [37] F. M. Harper, and J. A. Konstan,  “The Movielens Datasets: History and Context,”  Acm transactions on interactive intelligent systems, vol. 5, no. 4, pp. 1-19, 2015.
  * [38] G. Vallero, M. Deruyck, W. Joseph and M. Meo,  “Caching at the edge in high energy-efficient wireless access networks,”  IEEE International Conference on Communications (ICC), June 2020, pp. 1-7.
  * [39] C. Hong, F. Yi and Z. Huang,  “Deep-CSA: Deep Contrastive Learning for Dynamic Survival Analysis With Competing Risks,”  IEEE Journal of Biomedical and Health Informatics, vol. 26, no. 8, pp. 4248-4257, Aug. 2022.
  * [40] Y. Zhu, H. Li, Y. Liao, B. Wang, Z. Guan, H. Liu, and D. Cai,  “What to Do Next: Modeling User Behaviors by Time-LSTM,”  IJCAI, vol. 17, pp. 3602-3608, 2017.
  * [41] M. C. Cieslak, A. M. Castelfranco, V. Roncalli, P. H. Lenz, and D.K. Hartline,  “T-Distributed Stochastic Neighbor Embedding (T-SNE): A Tool for Eco-Physiological Transcriptomic Analysis,”  Marine genomics, vol. 51, p. 100723, 2020.
  * [42] D. Bahri, H. Jiang, Y. Tay, D. Metzler,  “Scarf: Self-Supervised Contrastive Learning using Random Feature Corruption,”  arXiv preprint arXiv:2106.15147, 2022.
  * [43] S. Bhandari, N. Ranjan, P. Khan, H. Kim, Y. S. Hong,  “Deep Learning-based Content Caching in the Fog Access Points,”  Electronics, vol. 10, no. 4, pp. 512, 2021.
  * [44] S. Li, J. Xu, M. van der Schaar and W. Li,  “Popularity-Driven Content Caching,”  IEEE International Conference on Computer Communications (INFOCOM), 2016, pp. 1-9.


39 | $\begin{pmatrix}1&0\\\ 0&2\end{pmatrix}$ | $\begin{pmatrix}4&0\\\ 0&1\end{pmatrix}$ | $\begin{pmatrix}1&0\\\ 1&4\end{pmatrix}$ | $\begin{pmatrix}1&0\\\ 7&8\end{pmatrix}$ | $\begin{pmatrix}2&0\\\ 1&2\end{pmatrix}$ | $\begin{pmatrix}1&0\\\ 3&8\end{pmatrix}$
40 | $\begin{pmatrix}1&0\\\ 0&3\end{pmatrix}$ | $\begin{pmatrix}3&0\\\ 0&1\end{pmatrix}$ | $\begin{pmatrix}1&0\\\ 1&6\end{pmatrix}$ | $\begin{pmatrix}1&0\\\ 2&3\end{pmatrix}$ | $\begin{pmatrix}1&0\\\ 4&6\end{pmatrix}$ | $\begin{pmatrix}2&0\\\ 2&3\end{pmatrix}$
41 | $\begin{pmatrix}1&0\\\ 0&3\end{pmatrix}$ | $\begin{pmatrix}3&0\\\ 0&1\end{pmatrix}$ | $\begin{pmatrix}1&0\\\ 1&9\end{pmatrix}$ | $\begin{pmatrix}1&0\\\ 2&3\end{pmatrix}$ | $\begin{pmatrix}1&0\\\ 4&9\end{pmatrix}$ | $\begin{pmatrix}1&0\\\ 7&9\end{pmatrix}$
42 | $\begin{pmatrix}1&0\\\ 0&3\end{pmatrix}$ | $\begin{pmatrix}6&0\\\ 0&1\end{pmatrix}$ | $\begin{pmatrix}1&0\\\ 1&3\end{pmatrix}$ | $\begin{pmatrix}1&0\\\ 2&3\end{pmatrix}$ | $\begin{pmatrix}3&0\\\ 1&2\end{pmatrix}$ | $\begin{pmatrix}3&0\\\ 0&2\end{pmatrix}$
43 | $\begin{pmatrix}1&0\\\ 0&3\end{pmatrix}$ | $\begin{pmatrix}9&0\\\ 0&1\end{pmatrix}$ | $\begin{pmatrix}1&0\\\ 1&3\end{pmatrix}$ | $\begin{pmatrix}1&0\\\ 2&3\end{pmatrix}$ | $\begin{pmatrix}3&0\\\ 1&3\end{pmatrix}$ | $\begin{pmatrix}3&0\\\ 2&3\end{pmatrix}$
44 | $\begin{pmatrix}1&0\\\ 0&6\end{pmatrix}$ | $\begin{pmatrix}3&0\\\ 0&1\end{pmatrix}$ | $\begin{pmatrix}1&0\\\ 1&3\end{pmatrix}$ | $\begin{pmatrix}1&0\\\ 2&3\end{pmatrix}$ | $\begin{pmatrix}1&0\\\ 3&6\end{pmatrix}$ | $\begin{pmatrix}2&0\\\ 0&3\end{pmatrix}$
45 | $\begin{pmatrix}1&0\\\ 0&9\end{pmatrix}$ | $\begin{pmatrix}3&0\\\ 0&1\end{pmatrix}$ | $\begin{pmatrix}1&0\\\ 1&3\end{pmatrix}$ | $\begin{pmatrix}1&0\\\ 2&3\end{pmatrix}$ | $\begin{pmatrix}1&0\\\ 3&9\end{pmatrix}$ | $\begin{pmatrix}1&0\\\ 6&9\end{pmatrix}$
46 | $\begin{pmatrix}1&0\\\ 0&4\end{pmatrix}$ | $\begin{pmatrix}2&0\\\ 0&1\end{pmatrix}$ | $\begin{pmatrix}1&0\\\ 1&4\end{pmatrix}$ | $\begin{pmatrix}1&0\\\ 2&4\end{pmatrix}$ | $\begin{pmatrix}1&0\\\ 3&4\end{pmatrix}$ |
47 | $\begin{pmatrix}1&0\\\ 0&4\end{pmatrix}$ | $\begin{pmatrix}2&0\\\ 0&1\end{pmatrix}$ | $\begin{pmatrix}1&0\\\ 1&4\end{pmatrix}$ | $\begin{pmatrix}1&0\\\ 3&4\end{pmatrix}$ | $\begin{pmatrix}1&0\\\ 2&8\end{pmatrix}$ | $\begin{pmatrix}1&0\\\ 6&8\end{pmatrix}$
48 | $\begin{pmatrix}1&0\\\ 0&8\end{pmatrix}$ | $\begin{pmatrix}2&0\\\ 0&1\end{pmatrix}$ | $\begin{pmatrix}1&0\\\ 1&4\end{pmatrix}$ | $\begin{pmatrix}1&0\\\ 3&4\end{pmatrix}$ | $\begin{pmatrix}1&0\\\ 2&4\end{pmatrix}$ | $\begin{pmatrix}1&0\\\ 4&8\end{pmatrix}$
49 | $\begin{pmatrix}1&0\\\ 0&4\end{pmatrix}$ | $\begin{pmatrix}2&0\\\ 0&1\end{pmatrix}$ | $\begin{pmatrix}1&0\\\ 1&8\end{pmatrix}$ | $\begin{pmatrix}1&0\\\ 3&4\end{pmatrix}$ | $\begin{pmatrix}1&0\\\ 2&4\end{pmatrix}$ | $\begin{pmatrix}1&0\\\ 5&8\end{pmatrix}$
50 | $\begin{pmatrix}1&0\\\ 0&4\end{pmatrix}$ | $\begin{pmatrix}2&0\\\ 0&1\end{pmatrix}$ | $\begin{pmatrix}1&0\\\ 1&4\end{pmatrix}$ | $\begin{pmatrix}1&0\\\ 7&8\end{pmatrix}$ | $\begin{pmatrix}1&0\\\ 2&4\end{pmatrix}$ | $\begin{pmatrix}1&0\\\ 3&8\end{pmatrix}$
51 | $\begin{pmatrix}1&0\\\ 0&3\end{pmatrix}$ | $\begin{pmatrix}3&0\\\ 0&1\end{pmatrix}$ | $\begin{pmatrix}1&0\\\ 1&3\end{pmatrix}$ | $\begin{pmatrix}1&0\\\ 5&6\end{pmatrix}$ | $\begin{pmatrix}1&0\\\ 2&6\end{pmatrix}$ | $\begin{pmatrix}2&0\\\ 1&3\end{pmatrix}$
52 | $\begin{pmatrix}1&0\\\ 0&3\end{pmatrix}$ | $\begin{pmatrix}3&0\\\ 0&1\end{pmatrix}$ | $\begin{pmatrix}1&0\\\ 1&3\end{pmatrix}$ | $\begin{pmatrix}1&0\\\ 8&9\end{pmatrix}$ | $\begin{pmatrix}1&0\\\ 2&9\end{pmatrix}$ | $\begin{pmatrix}1&0\\\ 5&9\end{pmatrix}$
53 | $\begin{pmatrix}1&0\\\ 0&4\end{pmatrix}$ | $\begin{pmatrix}4&0\\\ 0&1\end{pmatrix}$ | $\begin{pmatrix}1&0\\\ 1&4\end{pmatrix}$ | $\begin{pmatrix}1&0\\\ 3&4\end{pmatrix}$ | $\begin{pmatrix}1&0\\\ 2&4\end{pmatrix}$ | $\begin{pmatrix}2&0\\\ 1&2\end{pmatrix}$
54 | $\begin{pmatrix}1&0\\\ 0&5\end{pmatrix}$ | $\begin{pmatrix}5&0\\\ 0&1\end{pmatrix}$ | $\begin{pmatrix}1&0\\\ 1&5\end{pmatrix}$ | $\begin{pmatrix}1&0\\\ 4&5\end{pmatrix}$ | $\begin{pmatrix}1&0\\\ 2&5\end{pmatrix}$ | $\begin{pmatrix}1&0\\\ 3&5\end{pmatrix}$

## References

  * [1] É. Fouvry and P. Koymans, _Binary forms with the same value set I._ arXiv preprint:2404.11231.
  * [2] É. Fouvry and P. Koymans, _Binary forms with the same value set II. The case of ${\bf D}_{4}$._ arXiv preprint:2404.13952.
  * [3] É. Fouvry and M. Waldschmidt, _Number of integers represented by families of binary forms (I)._ Acta Arith. 209 (2023), 219–267.
  * [4] C. L. Stewart, _On the number of solutions of polynomial congruences and Thue equations._ J. Am. Math. Soc. 4 (1991), no. 4, 793–835.
  * [5] C. L. Stewart and S. Y. Xiao, _On the representation of integers by binary forms._ Math. Ann. 375 (2019), no. 1–2, 133–163.


# Te Vacancy-Driven Anomalous Transport in ZrTe5 and HfTe5

Elizabeth A. Peterson  Theoretical Division, Los Alamos National Laboratory,
Los Alamos, NM 87545, USA  Christopher Lane  Theoretical Division, Los Alamos
National Laboratory, Los Alamos, NM 87545, USA  Jian-Xin Zhu  Theoretical
Division, Los Alamos National Laboratory, Los Alamos, NM 87545, USA  Center
for Integrated Nanotechnologies, Los Alamos National Laboratory, Los Alamos,
NM 87545, USA

###### Abstract

In the search for experimental signatures of quantum anomalies, the layered
Dirac materials ZrTe5 and HfTe5 have received much attention for potentially
hosting a chiral anomaly. These materials exhibit a negative longitudinal
magnetoresistance (NLMR) that is taken as a signature of broken chiral
symmetry. The anomalous transport properties of ZrTe5 and HfTe5 are known to
strongly correlate with the presence of Te vacancies, prompting questions as
to the microscopic mechanism driving the NLMR. In this work, the effect of Te
vacancies on the electronic structure of ZrTe5 and HfTe5 is investigated via
first-principles calculations to garner insight into how they may modulate the
transport properties of these materials. While Te vacancies act as a source of
effective compressive strain, they also produce local changes to the
electronic structure that cannot be explained simply as volume effects. The
reorganization of the electronic structure near the Fermi energy indicates
that Te vacancies can rationalize both spectroscopic and transport
measurements that have remained elusive in prior first-principles studies.
These results show that Te vacancies contribute, at least in part, to the
anomalous transport properties of ZrTe5 and HfTe5 and offer a path towards
understanding the possibility of a chiral anomaly in these materials.

##  I Introduction

Theoretical predictions of quantum anomalies, symmetry breaking that occurs
when moving from classical field theories to quantum field theories,
originated in the fields of particle and cosmological physics Adler (1969);
Bell and Jackiw (1969); Kharzeev (2006, 2014). As experimental verification of
predicted quantum anomalies at the relevant length and energy scales of
particle and cosmological physics is generally challenging or simply
intractable, the possibility of observing quantum anomalies in experimentally
accessible condensed matter systems has engendered a great deal of excitement.
In the case of three-dimensional (3D) Dirac and Weyl semimetals, a chiral
anomaly is characterized by an imbalance in right- and left-handed chiral
fermions in the presence of an applied magnetic field. The resulting chiral
current may be induced by application of parallel electric and magnetic fields
producing an anomalous contribution to the conductivity Nielsen and Ninomiya
(1983); Fukushima  _et al._ (2008); Son and Spivak (2013); Kharzeev (2014);
Huang  _et al._ (2015); Xiong  _et al._ (2015); Arnold  _et al._ (2016); Zhang
_et al._ (2016); Li  _et al._ (2016a); Chi  _et al._ (2017); Nandy  _et al._
(2017); Nag and Nandy (2020); Nandy  _et al._ (2021); Xie  _et al._ (2022).
This manifests as a negative longitudinal magnetoresistance (NLMR). In recent
years the layered 3D Dirac materials ZrTe5 and HfTe5 have garnered attention
as potential material platforms that host a chiral anomaly Li  _et al._
(2016a); Zheng  _et al._ (2016); Chi  _et al._ (2017). These materials exhibit
a number of anomalous transport properties, including a Lifshitz transition as
a function of temperature and a negative longitudinal magnetoresistance Chen
_et al._ (2015a); Li  _et al._ (2016a); Liu  _et al._ (2016); Zheng  _et al._
(2016); Chi  _et al._ (2017); Galeski  _et al._ (2020); Zhang  _et al._
(2020); Xie  _et al._ (2022); Jo  _et al._ (2023). This NLMR is suggested to
be evidence of the chiral anomaly in these materials Li  _et al._ (2016a);
Zheng  _et al._ (2016); Chi  _et al._ (2017). However, NLMR is not sufficient
to prove the presence of a chiral anomaly, as it may arise from alternative
sources such as current jetting or defects Arnold  _et al._ (2016); dos Reis
_et al._ (2016); Li  _et al._ (2017). Lacking a robust microscopic description
of the source of the anomalous transport properties in ZrTe5 and HfTe5, it is
as yet unclear if these materials exhibit a genuine chiral anomaly. Further,
there is widespread controversy in the characterization of the topological
nature of these materials; different experiments report different topological
phases, including topological insulating and Dirac semimetallic states Weng
_et al._ (2014); Chen  _et al._ (2015b, a); Liu  _et al._ (2016); Li  _et al._
(2016a); Zheng  _et al._ (2016); Li  _et al._ (2016b); Wu  _et al._ (2016);
Moreschini  _et al._ (2016); Manzoni  _et al._ (2016).

The sample dependence of the topological Manzoni  _et al._ (2016); Fan  _et
al._ (2017); Zhang  _et al._ (2017); Monserrat and Narayan (2019); Jo  _et
al._ (2023) and transport Shahi  _et al._ (2018); Jiang  _et al._ (2023); Jo
_et al._ (2023) properties of ZrTe5 and HfTe5 are well documented and are
suspected to depend, at least in part, on sample volume. Topological phase
transitions from a strong topological insulator, to a Dirac semimetal, to a
weak topological insulator are theoretically predicted to occur with
increasing volume Manzoni  _et al._ (2016); Zhang  _et al._ (2017); Fan  _et
al._ (2017); Monserrat and Narayan (2019); Jo  _et al._ (2023). The Lifshitz
transition (experimentally characterized by a change from p-type to n-type
carriers and from metallic to semiconducting to metallic character with
decreasing temperature as measured through Hall and longitudinal conductivity,
respectively) that is widely featured in the literature on ZrTe5 and HfTe5
does not occur in all samples Shahi  _et al._ (2018); Jiang  _et al._ (2023);
Jo  _et al._ (2023). Puzzlingly, samples that are highly stoichiometric with a
1:5 ratio of Zr or Hf to Te generally do not exhibit strong signals of
anomalous transport behavior; in fact, Te deficiency appears to be critical to
the observation of the Lifshitz transition Shahi  _et al._ (2018); Jo  _et
al._ (2023). Recent experimental and theoretical results using angle-resolved
photoemission spectroscopy (ARPES) and first-principles calculations
rationalize the observed sample dependence of topological phases by reporting
that Te vacancies are a source of chemical pressure, or internal strain,
altering the volume of ZrTe5 and HfTe5 Jo  _et al._ (2023).

However, outstanding discrepancies between low temperature spectroscopic and
transport measurements of ZrTe5 and HfTe5 suggest that the role of Te
vacancies is more complex. Below the Lifshitz transition temperature, Hall
conductivity measurements indicate that the primary carriers in these
materials are electrons, suggesting the conduction band must be populated Liu
_et al._ (2016); Chi  _et al._ (2017); Shahi  _et al._ (2018). While most
reported ARPES measurements show that as temperature decreases the Fermi
energy shifts downward into the valence band at the $\Gamma$ point, suggesting
the primary carriers are holes Li  _et al._ (2016a); Wu  _et al._ (2016);
Moreschini  _et al._ (2016); one ARPES study observed that as temperature
decreases the Fermi energy shifts upward into the conduction band in the
neighborhood of the $\Gamma$ point suggesting the primary carriers are
electrons, consistent with Hall conductivity measurements Zhang  _et al._
(2017). This study went on to measure the full Brillouin zone (BZ),
identifying electron pockets far from the BZ center. Electronic structure
calculations on pristine ZrTe5 and HfTe5 do not predict any portion of the
conduction band dipping below the Fermi energy, an apparent contradiction.
Moreover, ARPES measurements generally suggest the Dirac point is centered at
the $\Gamma$ point while electronic structure calculations of the pristine
materials generally place the Dirac point along the $\Gamma$-Y high-symmetry
line (or $\Gamma$-Z depending on convention) Li  _et al._ (2016a); Wu  _et
al._ (2016); Moreschini  _et al._ (2016); Zhang  _et al._ (2017); Shahi  _et
al._ (2018).

The observation that the Lifshitz transition and NLMR occur only when Te
vacancies are present strongly indicates that Te vacancies serve as more than
just a source of effective strain in altering the electronic structure and
transport properties of ZrTe5 and HfTe5. However, the subtle effect of Te
vacancies on the electronic structure of these materials is poorly understood.
In this work first-principles density functional theory (DFT) calculations of
the electronic structure of ZrTe5 and HfTe5 in the presence of Te vacancies
offer insight into the role that Te vacancies play in modulating the
electronic structure of these materials. These calculations reveal that while
Te vacancies do indeed serve as a source of chemical pressure, or effective
strain, they also significantly modify the electronic structure of ZrTe5 and
HfTe5; understanding these modifications is critical to elucidating whether
the transport properties of these materials are a signature of the presence of
a chiral anomaly.

Figure 1: (a) The conventional crystal structure of HfTe5 characterized by Hf
ion (teal) polyhedron centers each coordinated to 8 Te ions (yellow). Layers
of HfTe8 polyhedra are formed by Hf-Te bonded chains along the $\mathbf{a}$
direction and Te-Te bonded chains along the $\mathbf{c}$ directions. Layers
are stacked along the $\mathbf{b}$ direction held together by vdW dispersion
forces. (b) The isostructural conventional crystal structure of ZrTe5 with Zr
ions in purple and Te ions in yellow. The three symmetrically distinct Te ion
sites are indicated. (c) Orthorhombic Brillouin zone with high-symmetry points
used in the band structure calculations shown in Figure 2-4 marked.

##  II Electronic Structure Calculations

| Source | a (Å) | b (Å) | c (Å)
---|---|---|---|---
ZrTe5 | ICSD | 3.987 | 14.530 | 13.724
DFT | 4.026 | 14.793 | 13.636
HfTe5 | ICSD | 3.968 | 14.455 | 13.691
DFT | 4.000 | 14.711 | 13.596
Table 1: Lattice constants of the experimental crystal structures of ZrTe5 and HfTe5 from the ICSD Fjellvåg and Kjekshus (1986); Furuseth  _et al._ (1973) and the fully geometrically relaxed crystal structures calculated using DFT with Grimme-D3 dispersion corrections. | V${}_{\textrm{Te}}$ Site | Ef (meV/atom) |  $\Delta$V ($\%$) |  $\Delta$a ($\%$) |  $\Delta$b ($\%$) |  $\Delta$c ($\%$)
---|---|---|---|---|---|---
ZrTe5 | 1 | 83.1 | -0.14$\%$ | -0.10$\%$ | 0.00$\%$ | -0.03$\%$
2 | 88.2 | -1.42 | -1.01 | -0.07 | -0.35
3 | 97.6 | -1.02$\%$ | -0.19 | 0.00 | -0.83
HfTe5 | 1 | 77.0 | -0.08 | -0.12 | +0.06 | -0.02
2 | 83.5 | -1.42 | -1.09 | -0.03 | -0.31
3 | 94.8 | -1.10 | -0.27 | +0.11 | -0.93
Table 2: Defect formation energies and changes in the volume and lattice
constants for ZrTe5 and HfTe5 with Te vacancies at each symmetrically distinct
site. Changes in the volume and lattice constants are for the fully relaxed
defect structure and reported relative to the volume and lattice constants of
the fully relaxed pristine crystal structure.

ZrTe5 and HfTe5 are layered materials with one-dimensional chains of HfTe8
polyhedra along the $\mathbf{a}$ axis connected by Te-Te bonds along the
$\mathbf{c}$ axis. These ac-planes are bound by vdW dispersion forces along
the $\mathbf{b}$ axis, as shown in Fig. 1(a). The lattice constants of ZrTe5
and HfTe5 are similar, as shown in Table 1. They have three symmetrically
distinct Te sites, as shown in Fig. 1(b). Prior experiments suggests that Te
vacancies form readily on sites 2 and 3 in ZrTe5 Shahi  _et al._ (2018),
however the calculated formation energies of Te vacancies on each site
indicate that a vacancy on site 1 is the most favorable, in agreement with
recent work Jo  _et al._ (2023). The DFT calculated formation energies for Te
vacancies on each site are all relatively high, ranging from 77-98 meV/atom,
significantly higher than average thermal energy $k_{b}T$ ($\sim$ 26 meV) at
room temperature, as seen in as shown in Table 2. Nonetheless, under high
temperature or non-equilibrium synthesis conditions a finite population of Te
vacancies would still be expected to form. Notably, the relative formation
energy differences between different vacancy sites are on the order of 10-20
meV/atom, meaning that a mixture of vacancy types should be present even
though site 1 is the most stable.

Figure 2: The DFT band structures near the Fermi energy of ZrTe5(a-d) and
HfTe5(e-h) in the fully relaxed geometries for the pristine crystal structure
(a,e) and the cases of a neutral Te vacancy in Te site 1 (b,f), Te site 2
(c,g), and Te site 3 (d,h). The band structures calculated both with (solid
lines) and without (dashed lines) spin-orbit coupling (SOC) are plotted to
illustrate how SOC affects the band structure around the Dirac point near the
Fermi energy.

###  II.1 Band Structures

The band structures of pristine ZrTe5 and HfTe5 are both characterized by a
Dirac point at the Fermi energy that opens into a narrow band gap between the
$\Gamma$ and Y high-symmetry points due to spin-orbit coupling (SOC) as shown
in Fig. 2(a,e). The PBE+SOC calculated indirect band gaps of pristine ZrTe5
and HfTe5 are 27 meV and 24 meV respectively. When Te vacancies are
introduced, the electronic structure calculated with SOC looks very similar
for both ZrTe5 and HfTe5. This is likely a reflection of the very close
similarity in lattice parameters between the two isostructural materials.

The PBE+SOC band structure for Te vacancies at site 1 (Fig. 2(b,f)) is quite
similar to the pristine cases except with (i) the direct band gap at the Dirac
point shifted from being along the $\Gamma$-Y high-symmetry line to being
exactly at the $\Gamma$ point and (ii) a decrease of the conduction band (CB)
energy below the Fermi level along the Z-U and R-T high-symmetry lines. In
ZrTe5, the Fermi energy still lies within the Dirac point gap while in HfTe5
the Fermi energy cuts through the valence band (VB) near the Dirac point at
$\Gamma$. Te vacancies at site 2 (Fig. 2(c,g)) shift the Fermi energy below
the VB top at the Dirac point and also shift the CB below the Fermi energy in
several more regions of the Brillouin zone (BZ). Te vacancies at site 3 (Fig.
2(d,h)) shift the CB and VB in opposite directions both towards the Fermi
energy.

Figure 3: The DFT+SOC band structures and densities of states for ZrTe5 in the
fully relaxed geometries for the pristine crystal structure (a,e) and a
neutral Te vacancy in Te site 1 (b,f), Te site 2 (c,g), and Te site 3 (d,h).
Species projected orbital contributions are plotted for Zr (purple) and Te
(yellow). The Dirac point is circled in each band structure to illustrate
shifts in its location relative to the Fermi energy and changes in gap size
and dispersion. The density of states for Te vacancies on each symmetrically
distinct site are plotted for the Zr polyhedron where the Te vacancy is
located (Near VTe, solid red line) and the Zr polyhedron farthest from the Te
vacancy (Far from VTe, black dashed line). Figure 4: Same as Fig. 3 for HfTe5.
The DFT+SOC band structures and densities of states for HfTe5 in the fully
relaxed geometries for the pristine crystal structure (a,e) and a neutral Te
vacancy in Te site 1 (b,f), Te site 2 (c,g), and Te site 3 (d,h). Species
projected orbital contributions are plotted for Hf (teal) and Te (yellow). The
Dirac point is circled in each band structure to illustrate shifts in its
location relative to the Fermi energy and changes in gap size and dispersion.
The density of states for Te vacancies on each symmetrically distinct site are
plotted for the Hf polyhedron where the Te vacancy is located (Near VTe, solid
red line) and the Hf polyhedron farthest from the Te vacancy (Far From VTe,
black dashed line).

In all cases, Te vacancies make the band structures of ZrTe5 and HfTe5
metallic, although least dramatically for Te vacancies at site 1, the most
favorable Te vacancy. The volume reduction caused by Te vacancies at site 3 is
the largest ($>1\%$), followed by site 2, while the volume reduction caused by
Te vacancies at site 1 is much less dramatic ($<0.1\%$)(see Table 2 for
details). This explains the trend in metallicity across the band structures of
each type of Te vacancy, with more effective compressive strain resulting in a
more metallic band structure.

Interestingly, at $\Gamma$, the direct band gap is larger in the presence of
Te vacancies than it is for the pristine crystal structures. This is
consistent with recent band structure calculations of HfTe5 under compressive
strain Jo  _et al._ (2023).

The band gap at the Dirac point appears both when ZrTe5 and HfTe5 are pristine
and when they host Te vacancies in our PBE+SOC band structures. As noted above
though, while the band gap for the pristine case is along the $\Gamma$-Y high-
symmetry line, the band gap shifts to being exactly at the $\Gamma$ point when
Te vacancies are present. Recent first-principles calculations on HfTe5 under
compressive strain do not report this qualitative change, indicating this is
not merely a volume effect Jo  _et al._ (2023). Experimental ARPES studies
consistently find the Dirac point centered at $\Gamma$ Li  _et al._ (2016a);
Moreschini  _et al._ (2016); Wu  _et al._ (2016); Zhang  _et al._ (2017),
contrary to DFT calculations of the pristine crystal structures. The suspected
omnipresence of Te vacancies of variable concentration may offer an
explanation for the experimental observation that the band gap is located at
the $\Gamma$ point.

###  II.2 Site Projected Densities of States

To better understand the source of the energy shifts and reordering of the
bands near the Fermi energy, the site-projected density of states is
calculated for ions in regions near the Te vacancy site and ions in regions
far from the Te vacancy site. Fig. 3(a-d) shows the band structures of
pristine ZrTe5 and ZrTe5 with Te vacancies plotted with orbital contributions
from each element indicated by purple and yellow dots for Zr and Te
respectively. The band gap at the Dirac point, circled in orange, is
consistently characterized by a mix of contributions from Zr and Te orbitals.
The conduction bands that shift below the Fermi energy are predominantly of Te
character. One notable exception is a flat band of Zr character that lies
right at the Fermi energy along the U-R high-symmetry line for Te vacancies at
site 3, which would be expected to have important implications for the
transport properties of this case. Nearly identical observations can be made
for the band structures of HfTe5 as shown in Fig. 4(a-d).

The local site projected density of states for the Zr and Te ions in the Zr
polyhedron furthest away from the Te vacancy (Far from V${}_{\textrm{Te}}$)
strongly resembles the density of states of pristine ZrTe5 (as shown in Fig.
3(e-h)), with the exception of a notable shift in the position of the band gap
for Te vacancies on site 2. Conversely, the local site projected density of
states for the ions in the Zr polyhedron where the Te vacancy is located (Near
V${}_{\textrm{Te}}$) indicates that orbitals near the Te vacancy are the
primary source of in-gap states. A similar conclusion is drawn from the site
projected densities of states both near and far from the Te vacancies for
HfTe5 (Fig. 4(e-h)).

##  III Discussion

The electronic structure calculations are of critical importance in
rationalizing the experimentally observed low temperature transport and
spectroscopic properties of ZrTe5 and HfTe5, as these DFT calculations
correspond to zero temperature. Near the $\Gamma$ point in the electronic
structure, the Fermi energy either remains in the band gap at the Dirac point
or cuts through the valence band (as in the case of Te vacancies at site 2).
This is consistent with most of the experimental ARPES observations near the
$\Gamma$ point Li  _et al._ (2016a); Wu  _et al._ (2016); Moreschini  _et al._
(2016). However, when considering the entire Brillouin zone, it becomes clear
that there is significant population of the conduction bands throughout
regions of the BZ that are further away from $\Gamma$ when Te vacancies are
present. This is consistent with more comprehensive ARPES measurements that
observed electron pockets along the Y-S high-symmetry line Zhang  _et al._
(2017). Combined, these calculations offer an explanation for the seemingly
conflicting results of ARPES measurements and Hall conductivity measurements.

These results support and help rationalize the observation that Te vacancies
contribute to the anomalous transport properties of ZrTe5 and HfTe5. While Te
vacancies do serve as a source of chemical pressure, or effective strain, by
reducing the volume, this does not account for the notable differences in the
local site projected density of states in proximity to Te vacancy sites. The
local density of states of regions far from the Te vacancies capture the
effects of volume modulation, including shifts of the Fermi energy relative to
the Dirac point and changes in the size of the band gap at the Dirac point.
The local density of states of regions near the Te vacancies capture the
effects of the Te vacancies themselves, most notably the introduction of
additional states at the Fermi level that fundamentally alter the band
structure and drive a transition to a metallic state.

##  IV Conclusions

In this work, first-principles DFT electronic structure calculations of the
layered Dirac materials ZrTe5 and HfTe5 with Te vacancies offer insight into
the sample dependence of their experimentally observed spectroscopic and
anomalous transport properties. Consistent with most low temperature ARPES
measurements, Te vacancies are shown to promote both semiconducting or hole-
doped behavior at the Dirac point, which is shifted to the $\Gamma$ point.
They further promote occupation of the conduction band by shifting segments to
energies below the Fermi energy which rationalizes the n-type conductivity at
low temperatures consistently reported by Hall measurements. Different
combinations of Te vacancies at different concentrations are expected to
produce a spectrum of electronic structure reorganizations accounting for
sample-dependent reports of gapped and metallic behavior. Te vacancies affect
the electronic structure and corresponding transport properties of ZrTe5 and
HfTe5 both as a general source of effective strain and, most importantly, by
introducing additional states near the Fermi energy localized in the vicinity
of the Te vacancy sites. The qualitative character of the Dirac point in the
DFT+SOC calculations is relatively robust to the introduction of Te vacancies,
only changing by shifting towards the $\Gamma$ point and exhibiting a larger
band gap. This suggests that the proposed mechanism of a chiral anomaly
occurring in ZrTe5 and HfTe5 deriving from the topological nature of Dirac
materials may still hold, however not via the simple picture of chiral
fermions originating strictly at the Dirac point of the pristine structures.
As the anomalous transport properties that suggest a chiral anomaly are tied
to the presence of Te vacancies, further analysis of the topological character
of the shifted Dirac point and the newly occupied conduction bands will reveal
their significance in potentially hosting a quantum anomaly in experimentally
accessible condensed matter systems.

##  V Calculation Details

First-principles calculations are performed using density functional theory
(DFT) with a plane-wave basis and projector augmented wave (PAW)
pseudopotentials Kresse and Joubert (1999) as implemented in the Vienna ab
initio simulation package (VASP) Kresse and Furthmüller (1996a, b).
Calculations are performed in the generalized gradient approximation (GGA) as
implemented by Perdew, Burke, and Ernzerhof (PBE) Perdew  _et al._ (1996) with
additional vdW dispersion forces approximately accounted for via the Grimme-D3
method Grimme  _et al._ (2010). The crystal structures of bulk ZrTe5 and HfTe5
are relaxed using a 600 eV energy cutoff and 20x8x8 $\Gamma$-centered k-mesh
until forces are converged to $<$ 1 meV/Å. This method results in good
agreement of the $\mathrm{a}$ and $\mathrm{c}$ lattice constants (Table 1).
The out-of-plane $\mathrm{b}$ lattice constants are only slightly
overestimated by 1.8$\%$ and 1.7$\%$ for ZrTe5 and HfTe5 respectively. The
band structure is calculated both with and without spin-orbit coupling (SOC)
using a 500 eV energy cutoff. The density of states is calculated with SOC.

Te vacancy calculations are performed on 2x1x1 supercells of bulk ZrTe5 and
HfTe5 with a single Te ion removed, a 2.5$\%$ Te vacancy concentration.
Vacancies are introduced at each of the three symmetrically distinct Te sites.
Geometry optimization is performed both for fixed lattice parameters with only
internal coordinates allowed to relax as well as full crystal structure
relaxation. For the supercell calculations, an energy cut-off of 500 eV and a
10x8x8 $\Gamma$-centered k-mesh are used. The band structure is calculated for
the fully relaxed crystal structures, to capture volume effects, both with and
without spin-orbit coupling (SOC) and the density of states is calculated with
SOC.

The defect formation energy is calculated using

$E_{f}=E_{tot}[\textrm{V}_{\textrm{Te}}]-E_{tot}[\textrm{pristine}]-\Sigma_{i}n_{i}\mu_{i}$
(1)

where $E_{tot}[\textrm{V}_{\textrm{Te}}]$ and $E_{tot}[\textrm{pristine}]$ are
the DFT calculated total energies of the crystal structure with a Te vacancy
and the pristine crystal structure. $n_{i}$ is the number of atoms of species
$i$ added or removed (in this case one Te atom) and $\mu_{i}$ is the chemical
potential of that species, here taken from a DFT calculation of bulk Te in the
trigonal P3121 space group as tabulated in the ICSD Adenis  _et al._ (1989).
The standard additional terms to correct for charge effects are neglected
because only neutral Te vacancies are considered Van de Walle and Neugebauer
(2004).

##  VI Acknowledgements

This work was supported by the U.S. DOE NNSA under Contract No.
89233218CNA000001. It was supported by the LANL LDRD Program, and in part by
the Center for Integrated Nanotechnologies, an Office of Science User Facility
operated by the U.S. Department of Energy (DOE) Office of Science, in
partnership with the LANL Institutional Computing Program for computational
resources. Additional computations were performed at the National Energy
Research Scientific Computing Center (NERSC), a U.S. Department of Energy
Office of Science User Facility located at Lawrence Berkeley National
Laboratory, operated under Contract No. DE-AC02-05CH11231 using NERSC award
ERCAP0020494.

## References

  * Adler (1969) S. L. Adler, Axial-Vector Vertex in Spinor Electrodynamics, Phys. Rev. 177, 2426 (1969).
  * Bell and Jackiw (1969) J. S. Bell and R. Jackiw, A PCAC puzzle: $\pi^{0}$ $\rightarrow$ $\gamma\gamma$ in the $\sigma$-model, Il Nuovo Cimento A 60, 47 (1969).
  * Kharzeev (2006) D. Kharzeev, Parity violation in hot QCD: Why it can happen, and how to look for it, Physics Letters B 633, 260 (2006).
  * Kharzeev (2014) E. Kharzeev, Dmitri, The Chiral Magnetic Effect and anomaly-induced transport, Progress in Particle and Nuclear Physics 75, 133 (2014).
  * Nielsen and Ninomiya (1983) H. B. Nielsen and M. Ninomiya, The Adler-Bell-Jackiw anomaly and Weyl fermions in a crystal, Physics Letters B 130, 389 (1983).
  * Fukushima  _et al._ (2008) K. Fukushima, D. E. Kharzeev, and H. J. Warringa, Chiral magnetic effect, Phys. Rev. D 78, 074033 (2008).
  * Son and Spivak (2013) D. T. Son and B. Z. Spivak, Chiral anomaly and classical negative magnetoresistance of Weyl metals, Phys. Rev. B 88, 104412 (2013).
  * Huang  _et al._ (2015) X. Huang, L. Zhao, Y. Long, P. Wang, D. Chen, Z. Yang, H. Liang, M. Xue, H. Weng, Z. Fang, X. Dai, and G. Chen, Observation of the Chiral-Anomaly-Induced Negative Magnetoresistance in 3D Weyl Semimetal TaAs, Phys. Rev. X 5, 031023 (2015).
  * Xiong  _et al._ (2015) J. Xiong, S. K. Kushwaha, T. Liang, J. W. Krizan, M. Hirschberger, W. Wang, R. J. Cava, and N. P. Ong, Evidence for the chiral anomaly in the Dirac semimetal Na3Bi, Science 350, 413 (2015), https://www.science.org/doi/pdf/10.1126/science.aac6089 .
  * Arnold  _et al._ (2016) F. Arnold, C. Shekhar, S.-C. Wu, Y. Sun, R. D. dos Reis, N. Kumar, M. Naumann, M. O. Ajeesh, M. Schmidt, A. G. Grushin, J. H. Bardarson, M. Baenitz, D. Sokolov, H. Borrmann, M. Nicklas, C. Felser, E. Hassinger, and B. Yan, Negative magnetoresistance without well-defined chirality in the Weyl semimetal TaP, Nat. Commun. 7, 10.1038/ncomms11615 (2016).
  * Zhang  _et al._ (2016) C.-L. Zhang, S.-Y. Xu, I. Belopolski, Z. Yuan, Z. Lin, B. Tong, G. Bian, N. Alidoust, C.-C. Lee, S.-M. Huang, T.-R. Chang, G. Chang, C.-H. Hsu, H.-T. Jeng, M. Neupane, D. S. Sanchez, H. Zheng, J. Wang, H. Lin, C. Zhang, H.-Z. Lu, S.-Q. Shen, T. Neupert, M. Z. Hasan, and S. Jia, Signatures of the Adler-Bell-Jackiw chiral anomaly in a Weyl fermion semimetal, Nat. Commun. 7, 10.1038/ncomms10735 (2016).
  * Li  _et al._ (2016a) Q. Li, D. E. Kharzeev, C. Zhang, Y. Huang, I. Pletikosić, A. V. Fedorov, R. D. Zhong, J. A. Schneeloch, G. D. Gu, and T. Valla, Chiral magnetic effect in ZrTe5, Nature Phys. 12, 550 (2016a).
  * Chi  _et al._ (2017) H. Chi, C. Zhang, G. Gu, D. E. Kharzeev, X. Dai, and Q. Li, Lifshitz transition mediated electronic transport anomaly in bulk ZrTe5, New J. Phys. 19, 10.1088/1367-2630/aa55a3 (2017).
  * Nandy  _et al._ (2017) S. Nandy, G. Sharma, A. Taraphder, and S. Tewari, Chiral Anomaly as the Origin of the Planar Hall Effect in Weyl Semimetals, Phys. Rev. Lett. 119, 176804 (2017).
  * Nag and Nandy (2020) T. Nag and S. Nandy, Magneto-transport phenomena of type-I multi-Weyl semimetals in co-planar setups, Journal of Physics: Condensed Matter 33, 075504 (2020).
  * Nandy  _et al._ (2021) S. Nandy, C. Zeng, and S. Tewari, Chiral anomaly induced nonlinear Hall effect in semimetals with multiple Weyl points, Phys. Rev. B 104, 205124 (2021).
  * Xie  _et al._ (2022) Q. Xie, C. Wang, S. Yan, L. Chen, J. Zheng, and W. W., Quantum oscillations and stacked quantum Hall effect in HfTe5, Appl. Phys. Lett. 120, 10.1063/5.0084727 (2022).
  * Zheng  _et al._ (2016) G. Zheng, J. Lu, X. Zhu, W. Ning, Y. Han, H. Zhang, J. Zhang, C. Xi, J. Yang, H. Du, K. Yang, Y. Zhang, and M. Tian, Transport evidence for the three-dimensional Dirac semimetal phase in $\mathrm{ZrT}{\mathrm{e}}_{5}$, Phys. Rev. B 93, 115414 (2016).
  * Chen  _et al._ (2015a) R. Y. Chen, S. J. Zhang, J. A. Schneeloch, C. Zhang, Q. Li, G. D. Gu, and N. L. Wang, Optical spectroscopy study of the three-dimensional Dirac semimetal ${\mathrm{ZrTe}}_{5}$, Phys. Rev. B 92, 075107 (2015a).
  * Liu  _et al._ (2016) Y. Liu, X. Yuan, C. Zhang, Z. Jin, Q. Narayan, C. Luo, Z. Chen, L. Yang, J. Zou, X. Wu, S. Sanvito, Z. Xia, L. Li, Z. Wang, and F. Xiu, Zeeman splitting and dynamical mass generation in Dirac semimetal ZrTe5, Nature Commun. 7, 10.1038/ncomms1251 (2016).
  * Galeski  _et al._ (2020) S. Galeski, X. Zhao, R. Wawrzyńczak, T. Meng, T. Förster, P. M. Lozano, S. Honnali, N. Lmaba, T. Ehmcke, A. Markou, Q. Li, G. Gu, W. Zhu, J. Wosnitza, C. Felser, G. F. Chen, and J. Gooth, Unconventional Hall response in the quantum limit of HfTe5, Nature Commun. 11, 10.1038/s41467-020-19773-y (2020).
  * Zhang  _et al._ (2020) W. Zhang, P. Wang, G. Gu, X. Wu, and L. Zhang, Negative longitudinal magnetothermopower in the topological semimetal ${\mathrm{ZrTe}}_{5}$, Phys. Rev. B 102, 115147 (2020).
  * Jo  _et al._ (2023) N. H. Jo, O. A. Ashour, Z. Shu, C. Jozwiak, A. Bostwick, S. H. Ryu, K. Sun, T. Kong, S. M. Griffin, and E. Rotenberg, On the effects of strain, defects, and interactions on the topological properties of HfTe5 (2023), arXiv:2303.10836 [cond-mat.mtrl-sci] .
  * dos Reis  _et al._ (2016) R. D. dos Reis, M. O. Ajeesh, N. Kumar, F. Arnold, C. Shekhar, M. Naumann, M. Schmidt, M. Nicklas, and E. Hassinger, On the search for the chiral anomaly in Weyl semimetals: the negative longitudinal magnetoresistance, New Journal of Physics 18, 085006 (2016).
  * Li  _et al._ (2017) Y. Li, Z. Wang, P. Li, X. Yang, Z. Shen, F. Sheng, X. Li, Y. Lu, Y. Zheng, and Z.-A. Xu, Negative magnetoresistance in Weyl semimetals NbAs and NbP: Intrinsic chiral anomaly and extrinsic effects, Front. Phys. 12, 10.1007/s11467-016-0636-8 (2017).
  * Weng  _et al._ (2014) H. Weng, X. Dai, and Z. Fang, Transition-Metal Pentatelluride $\mathrm{ZrTe}{}_{5}$ and $\mathrm{HfTe}{}_{5}$: A Paradigm for Large-Gap Quantum Spin Hall Insulators, Phys. Rev. X 4, 011002 (2014).
  * Chen  _et al._ (2015b) R. Y. Chen, Z. G. Chen, X.-Y. Song, J. A. Schneeloch, G. D. Gu, F. Wang, and N. L. Wang, Magnetoinfrared Spectroscopy of Landau Levels and Zeeman Splitting of Three-Dimensional Massless Dirac Fermions in ${\mathrm{ZrTe}}_{5}$, Phys. Rev. Lett. 115, 176404 (2015b).
  * Li  _et al._ (2016b) X.-B. Li, W.-K. Huang, Y.-Y. Lv, K.-W. Zhang, C.-L. Yang, B.-B. Zhang, Y. B. Chen, S.-H. Yao, J. Zhou, M.-H. Lu, L. Sheng, S.-C. Li, J.-F. Jia, Q.-K. Xue, Y.-F. Chen, and D.-Y. Xing, Experimental Observation of Topological Edge States at the Surface Step Edge of the Topological Insulator ${\mathrm{ZrTe}}_{5}$, Phys. Rev. Lett. 116, 176803 (2016b).
  * Wu  _et al._ (2016) R. Wu, J.-Z. Ma, S.-M. Nie, L.-X. Zhao, X. Huang, J.-X. Yin, B.-B. Fu, P. Richard, G.-F. Chen, Z. Fang, X. Dai, H.-M. Weng, T. Qian, H. Ding, and S. H. Pan, Evidence for Topological Edge States in a Large Energy Gap near the Step Edges on the Surface of ${\mathrm{ZrTe}}_{5}$, Phys. Rev. X 6, 021017 (2016).
  * Moreschini  _et al._ (2016) L. Moreschini, J. C. Johannsen, H. Berger, J. Denlinger, C. Jozwiak, E. Rotenberg, K. S. Kim, A. Bostwick, and M. Grioni, Nature and topology of the low-energy states in ${\mathrm{ZrTe}}_{5}$, Phys. Rev. B 94, 081101 (2016).
  * Manzoni  _et al._ (2016) G. Manzoni, L. Gragnaniello, G. Autès, T. Kuhn, A. Sterzi, F. Cilento, M. Zacchigna, V. Enenkel, I. Vobornik, L. Barba, F. Bisti, P. Bugnon, A. Magrez, V. N. Strocov, H. Berger, O. V. Yazyev, M. Fonin, F. Parmigiani, and A. Crepaldi, Evidence for a Strong Topological Insulator Phase in ${\mathrm{ZrTe}}_{5}$, Phys. Rev. Lett. 117, 237601 (2016).
  * Fan  _et al._ (2017) Z. Fan, Q.-L. Liang, Y. B. Chen, S.-H. Yao, and J. Zhou, Transition between strong and weak topological insulator in ZrTe5 and HfTe5, Sci. Rep. 7, 10.1038/srep45667 (2017).
  * Zhang  _et al._ (2017) Y. Zhang, C. Wang, L. Yu, G. Liu, A. Liang, J. Huang, S. Nie, X. Sun, Y. Zhang, B. Shen, J. Liu, W. Hongming, L. Zhao, G. Chen, X. Jia, C. Hu, Y. Ding, W. Zhao, Q. Gao, C. Li, S. He, L. Zhao, F. Zhang, S. Zhang, F. Yang, Z. Wang, Q. Peng, X. Dai, Z. Fang, Z. Xu, C. Chen, and X. J. Zhou, Electronic evidence of temperature-induced Lifshitz transition and topological nature in ZrTe5, Nature Commun. 8, 10.1038/ncomms15512 (2017).
  * Monserrat and Narayan (2019) B. Monserrat and A. Narayan, Unraveling the topology of ${\mathrm{ZrTe}}_{5}$ by changing temperature, Phys. Rev. Res. 1, 033181 (2019).
  * Shahi  _et al._ (2018) P. Shahi, D. J. Singh, J. P. Sun, L. X. Zhao, G. F. Chen, Y. Y. Lv, J. Li, J.-Q. Yan, D. G. Mandrus, and J.-G. Cheng, Bipolar Conduction as the Possible Origin of the Electronic Transition in Pentatellurides: Metallic vs Semiconducting Behavior, Phys. Rev. X 8, 021055 (2018).
  * Jiang  _et al._ (2023) Y. Jiang, T. Zhao, L. Zhang, Q. Chen, H. Zhou, M. Ozerov, D. Smirnov, and Z. Jiang, Revealing temperature evolution of the Dirac band in ${\mathrm{ZrTe}}_{5}$ via magnetoinfrared spectroscopy, Phys. Rev. B 108, L041202 (2023).
  * Fjellvåg and Kjekshus (1986) H. Fjellvåg and A. Kjekshus, Structural properties of ZrTe5 and HfTe5 as seen by powder diffraction, Solid State Communications 60, 91 (1986).
  * Furuseth  _et al._ (1973) S. Furuseth, L. Brattås, and A. Kjekshus, The Crystal Structure of HfTe5, Acta Chem. Scand. 27, 2367 (1973).
  * Kresse and Joubert (1999) G. Kresse and D. Joubert, From ultrasoft pseudopotentials to the projector augmented-wave method, Phys. Rev. B 59, 1758 (1999).
  * Kresse and Furthmüller (1996a) G. Kresse and J. Furthmüller, Efficiency of ab-initio total energy calculations for metals and semiconductors using a plane-wave basis set, Computational Materials Science 6, 15 (1996a).
  * Kresse and Furthmüller (1996b) G. Kresse and J. Furthmüller, Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set, Phys. Rev. B 54, 11169 (1996b).
  * Perdew  _et al._ (1996) J. P. Perdew, K. Burke, and M. Ernzerhof, Generalized Gradient Approximation Made Simple, Phys. Rev. Lett. 77, 3865 (1996).
  * Grimme  _et al._ (2010) S. Grimme, J. Antony, S. Ehrlich, and H. Krieg, A consistent and accurate ab initio parametrization of density functional dispersion correction (DFT-D) for the 94 elements H-Pu, The Journal of Chemical Physics 132, 154104 (2010), https://pubs.aip.org/aip/jcp/article-pdf/doi/10.1063/1.3382344/15684000/154104_1_online.pdf .
  * Adenis  _et al._ (1989) C. Adenis, V. Langer, and O. Lindqvist, Reinvestigation of the structure of tellurium, Acta Crystallographica Section C 45, 941 (1989).
  * Van de Walle and Neugebauer (2004) C. G. Van de Walle and J. Neugebauer, First-principles calculations for defects and impurities: Applications to III-nitrides, Journal of Applied Physics 95, 3851 (2004).