Large Language Model Enhanced Machine Learning Estimators for Classification

Yuhang Wu1, Yingfei Wang2, Chu Wang3, and Zeyu Zheng1 1Department of
Industrial Engineering and Operations Research, University of California
Berkeley
2Foster School of Business, University of Washington
3Amazon

Pre-trained large language models (LLM) have emerged as a powerful tool for
simulating various scenarios and generating output given specific instructions
and multimodal input. In this work, we analyze the specific use of LLM to
enhance a classical supervised machine learning method for classification
problems. We propose a few approaches to integrate LLM into a classical
machine learning estimator to further enhance the prediction performance. We
examine the performance of the proposed approaches through both standard
supervised learning binary classification tasks, and a transfer learning task
where the test data observe distribution changes compared to the training
data. Numerical experiments using four publicly available datasets are
conducted and suggest that using LLM to enhance classical machine learning
estimators can provide significant improvement on prediction performance.

##  1 Introduction

Classification is a fundamental task in supervised machine learning, common
across a wide array of applications. It involves training a model on a dataset
where each instance is assigned a specific class label. When presented with a
new, unlabeled instance, the trained model is expected to accurately predict
the instance’s class.

To illustrate, consider a classification problem aimed at predicting whether a
customer will find a product relevant online. In this problem, each instance
involves two-fold information on (i) customer needs, e.g., reflected by the
search queries of the customer and other attributes, and (ii) product
attributes, including title, description, and/or images of the product. In the
training dataset, each instance is labelled with one of the two classifying
categories: “relevant” or “not relevant”.

There are a wide range of well-established machine learning methods for such
classification tasks, ranging from logistic regression, tree-based methods, to
neural-networks, among many other; see Efron and Hastie (2021). These methods
are typically trained and calibrated using the training dataset. Once trained,
they serve as estimators to predict the labels of new, unlabeled instances.

On the other hand, the emergence of pre-trained large language models (LLMs)
offers an additional approach to these classification tasks, capable of
functioning as estimators with or without additional fine-tuning on task-
specific data. For instance, GPT Radford et al. (2018) by OpenAI can serve as
an estimator by simply taking an instance (comprising customer query and
product information) as input and generating a prediction on whether the
product is relevant to the customer’s needs. Moreover, prompt engineering and
fine-tuning can further enhance the model’s performance on specific
classification tasks.

This work is motivated by the following questions. Can we integrate LLM to
classical machine learning methods to significantly enhance the performance,
compared to the separate use of LLM or the separate use of a classical machine
learning model? What are the different ways of doing such integration, and how
do they perform compared to benchmarks?

We deliver the following results and analysis in this work.

  1. 1.

We analyze the linear combination of an LLM model and a machine learning (ML)
model. We observe that generally LLM predictions are more reliable than that
of a machine learning model on borderline data. We then develop an adaptive
weighted linear combination of LLM and ML to further enhance the performance
via the heavier use of LLM on those regions where ML shows less confidence.

  2. 2.

By treating LLM predictions as additional group information, we apply model
calibration methods to classical machine learning models. This method is
straightforward to carry out and can be used upon any classical machine
learning model.

  3. 3.

We consider the integration of LLM and ML on a transfer learning
classification tasks. For the transfer learning task with covariate shift, we
improve the machine learning model by augmenting the training data with
additional samples from the target distribution, where the labels of those
samples are generated by LLM. We then train a machine learning model on the
augmented dataset to obtain better performance on the target distribution.

  4. 4.

We illustrate the empirical performances of our methods on four public
datasets, including tasks such as relevance prediction, emotion recognition
and hate speech detection. Numerical results show that all our methods perform
better than only using LLM or only using a classical machine learning model.

We would like to add some discussions before proceeding to the main sections
of this work. The advantage of leveraging LLMs to enhance classical ML methods
comes from two primary sources. Firstly, LLMs can serve as a variance
reduction tool in addition to a classical machine learning model trained on
the dataset. The use of LLMs in this context draws a close analogy to the
method of control variates as a variance reduction tool in simulation
literature; refer to Asmussen and Glynn (2007). We also note that the use of
LLM in a classical data-driven method creates additional needs for input data
uncertainty analysis, and leave that for future discussions; see Song et al.
(2014), Feng and Song (2019). Secondly, LLMs can enhance model accuracy by
leveraging their knowledge on a broader range of data; see Møller et al.
(2023), Gao et al. (2023), Chen et al. (2024). This aspect is particularly
beneficial for improving classical ML models in transfer learning tasks.

##  2 Problem formulation and notation

In this work, we focus on binary classification problems, using the task of
predicting relevance as a representative example to illustrate our approach.
The mathematical formulation for this task is consistent with other tasks, and
we will further discuss in the numerical experiments section. The
classification task is formulated as follows. The training data are given by
$\left\\{(\text{query}_{i},\text{product}_{i},\text{label}_{i})\right\\}_{i=1}^{n}$.
Here, $\text{query}_{i}$ consists of the searching content of a customer, e.g.
“modern outdoor furniture”. Next, $\text{product}_{i}$ consists of a set of
product information that includes product description, product features and
possibly product image. The $\text{label}_{i}$ denoted the true relevance
label of “relevant” or “irrelevant”, usually obtained by manual annotation
validated by several independent human annotators. For classification tasks,
embeddings are utilized to derive embedded vectors of $\text{query}_{i}$ and
$\text{product}_{i}$, denoted by
$v_{i1}\in\mathbb{R}^{d_{1}},v_{i2}\in\mathbb{R}^{d_{2}}$. The feature
covariate $x_{i}$ is then given by concatenating $v_{i1}$ and $v_{i2}$ as
$x_{i}=\mathrm{concat}(v_{i1},v_{i2})\in\mathbb{R}^{d}$, where
$d=d_{1}+d_{2}$. In text-based classification, it’s also common to concatenate
text strings using a special token before vectorization. We do not consider
the fine-tuning of LLM at this stage, and simply consider the integration of a
simple machine learning model and a pre-trained LLM model without fine tuning.

The training dataset is then given by
$\mathcal{D}_{\mathrm{train}}=\left\\{(x_{i},y_{i})\right\\}_{i=1}^{n}$, where
$x_{i}\,(1\leq i\leq n)\in\mathbb{R}^{d}$ are features and $y_{i}\,(1\leq
i\leq n)\in\\{0,1\\}$ are the relevance label. Here it is set that $y_{i}=0$
if $\text{label}_{i}$ is relevant, and $y_{i}=1$ for the irrelevant label.
This formulation is generic for a wide range of binary classification
problems.

One can train a classical machine learning model $\hat{f}_{n}(\cdot)$ on
$\mathcal{D}_{\mathrm{train}}$, where $\hat{f}_{n}:\mathbb{R}^{d}\to[0,1]$
maps the feature $x_{i}$ to a score between $0$ and $1$, where the score
represents the chance of irrelevance. If $\hat{y}_{i}=\hat{f}_{n}(x_{i})\leq
0.5$, one predicts “relevant” and otherwise “irrelevant”. The test set is
given by $\mathcal{D}_{\mathrm{test}}=\\{(x_{i},y_{i})\\}_{i=n+1}^{n+m}$ and a
trained machine learning model is evaluated by
${L}_{\mathrm{test}}=\frac{1}{m}\sum_{i=n+1}^{n+m}l\left(y_{i},\hat{y}_{i}\right)$
for some loss function $l(\cdot,\cdot)$. In addition to the classical machine
learning model, in this work we also use a pre-trained large language model
(LLM) to predict $y_{i}$. The LLM is given the input of a query-product pair
and asked to output a score between 0 and 1 representing how likely the LLM
evaluates the query and product as irrelevant. We denote the output of LLM for
the $i$-th pair as $z_{i}\in[0,1]$ for $i=1,\cdots,n+m$.

##  3 Integrating LLM into a Classical ML Estimator

One naive way to utilize LLM is to treat it as a trained model, so
$z_{i}\,(1\leq i\leq n+m)$ are just estimators of $y_{i}\,(1\leq i\leq n+m)$,
in addition to the classical trained machine learning (ML) model, e.g.,
logistic regression. Then a first thought is to combine the two models (LLM
and classical ML) to create a better ensemble model (noted as LLM-ML in this
work). This is a standard combination task and can be executed by several
methods. We illustrate such combination task by examining a simple linear
combination of ML and LLM estimators. We first discuss a straightforward
approach using a fixed constant weight. We then demonstrate through empirical
data analysis that the performance can be enhanced by employing an adaptive
weighting strategy, which adjusts the weights according to the prediction
scores generated by the ML estimator. Such adaptive weighting strategy carries
the thought that LLM may be more heavily used for instances where the
classical ML is not that certain of.

###  3.1 Linear combination of ML and LLM estimators

To combine the ML estimator and LLM estimator, one simple method is to
consider the linear combination with constant weight, i.e.,

$\hat{y}^{\mathrm{Linear}}_{i}=\alpha\hat{y}_{i}+(1-\alpha)z_{i},$ (1)

with some weight $\alpha\in[0,1]$.

The choice of $\alpha$ can be determined through standard cross validation
procedure as follows. By randomly splitting $\mathcal{D}_{\mathrm{train}}$ as
$\cup_{i=j}^{k}\mathcal{D}_{\mathrm{train}}^{j}$ for some
$k\in\mathbb{N}_{+}$, we can train $k$ machine learning models
$\hat{f}^{j}_{n}(\cdot)$ for $j=1,\cdots,k$, where $\hat{f}^{j}_{n}(\cdot)$ is
trained on $\cup_{i=1,i\neq j}^{k}\mathcal{D}_{\mathrm{train}}^{i}$. Then the
best weight $\alpha^{*}$ is given by

$\hat{\alpha}_{n}=\arg\min_{\alpha\in[0,1]}\frac{1}{n}\sum_{j=1}^{k}\sum_{(x_{i},y_{i})\in\mathcal{D}_{\mathrm{train}}^{j}}l\left(\alpha\hat{f}^{j}_{n}(x_{i})+(1-\alpha)z_{i},y_{i}\right).$
(2)

If we omit the superscript $j$ and denote $\hat{f}^{j}_{n}(x_{i})$ by
$\hat{y}_{i}^{\mathrm{cv}}$, for the $l_{2}$ loss $l(x,y)=(x-y)^{2}$, this
procedure is equivalent to regress $\\{y_{i}-z_{i}\\}_{i=1}^{n}$ on
$\\{\hat{y}_{i}^{\mathrm{cv}}-z_{i}\\}_{i=1}^{n}$, which gives
$\hat{\alpha}_{n}^{l_{2}}=\frac{\sum_{i=1}^{n}(y_{i}-z_{i})(\hat{y}_{i}^{\mathrm{cv}}-z_{i})}{(\hat{y}_{i}^{\mathrm{cv}}-z_{i})^{2}}$.
Informally speaking, if we assume
$(x_{i},z_{i},y_{i})\stackrel{{\scriptstyle\mathrm{i.i.d.}}}{{\sim}}P$ for
some distribution $P$ on $\mathbb{R}^{d}\times[0,1]^{2}$, and assume that
$\hat{f}_{n}$ converges to some $f_{0}$ as $n\to+\infty$, then we will have

$\hat{\alpha}_{n,l_{2}}\to\alpha^{*}=\frac{\mathbb{E}_{P}\left[(Y-Z)(f_{0}(X)-Z)\right]}{\mathbb{E}_{P}\left(f_{0}(X)-Z\right)^{2}},\text{
as }n\to+\infty,$

which is the minimizer of the mean-squared error
$\mathrm{MSE}_{\mathrm{L}}(\alpha)=\mathbb{E}\left[\alpha
f_{0}(X)+(1-\alpha)Z-Y\right]^{2},$ so
$\mathrm{MSE}_{\mathrm{L}}(\alpha^{*})\leq\mathrm{MSE}_{\mathrm{L}}(0)=\mathbb{E}(Z-Y)^{2}$
and
$\mathrm{MSE}_{\mathrm{L}}(\alpha^{*})\leq\mathrm{MSE}_{\mathrm{L}}(1)=\mathbb{E}(f_{0}(X)-Y)^{2}$,
which means the asymptotic behaviour of the LLM-ML ensemble estimator is no
worse than two baseline estimators.

###  3.2 Linear combination with adaptive weights

One can improve this method by taking $\alpha$ to be dependent on
$\hat{y}_{i}$ as $\alpha(\hat{y}_{i})$. Consider the following visualization
results in Figure 1. Here t-SNE Hinton and Roweis (2002) is a statistical
method for visualizing high-dimensional data by giving each sample a location
in a two-dimensional map. We can see the classical machine learning model
sometimes does not do well on the boundary of relevant samples and irrelevant
samples, and many of such errors are made when the prediction score
$\hat{y}_{i}$ is close to $0.5$. On the contrary, the LLM model shows
different patterns, in the sense that its accuracy is relatively stable even
for those borderline data. We refer to the red circles in the figure for some
illustration. This motivates us to take $\alpha(\hat{y}_{i})$ to be larger
when $\hat{y}_{i}$ is close to $0$ or $1$, and if $\hat{y}_{i}$ is close to
$0.5$, which indicates $x_{i}$ is likely to located on the boundary, we set it
to be smaller.

Figure 1: t-SNE visualization of an illustration data set

In this case, the LLM-ML estimator is given as an adaptive weighted linear
combination of two estimator:

$\hat{y}^{\mathrm{AL}}_{i}=\alpha(\hat{y}_{i})\cdot\hat{y}_{i}+\left(1-\alpha(\hat{y}_{i})\right)\cdot
z_{i}.$ (3)

We would expect $\alpha(\cdot):[0,1]\to[0,1]$ to be a decreasing function on
$[0,0.5]$, and then increases on $[0.5,1]$. By taking a predetermined
hypothesis space $\mathcal{A}\subset\\{\alpha|\alpha:[0,1]\to[0,1]\\}$ for all
possible choices of $\alpha(\cdot)$ we are interested, the choice of
$\alpha(\cdot)$ can similarly be determined through cross validation as what
we just did for a fixed $\alpha$. To see this, we can now modify (2) as:

$\hat{\alpha}_{n}=\arg\min_{\alpha\in\mathcal{A}}\frac{1}{n}\sum_{i=1}^{n}l\left(\alpha\left(\hat{y}_{i}^{\mathrm{cv}}\right)\cdot\hat{y}_{i}^{\mathrm{cv}}+\left(1-\alpha\left(\hat{y}_{i}^{\mathrm{cv}}\right)\right)\cdot
z_{i},y_{i}\right),$ (4)

where $\hat{y}_{i}^{\mathrm{cv}}=\hat{f}^{j}_{n}(x_{i})$ and by slightly
abusing the notation $\hat{\alpha}_{n}$ is now a function in $\mathcal{A}$.

We now consider the choice of $\mathcal{A}$. A large $\mathcal{A}$ such as all
continuous functions from $[0,1]$ to $[0,1]$
$\\{\alpha|\alpha:[0,1]\to[0,1],\alpha\in C\\}$ will make the optimization
problem in (4) intractable and can easily lead to overfit, but a simple
$\mathcal{A}$ such as all constant functions on $[0,1]$ will just give a fixed
weight and cannot fully use the strength of LLM. In this work, we consider the
class of piecewise constant functions, i.e.

$\mathcal{A}_{r}=\\{\alpha|\alpha(x)=\sum_{i=1}^{r}w_{i}\mathbf{1}_{\\{x\in
I_{i}\\}}\text{ with }w_{i}\in[0,1](i=1,\cdots,r)\\},$ (5)

here $r\in\mathbb{N}_{+}$ is a pre-determined positive integer and
$I_{i}(1\leq i\leq r)$ is a partition of $[0,1]$. For simplicity, we take
$I_{i}=[\frac{i-1}{r},\frac{i}{r})$ for $1\leq i\leq r-1$ and
$I_{r}=[\frac{r-1}{r},1]$. In this case, (4) can now be formulated as:

$\displaystyle\hat{\alpha}_{n,r}$
$\displaystyle=\arg\min_{\alpha\in\mathcal{A}_{r}}\frac{1}{n}\sum_{i=1}^{n}l\left(\alpha\left(\hat{y}_{i}^{\mathrm{cv}}\right)\cdot\hat{y}_{i}^{\mathrm{cv}}+\left(1-\alpha\left(\hat{y}_{i}^{\mathrm{cv}}\right)\right)\cdot
z_{i},y_{i}\right)$ (6)
$\displaystyle=\arg\min_{w_{1},\cdots,w_{r}}\frac{1}{n}\sum_{j=1}^{r}\sum_{i:\hat{y}_{i}^{\mathrm{cv}}\in
I_{j}}l\left(w_{j}\hat{y}_{i}^{\mathrm{cv}}+\left(1-w_{j}\right)z_{i},y_{i}\right).$

This is equivalent to solve $r$ optimization problems for samples with indices
in $\\{i|\hat{y}_{i}^{\mathrm{cv}}\in I_{j}\\}(j=1,\cdots,r)$ separately, and
each of them is what we just discussed in last part. If we are using $l_{2}$
loss then clearly this piecewise constant weight function gives an MSE no
larger than any fixed constant weight.

Algorithm 1 Adaptive weighted linear combination of LLM and ML

1: Input:  $\mathcal{D}_{\mathrm{train}}=\\{(x_{i},z_{i},y_{i})\\}_{i=1}^{n}$,
Number of pieces $r$, and a partition of $[0,1]$ $I_{i}(1\leq i\leq r)$;

2: Train a ML model $f(\cdot)$ on $\\{(x_{i},y_{i})\\}_{i=1}^{n}$ and
construct cross validation predictions
$\\{(x_{i},\hat{y}_{i}^{\mathrm{cv}})\\}_{i=1}^{n}$.

3: Solve subproblems $w_{j}=\arg\min_{w}\sum_{i:\hat{y}_{i}^{\mathrm{cv}}\in
I_{j}}l\left(w\hat{y}_{i}^{\mathrm{cv}}+\left(1-w\right)z_{i},y_{i}\right)$
for $j=1,\cdots,r$.

4: Output: the weight function
$\alpha(x)=\sum_{i=j}^{r}w_{j}\mathbf{1}_{\\{x\in I_{j}\\}}$.

##  4 Calibration with LLM

In this section, we consider calibrating a classical ML model (e.g., logistic
regression) with LLM. We start by introducing some backgrounds on calibration,
then we discuss how to use LLM to calibrate a classical ML model.

###  4.1 Background on calibration

Suppose our data are realizations of $(X,Y)\sim\mathcal{D}$ with
$Y\in\\{0,1\\}$, and the model is $f$. Ideally, we hope to find a model such
that

$f(x)=P(Y=1|X=x).\text{ (Informally)}$ (7)

However, this is generally unrealistic and there are at least two reasons.
First, for each fixed $x$, the label is already determined, and there is no
randomness so the probability in the RHS is not well-defined. Second, the
domain of $X$ in our scenario consists of pairs of query and item, and each of
them comes from an extremely large set of English words, so it is generally
impossible to find such a function with limited samples. Calibration is a
tractable condition that can be viewed as a coarsening of the condition above.

###### Definition 4.1

(Calibration) For $(X,Y)\sim\mathcal{D}$, the bias of $f$ at the $p$-th level
set is defined as

$\Delta_{p}(f)\triangleq\mathbb{E}_{\mathcal{D}}[Y-f(X)|f(X)=p].$ (8)

If $\Delta_{p}(f)\equiv 0$, we say that $f$ is calibrated w.r.t.
$\mathcal{D}$. We omit the subscript $\mathcal{D}$ when there is no ambiguity.

Intuitively, calibration can be considered as a minimal condition for a model
$f$ to be good, in the sense that the probability of $Y=1$ conditional on
$f(x)=p$ is indeed $p$. For each model $f$, if it is not calibrated, we can
define $\hat{f}(x)=f(x)+\Delta_{f(x)}(f)$ and it is straightforward to verify
that $\hat{f}$ is calibrated. However, it is infeasible to conditional on
$f(X)=p$ in practice with finite samples, and a simple strategy is to do
discretization on a uniform grid and calibrate $f$ on the grid. Specifically,
we define the uniform grid $[\frac{1}{M}]=\\{\frac{i}{m}\\}_{i=0}^{M}$ for
some $M\in\mathbb{N}_{+}$, and define

$f^{\prime}(x)=\arg\min_{u\in[\frac{1}{M}]}|f(x)-u|.$ (9)

This rounds the output of $f$ to the grid $[\frac{1}{M}]$ and can be viewed as
a discretization of $[0,1]$. Now, for our data
$\\{(x_{i},y_{i})\\}_{i=1}^{n}$, we can construct the calibration of $f$ as

$\hat{f}^{(M)}(x)=f(x)+\hat{\Delta}_{f^{\prime}(x)}(f),$ (10)

where $\hat{\Delta}_{f^{\prime}(x)}(f)$ can be estimated through data by
calculating the mean of $y_{i}-f(x_{i})$ conditional on the value of
$f^{\prime}(x_{i})$. We refer to Roth (2022) for more theoretical results of
$\hat{f}^{(M)}$.

###  4.2 Calibrate ML model with LLM

Above calibration procedure only involves the model $f$, and we now integrate
LLM into this procedure. In this scenario, we assume our data are realizations
of $(X,Z,Y)\sim P$, where $Z\in[0,1]$ is the output of LLM. We have a ML model
$f$ and we want to use it to construct a improved model, which takes $X$ and
$Z$ as inputs and generates a prediction for $Y$. We start by introducing an
enhanced version of Definition 4.1, which is also known as the “multi-
accuracy” condition in Kim et al. (2019).

###### Definition 4.2

(Multi-accuracy with LLM) Suppose $(X,Z,Y)\sim P$, A model $f$ is multi-
accurate w.r.t. $Z$ if

$\mathbb{E}_{P}\left[Y-f(X)|Z\right]=0.$ (11)

If we treat $Z$ as a covariate, then Definition 4.2 can be viewed as a
condition for a model $f$ to be consistent conditional on the covariate. Now,
with both Definition 4.1 and 4.2, we want to construct some
$\hat{f}^{M,M^{\prime}}$ from $f$ such that $\hat{f}^{M,M^{\prime}}$ is
calibrated and multi-accurate. Again it is infeasible to conditional on $Z$ if
$Z$ is continuously distributed with finite samples, so we consider a
discretization similar as what we did in (9). Let
$[\frac{1}{M^{\prime}}]=\\{\frac{i}{m}\\}_{i=0}^{M}$ for some
$M^{\prime}\in\mathbb{N}_{+}$, and define

$S_{p,q}(f)=\left\\{(x,z)\Big{|}\arg\min_{u\in[\frac{1}{M}]}|u-f(x)|=p,\arg\min_{v\in[\frac{1}{M^{\prime}}]}|v-z|=q\right\\}$
(12)

for $p=\frac{i}{M}(0\leq i\leq M)$ and $q=\frac{j}{M^{\prime}}(0\leq j\leq
M^{\prime})$. Also define

$\Delta_{p,q}(\hat{f},f)\triangleq\mathbb{E}\left[Y-\hat{f}(X)|S_{p,q}(f)\right].$
(13)

Then one simple way to construct some $\hat{f}^{M,M^{\prime}}$ is to make sure
$\Delta_{p,q}(\hat{f}^{M,M^{\prime}},f)=0$ for all $p,q$. Informally speaking,
if $\Delta_{p,q}(\hat{f}^{M,M^{\prime}},f)=0$ for all $p,q$ then
$\hat{f}^{M,M^{\prime}}$ is calibrated and multi-accurate up to some
discretization error. This naturally leads to the following choice:

$\hat{f}^{M,M^{\prime}}_{1}(x)=f(x)+\hat{\Delta}_{p,q}(f,f),\text{ if
}(x,z)\in S_{p,q}(f),$ (14)

where

$\hat{\Delta}_{p,q}(f,f)=\frac{\sum_{1\leq i\leq n,(x_{i},z_{i})\in
S_{p,q}(f)}\left(y_{i}-f(x_{i})\right)}{\sum_{1\leq i\leq n,(x_{i},z_{i})\in
S_{p,q}(f)}1}$ (15)

is the empirical mean of $Y-f(X)$ conditional on $S_{p,q}(f)$ with samples in
the training set. Then, if we assume
$\\{(x_{i},z_{i},y_{i})\\}_{i=1}^{n+m}\stackrel{{\scriptstyle\text{i.i.d.}}}{{\sim}}P$,
we have

$\displaystyle\Delta_{p,q}(\hat{f}^{M,M^{\prime}}_{1},f)$ (16)
$\displaystyle=$
$\displaystyle\mathbb{E}\left[Y-\hat{f}^{M,M^{\prime}}(X)\big{|}S_{p,q}(f)\right]$
$\displaystyle=$
$\displaystyle\mathbb{E}\left[Y-f(X)|S_{p,q}(f)\right]-\mathbb{E}\left[\hat{\Delta}_{p,q}(f,f)|S_{p,q}(f)\right]$
$\displaystyle=$ $\displaystyle 0,$

which is exactly what we want.

While (14) gives one way to construct a function that is both calibrated and
multi-accurate up to some discretization error, the problem is that it
involves calculating $(M+1)\times(M^{\prime}+1)$ conditional means for
$S_{p,q}(f)$. When $M$ and $M^{\prime}$ are large, this quantity may be of the
same order of the sample size $n$ or even larger, making the estimation of
conditional means not precise enough. This motivates the following choice. We
take

$\hat{f}^{M,M^{\prime}}_{2}(x)=f(x)+\hat{w}_{i,0}+\hat{w}_{j,1},\text{ if
}(x,z)\in S_{\frac{i}{M},\frac{j}{M^{\prime}}}(f),$ (17)

where $\hat{w}_{i,0}(0\leq i\leq M)$ and $\hat{w}_{j,1}(0\leq j\leq
M^{\prime})$ are determined by the following optimization problem:

$\hat{w}_{0,0},\cdots,\hat{w}_{M,0},\hat{w}_{0,1},\cdots,\hat{w}_{M^{\prime},1}=\arg\min_{\begin{subarray}{c}w_{0,0},\cdots,w_{M,0},\\\
w_{0,1},\cdots,w_{M^{\prime},1}\end{subarray}}\sum_{i=1}^{n}\left(\hat{f}^{M,M^{\prime}}_{2}(x_{i})-y_{i}\right)^{2}.$
(18)

Since this is a least square linear regression problem, the solution always
exists. For any solution, the optimality conditions give:

$\sum_{f^{\prime}(x_{i})=\frac{j}{M}}\left(\hat{f}^{M,M^{\prime}}_{2}(x_{i})-y_{i}\right)=0,\quad\forall~{}0\leq
j\leq M,$ (19)

so
$\mathbb{E}\left[Y-\hat{f}^{M,M^{\prime}}_{2}(X)\Big{|}f^{\prime}(X)\right]=0$,
which implies $\hat{f}^{M,M^{\prime}}_{2}$ is calibrated up to some
discretization error. Similarly, it is also multi-accurate up to some
discretization error. The construction of $\hat{f}^{M,M^{\prime}}_{2}$ only
requires the estimation of $M+M^{\prime}+2$ parameters, which is much smaller
than $(M+1)\times(M^{\prime}+1)$ parameters for $\hat{f}^{M,M^{\prime}}_{1}$,
and can thus mitigate issues such as large estimation error and overfitting in
$\hat{f}^{M,M^{\prime}}_{1}$.

Algorithm 2 Calibrate ML model with LLM

1: Input:  $\mathcal{D}_{\mathrm{train}}=\\{(x_{i},z_{i},y_{i})\\}_{i=1}^{n}$,
$M,M^{\prime}\in\mathbb{N}_{+}$;

2: Train a ML model $f(\cdot)$ on $\\{(x_{i},y_{i})\\}_{i=1}^{n}$.

3: Output: the calibration of $f(\cdot)$ given by
$\hat{f}_{1}^{M,M^{\prime}}(\cdot)$ in (14) or
$\hat{f}_{2}^{M,M^{\prime}}(\cdot)$ in (17).

##  5 Transfer learning with LLM

Recall that in Figure 1 we saw that LLM can perform more stable than a
classical ML model, this motivates us to apply LLM to transfer learning, where
the training distribution of the model can be different from that of its
application. Such transfer learning tasks can be challenging to a classical ML
model.

Again consider the relevance label prediction task, suppose that there are no
or very few samples in the labeled training set that contains products related
to “bed”, but the test datasets contain a large number of products related to
“bed”. This can happen in applications where a new type of product is
introduced to the platform. We would expect that the classical ML model may
not do well on the test datasets because the model has not seen much
information about “bed”. One possible remedy is to utilize LLM to augment the
training set for the bed category by adding data labeled with predictions from
the LLM. This can be formulated as follows. Suppose we have a training set
$\mathcal{D}_{\mathrm{train}}=\\{(x_{i},y_{i})\\}_{i=1}^{n}\stackrel{{\scriptstyle\text{i.i.d.}}}{{\sim}}P_{1}$,
and we train a ML model $\hat{f}_{n}\in\mathcal{F}$ on
$\mathcal{D}_{\mathrm{train}}$ for some model class $\mathcal{F}$. We want to
evaluate its performance by its expected loss on $P_{2}$ w.r.t. some loss
function $l$: $\mathbb{E}_{(X,Y)\sim P_{2}}\left[l(\hat{f}_{n}(X),Y)\right]$.
Here we assume $P_{1}=P_{X,1}\times P_{Y|X}$ and $P_{2}=P_{X,2}\times
P_{Y|X}$, where $P_{X,l}~{}(l=1,2)$ are two distributions of $X$ on its domain
$\mathcal{X}$ and the conditional distribution of $Y|X$ is given by $P_{Y|X}$
and is the same for $P_{1}$ and $P_{2}$. This scenario is usually referred as
covariate shift Sugiyama et al. (2007). While we do not have other labeled
training data, we are accessible to sample more $X$ as we want, that is to
say, we can augment the training set by adding
$\\{(x_{i},z_{i})\\}_{i=n+1}^{n+m}$ for some $x_{i}$ and $m\in\mathbb{N}_{+}$
as we want. With the augmented training set, we can train a new ML model
$\hat{f}_{n,m}$ as follows:

$\hat{f}_{n,m}=\arg\min_{f\in\mathcal{F}}\frac{1}{n+m}\sum_{i=1}^{n}l\left(f(x_{i}),y_{i}\right)+\frac{1}{n+m}\sum_{i=n+1}^{n+m}l_{0}(f(x_{i}),z_{i}).$
(20)

Here $l(\cdot,\cdot)$ is the loss we are interested, and $l_{0}(\cdot,\cdot)$
is some weak-supervised loss that can be viewed as a relaxation of
$l(\cdot,\cdot)$. The rational is that the label of
$\\{x_{i}\\}_{i=n+1}^{n+m}$ are generated from LLM and may be biased and
noisy, so a relaxed loss such as
$l_{0}(x,y)=\min_{\varepsilon\in[-a,a]}l\left(x,y+\varepsilon\right)$ for some
$a\geq 0$ gives looser penalization on those samples.

Regarding the choice of $\\{x_{i}\\}_{i=n+1}^{n+m}$, suppose we sample them
from some distribution $P_{X,3}$. Let the probability density function of
$P_{X,3}$ be $p_{i}(x)$ for $i=1,2,3$, then in the ideal case we would expect

$\frac{n}{n+m}p_{1}(x)+\frac{m}{n+m}p_{3}(x)\equiv p_{2}(x),$ (21)

which then gives a choice of $P_{3}$ as

$p_{3}(x)=p_{2}(x)+\frac{n}{m}\left(p_{2}(x)-p_{1}(x)\right).$ (22)

However, in order to make sure $p_{3}(x)\geq 0$, we need
$m\geq\sup_{x}\frac{p_{1}(x)-p_{2}(x)}{p_{2}(x)}n$. When the quantity in the
RHS is large, we will need lots of labels generated from LLM, which may not be
a reasonable choice. In the extreme scenario that $p_{2}(x)=0$ and
$p_{1}(x)>0$, it is impossible to find a $p_{3}(x)\in[0,1]$ satisfying (22).
Thus, in these scenarios, a heuristic choice such as
$p_{3}(x)\propto\max\left(p_{2}(x)+\frac{n}{m}\left(p_{2}(x)-p_{1}(x)\right),0\right)$
or simply $p_{3}(x)=p_{2}(x)$ may be better. In Section 6.3, we illustrate our
method in the scenario that the supports of $p_{1}(\cdot)$ and $p_{2}(\cdot)$
are disjoint, and we show that a naive choice of $p_{3}(x)=p_{2}(x)$ can
already significantly improve the performance of the classical machine
learning model.

Algorithm 3 Transfer learning with LLM

1: Input:  $\mathcal{D}_{\mathrm{train}}=\\{(x_{i},z_{i},y_{i})\\}_{i=1}^{n}$
and a training distribution $P_{1}$, a target distribution $P_{2}$,
$m\in\mathbb{N}_{+}$;

2: Choose a sampling distribution $P_{X,3}$ based on $P_{1}$ and $P_{2}$,
sample $\\{x_{i}\\}_{i=n+1}^{n+m}$.

3: Label $\\{x_{i}\\}_{i=n+1}^{n+m}$ with LLM to obtain
$\\{(x_{i},z_{i})\\}_{i=n+1}^{n+m}$

4: Train a ML model $\hat{f}_{n,m}$ on
$\\{(x_{i},y_{i})\\}_{i=1}^{n}\cup\\{(x_{i},z_{i})\\}_{i=n+1}^{n+m}$.

5: Output: $\hat{f}_{n,m}$.

##  6 Numerical experiments

We now illustrate the performances of the proposed LLM-ML integration methods
through numerical experiments with public datasets. In Section 6.1, we provide
the information of the datasets in use. In Section 6.2, we compare the
performances of the proposed methods in Section 3 and 4. In Section 6.3, we
consider the transfer learning task in Section 5. Our code is available at
https://github.com/wyhArturia/llm_enhanced_ml.

###  6.1 Datasets

We illustrate our methods on four public datasets. For predicting relevance
labels, we use the Wayfair Annotation DataSet (WANDS) Chen et al. (2022).
WANDS is a discriminative, reusable, and fair human-labeled dataset. It is one
of the biggest publicly available search relevance dataset and is effective in
evaluating and discriminating between different models. It consists of $480$
queries, $42994$ products, and $233000$ annotated query-product relevance
labels. We also use the following three public datasets: Yelp’s dataset
(https://www.yelp.com/dataset), Emotion dataset Saravia et al. (2018) and Hate
speech dataset De Gibert et al. (2018). All of them are NLP classification
tasks and have similar formulations as the relevance label prediction task as
described in the manuscript.

###  6.2 Combined estimators and calibration with LLM

We first illustrate the LLM-ML combined estimators in Section 3 and the
calibrated estimators in Section 4 on four datasets. We randomly split each
dataset into a training set and a testing set, and we compare the performances
of following methods on testing sets:

  1. 1.

The large language model method (LLM): directly use GPT-3.5-Turbo-Instruct to
classify the input instance. The prompt for each task consists of step-by-step
instructions such as key indicators, evaluation rules, examples and quality
checks. All prompts are provided in our code.

  2. 2.

The machine learning method (ML): for each training set which consists of word
vectors and labels, we train a logistic regression model and apply it to the
testing set for predictions. We use this simple machine learning model for
illustration purposes, as we are mainly focus on how to use LLM to improve the
ML model instead of the performance of the ML model.

  3. 3.

Naive linear combination of LLM and ML methods (Linear): our method in Section
3.1 that uses a linear combination of LLM and ML estimators with a fixed
weight. The weight is estimated through least square method on the training
set.

  4. 4.

Adaptive weighted linear combination (AdaLinear): our method in Section 3.2
that uses a linear combination of LLM and ML estimators with a piecewise
constant weight function. The weight function is estimated through least
square method on the training set. The number of pieces $r$ is also tuned on
the training set, but its value will not affect the performance too much as
long as it is not too small ($r=1$) or too large ($r>20$).

  5. 5.

Calibrated estimator with LLM (Calibration): our method in Section 4.2 that
calibrates the ML estimator with LLM predictions. Since our LLM predictions
are binary, we take $M^{\prime}=2$. The value of $M$ is tuned on the training
set. We use the construction in (14) for simplicity.

The results are given as follows:

Table 1: Test accuracy of different methods on four datasets | LLM | ML | Linear | AdaLinear$(r)$ | Calibration$(M)$
---|---|---|---|---|---
WANDS | $0.775$ | $0.803$ | $0.838$ | $\mathbf{0.846}(4)$ | $0.840(10)$
Yelp | $0.724$ | $0.691$ | $0.739$ | $0.742(4)$ | $\mathbf{0.743}(20)$
Emotion | $0.759$ | $0.799$ | $0.806$ | $\mathbf{0.812}(10)$ | $0.806(10)$
Hate | $0.672$ | $0.717$ | $0.720$ | $0.724(4)$ | $\mathbf{0.731}(12)$

From the above results we can see that, on all four datasets, our three
methods (Linear, AdaLinear and Calibration) that utilizing both ML and LLM
estimators perform better than only use LLM or ML. The AdaLinear method
performs better than the Linear method, which is reasonable as Linear is a
special case of AdaLinear. AdaLinear and Calibration have comparable
performances, indicating that both model ensemble and calibration are
reasonable ways to enhance ML models with LLM.

###  6.3 Transfer learning with LLM

We now consider the task of transfer learning with LLM in Section 5. In order
to illustrate the performance of our method, we manually construct a table
dataset and a bed dataset from the WANDS dataset and are provided in our code.
We compare the performances of following methods on testing sets:

  1. 1.

The large language model method (LLM): directly use GPT-3.5-Turbo to classify
the input instances, similar as that in Section 6.2.

  2. 2.

The machine learning model for the table dataset (MLTable): we train a
logistic regression model on $1000$ word vectors in the table dataset. We
would expect this ML model works well on the table dataset, but may not do
well on the bed dataset.

  3. 3.

Naive linear combination of LLM and ML methods (Linear): our method in Section
3.1 that uses a linear combination of LLM and MLTable methods with a fixed
weight, similar as that in Section 6.2. We use this as one baseline method.
The reason we do not use an adaptive weight function is that, we can only
calculate the weight function on the training set, which can lead to
overfitting as the target distribution is different from the training set.

  4. 4.

The transfer learning method with LLM discussed in Section 5 (Transfer$(m)$):
we randomly add $m$ bed samples and corresponding labels generated by LLM into
our table dataset, then train a logistic regression model on all $10000+m$
word vectors. With the notations in Section 5, this is equivalent to take
$p_{3}(x)=p_{2}(x)$.

We uses a table testing set and a bed testing set to evaluate these methods.
The results are given as follows:

Table 2: Test accuracy of four methods on two datasets | LLM | MLTable | Linear | Transfer$(5000)$ | Transfer$(10000)$ | Transfer$(15000)$
---|---|---|---|---|---|---
Table | $0.738$ | $0.915$ | $\mathbf{0.923}$ | $0.910$ | $0.910$ | ${0.911}$
Bed | $0.753$ | $0.721$ | $0.737$ | $0.784$ | $\mathbf{0.786}$ | $\mathbf{0.786}$

We can see that the Transfer method always performs better than LLM and
MLTable on the bed dataset with a small sacrifice on its performance on the
table dataset, regardless of the choice of $m$. However, the choice of $m$ is
subtle and can still affect the performance of our method. The linear
combination of LLM and MLTable can perform well on the table dataset, but is
not robust to the covariate shift and is not as good as our Transfer method on
the bed dataset.

##  7 Conclusion

We conclude the work by discussing some unexplored aspects of this work. Our
goal is to show that LLM can enhance a classification machine learning method
for classification problems. We have selected the standard logistic regression
method as the benchmark machine learning method. We believe that the exact
size of improvement brought by LLM will vary by choosing different machine
learning methods as benchmark. Also, the performance of LLM can vary by using
different prompts. We did not push on prompt engineering to figure out the
best way to design prompts to max out the LLM potential. The proposed
approaches in this work may be further expanded by integrating prompt
engineering and fine tuning with certain training data.

## References

  * Asmussen and Glynn (2007) Asmussen S, Glynn PW (2007) _Stochastic simulation: algorithms and analysis_ , volume 57 (Springer).
  * Chen et al. (2022) Chen Y, Liu S, Liu Z, Sun W, Baltrunas L, Schroeder B (2022) Wands: Dataset for product search relevance assessment. _European Conference on Information Retrieval_ , 128–141 (Springer).
  * Chen et al. (2024) Chen Z, Mao H, Li H, Jin W, Wen H, Wei X, Wang S, Yin D, Fan W, Liu H, et al. (2024) Exploring the potential of large language models (llms) in learning on graphs. _ACM SIGKDD Explorations Newsletter_ 25(2):42–61.
  * De Gibert et al. (2018) De Gibert O, Perez N, García-Pablos A, Cuadros M (2018) Hate speech dataset from a white supremacy forum. _arXiv preprint arXiv:1809.04444_ .
  * Efron and Hastie (2021) Efron B, Hastie T (2021) _Computer age statistical inference, student edition: algorithms, evidence, and data science_ , volume 6 (Cambridge University Press).
  * Feng and Song (2019) Feng BM, Song E (2019) Efficient input uncertainty quantification via green simulation using sample path likelihood ratios. _2019 Winter Simulation Conference (WSC)_ , 3693–3704 (IEEE).
  * Gao et al. (2023) Gao Y, Sheng T, Xiang Y, Xiong Y, Wang H, Zhang J (2023) Chat-rec: Towards interactive and explainable llms-augmented recommender system. _arXiv preprint arXiv:2303.14524_ .
  * Hinton and Roweis (2002) Hinton GE, Roweis S (2002) Stochastic neighbor embedding. _Advances in neural information processing systems_ 15\.
  * Kim et al. (2019) Kim MP, Ghorbani A, Zou J (2019) Multiaccuracy: Black-box post-processing for fairness in classification. _Proceedings of the 2019 AAAI/ACM Conference on AI, Ethics, and Society_ , 247–254.
  * Møller et al. (2023) Møller AG, Aarup Dalsgaard J, Pera A, Aiello LM (2023) Is a prompt and a few samples all you need? using gpt-4 for data augmentation in low-resource classification tasks. _arXiv e-prints_ arXiv–2304.
  * Radford et al. (2018) Radford A, Narasimhan K, Salimans T, Sutskever I, et al. (2018) Improving language understanding by generative pre-training .
  * Roth (2022) Roth A (2022) Uncertain: Modern topics in uncertainty estimation.
  * Saravia et al. (2018) Saravia E, Liu HCT, Huang YH, Wu J, Chen YS (2018) CARER: Contextualized affect representations for emotion recognition. _Proceedings of the 2018 Conference on Empirical Methods in Natural Language Processing_ , 3687–3697 (Brussels, Belgium: Association for Computational Linguistics), URL http://dx.doi.org/10.18653/v1/D18-1404.
  * Song et al. (2014) Song E, Nelson BL, Pegden CD (2014) Advanced tutorial: Input uncertainty quantification. _Proceedings of the Winter Simulation Conference 2014_ , 162–176 (IEEE).
  * Sugiyama et al. (2007) Sugiyama M, Krauledat M, Müller KR (2007) Covariate shift adaptation by importance weighted cross validation. _Journal of Machine Learning Research_ 8(5).


# Semiclassical perspective on Landau levels and Hall conductivity in an
anisotropic Cubic Dirac Semi-Metal and the peculiar case of star-shaped
classical orbits

Ahmed Jellal  Laboratory of Theoretical Physics, Faculty of Sciences, Chouaïb
Doukkali University, PO Box 20, 24000 El Jadida, Morocco  Canadian Quantum
Research Center, 204-3002 32 Ave Vernon, BC V1T 2L7, Canada  Hocine Bahlouli
Physics Department, King Fahd University of Petroleum $\&$ Minerals, Dhahran
31261, Saudi Arabia  Michael Vogl  Physics Department, King Fahd University of
Petroleum $\&$ Minerals, Dhahran 31261, Saudi Arabia  Interdisciplinary
Research Center (IRC) for Intelligent Secure Systems, KFUPM, Dhahran, Saudi
Arabia

###### Abstract

We study an anisotropic cubic Dirac semi-metal subjected to a constant
magnetic field. In the case of an isotropic dispersion in the $x$-$y$ plane,
with parameters $v_{x}=v_{y}$, it is possible to find exact Landau levels,
indexed by the quantum number $n$, using the typical ladder operator approach.
Interestingly, we find that the lowest energy level (the zero energy state in
the case $k_{z}=0$) has a degeneracy that is three times that of other states.
This degeneracy manifests in the Hall conductivity as a step at zero chemical
potential that is 3/2 the size of other steps. Moreover, as $n\to\infty$ we
find energies $E_{n}\propto n^{3/2}$, which means the $n$-th step as a
function of chemical potential roughly occurs at a value $\mu\propto n^{3/2}$.
We propose that these exciting features could be used to identify cubic Dirac
semi-metals experimentally. Subsequently, we analyze the anisotropic case
$v_{y}=\lambda v_{x}$ with $\lambda\neq 1$. First, we consider a perturbative
treatment around $\lambda\approx 1$ and find that energies $E_{n}\propto
n^{3/2}$ still holds as $n\to\infty$. To gain further insight into the Landau
level structure for a maximum anisotropy, we turn to a semi-classical
treatment that reveals interesting star-shaped orbits in phase space that
close at infinity. This property is a manifestation of weakly localized
states. Despite being infinite in length, these orbits enclose a finite phase
space volume and permit finding a simple semi-classical formula for the
energy, which again has the form as above. Our findings suggest that both
isotropic and anisotropic cubic Dirac semi-metals should leave similar
experimental imprints.

##  I  Introduction

Semi-classical approximations typically are achieved by taking the
$\hbar\rightarrow 0$ limit of quantum mechanics [1, 2]. More precisely, this
is the limit at which the characteristic action $S$ is large compared to the
reduced Planck constant $\hbar$, that is $S\gg\hbar$ [3, 4]. Within this
limit, capturing much of the rich and exotic behavior of quantum systems and
preserving the intuitive explanations that classical mechanics provides is
possible. More than that, semi-classical physics has allowed us to gain deep
insights into the interplay between classical chaos and quantum mechanics [5].
A prototypical example of this is the effect of quantum scarring that occurs
with wavefunctions in quantum billiards that can be traced back to the stable
orbits in a corresponding classical chaotic system, where such orbits leave
their imprints on the quantum wavefunction [6]. An even more important reason
for the interest in semi-classical methods is that they provide a non-
perturbative approximation scheme, which is analytically accessible -
something that is exceptional.

This strength of the approach can be seen most lucidly in that the semi-
classical method can solve problems approximately that ordinary methods such
as perturbation theory cannot. For instance, if one treats the harmonic
oscillator potential as a perturbation around the free particle, then
perturbation theory cannot produce the well-known discrete energy levels.
However, in this case, the semi-classical method is impressive because it
produces exact results for the discrete energy levels [7, 8].

This work aims to apply semi-classical methods not to the ordinary Schrödinger
equation with a quadratic dispersion, whose implications are well-known.
Instead, we want to see its consequences in the context of more exotic phases
of matter - so-called semi-metals.

After the discovery of graphene [9] (a single layer of carbon atoms arranged
in a honeycomb pattern), researchers discovered many exciting transport
properties. A famous example is Klein tunneling [10] showing that an electron
under normal incidence will tunnel through a barrier with 100% probability.
These features can be traced back to graphene’s peculiar band structure, where
the conduction and valence bands touch at isolated points. This property made
graphene an early example of a semi-metal -i.e., a material where the overlap
between conduction and valence bands is not precisely zero like in an
insulator or semiconductor and not of finite measure like in a conductor.
Instead, while the overlap is non-zero, it has zero measure. Motivated by this
discovery, much interest has emerged in similar classes of materials. One
example is Dirac semi-metals [11, 12], which were the first experimentally
confirmed material in this class. Dirac semi-metals can be considered a three-
dimensional analog of graphene [11, 13] because, much like graphene, it has a
linear dispersion near isolated band touching points. Dirac semi-metals are
interesting because the exciting linear band touching is robust. Indeed, Dirac
cones do not appear due to accidental band touchings; instead, they are
protected by crystalline symmetry and spin-orbit coupling [14, 15] and are
therefore very stable. Of course, there is a plethora of other kinds of semi-
metals with robust band features, such as nodal line semi-metals where bands
touch along a one-dimensional line [16, 17] or also semi-metals with band
touchings of higher than linear order - such as quadratic or cubic dispersion
[18, 12, 19, 20, 21].

In this work, we will focus our investigation on the exotic case of an
anisotropic cubic Dirac metal [12, 19] subjected to a constant magnetic field.
We will do so by using both fully quantum and semi-classical methods. Semi-
classical results are compared to fully quantum results. Both approaches offer
exciting insights and permit us to speculate about how one might
experimentally identify a cubic Dirac semi-metal, whether it is isotropic or
anisotropic.

We organized our work as follows. In section II, we introduce the Hamiltonian
of the system. In section III, we study the Hamiltonian for two cases that
allow for an exact quantum mechanical analytical solution - these cases are
isotropic in the $x$-$y$ plane. First, we consider the purely two-dimensional
limit. Here, the material is thin, and we may restrict ourselves to momentum
$k_{z}=0$ in the $z$-direction for low energies. Next, we consider the case
$k_{z}\neq 0$, where the material is thick enough so that not all momenta in
the $z-$ direction are frozen out. In section IV, we consider perturbatively
the effect of the anisotropy by introducing a small anisotropy parameter
$\lambda\ll 1$. In section V, we recast the problem in a semi-classical
language, which allows us to investigate classical orbits and study the
differences between exact quantum and semi-classical results. We study the
maximally anisotropic limit, which is inaccessible to a fully quantum
analytical treatment. We find classical orbits and solutions for energy levels
that are exact in the semi-classical limit and compare them to numerical
results. Lastly, in Sec. VI we draw our conclusions

##  II Anisotropic cubic Dirac semi-metals

We take the Hamiltonian describing an anisotropic cubic Dirac semi-metal [12]

$H(k)=\hbar\left(\begin{array}[]{cc}h(k)&0\\\ 0&-h(k)\end{array}\right),$ (1)

as the starting point of our discussion. Such a Hamiltonian, for instance, is
expected to be realizable in quasi-one-dimensional molybdenum monochalcogenide
compounds [12]. For such a material the basis set that was used to write the
above Hamiltonian is
$|{\Psi}\rangle=(|A,\uparrow\rangle,|B,\uparrow\rangle,|A,\downarrow\rangle,|B,\downarrow\rangle)$,
with $A/B$ and $\uparrow/\downarrow$ designating the orbital and spin degrees
of freedom, respectively. Of course, there is the possibility that similar
Hamiltonians may be realizable in other types of materials and also with
different degrees of freedom (such as sublattice degrees of freedom), giving
rise to an identical matrix structure. For this reason, our results remain
general in what follows, so we will not explicitly use the basis set
introduced above.

The Hamiltonian (1) for our case of a cubic Dirac semimetal consists of two
blocks with the form

$h(k)=v_{x}(\hat{k}^{3}_{+}+\hat{k}^{3}_{-})\sigma_{x}+iv_{y}(\hat{k}^{3}_{+}-\hat{k}^{3}_{-})\sigma_{y}+v_{z}k_{z}\sigma_{z},$
(2)

where $\hat{k}_{\pm}=\hat{k}_{x}\pm i\hat{k}_{y}$ are momentum operators with
$\hat{k}_{i}=-i\hbar\partial_{i}$. Finally, $\sigma_{i}$ are Pauli matrices,
and $v_{x,y,z}$ are real-valued and independent coefficients (they are not
velocities) with dimension $T^{-1}$, $T$ stands for time.

A constant magnetic field $B$ can be introduced into our description with a
convenient choice of the Landau gauge. That is, it is introduced via a vector
potential $\mathbf{A}=(0,Bx,0)$ using the minimal substitution procedure
${\mathbf{k}}\to{\mathbf{k}}-\frac{e}{\hbar c}\mathbf{A}$. In the unit system
($e=c=\hbar=1$), we therefore replace canonical momentum $\hat{p}_{i}$
operators by kinetic momentum operators $\hat{\pi}_{i}$ as indicated below

$\hat{k}_{\pm}\to\hat{\pi}_{\pm}=\hat{k}_{x}\pm i(\hat{k}_{y}-B\hat{x}).$ (3)

In the following section, we will first study this Hamiltonian’s energy levels
- so-called Landau levels - fully quantum mechanically, which is only possible
in specific high symmetry situations.

##  III Exact energy levels for the isotropic case $v_{x}=v_{y}$

###  III.1 Isotropic case with $k_{z}=0$

The Hamiltonian describing an isotropic cubic Dirac semi-metal with
$v_{x}=v_{y}=v$ and the $z$-direction frozen out as $k_{z}=0$ (see the
appendix of [19] for a more detailed discussion on how this happens for a thin
material), is given by

$H(\hat{\mathbf{k}},x)=2v\left(\begin{array}[]{cccc}0&\hat{\pi}^{3}_{+}&0&0\\\
\hat{\pi}^{3}_{-}&0&0&0\\\ 0&0&0&-\hat{\pi}^{3}_{+}\\\
0&0&-\hat{\pi}^{3}_{-}&0\end{array}\right).$ (4)

To solve the eigenvalue problem, we can separate variables and write the
eigenspinors as plane waves in the $y$-direction. This is due to the fact that
$[H(k),\hat{k}_{y}]=0$ implies conservation of momentum along the
$y$-direction $\Psi(x,y)=e^{ik_{y}y}\Psi(x)$ with
$\Psi(x)=\left[\Psi_{1}(x),\Psi_{2}(x),\Psi_{3}(x),\Psi_{4}(x)\right]^{T}$. We
may now define operators

$\displaystyle\hat{a}=\dfrac{1}{\sqrt{2B}}\left[\hat{k}_{x}+i(k_{y}-B\hat{x})\right],$
(5)
$\displaystyle\hat{a}^{\dagger}=\dfrac{1}{\sqrt{2B}}\left[\hat{k}_{x}-i(k_{y}-B\hat{x})\right],$
(6)

that satisfy the commutation relation $[\hat{a},\hat{a}^{\dagger}]=\mathbb{I}$
and can therefore be interpreted as ladder operators. Using this notation, the
Hamiltonian becomes

$H=\omega_{c}\left(\begin{array}[]{cccc}0&\hat{a}^{3}&0&0\\\ \hat{a}^{\dagger
3}&0&0&0\\\ 0&0&0&-\hat{a}^{3}\\\ 0&0&-\hat{a}^{\dagger
3}&0\end{array}\right),$ (7)

where $\omega_{c}=2v\left(2eB\right)^{3/2}$ is called cyclotron frequency
(recall $v$ is not a velocity). The eigenvalue problem
$H(k)\Psi(x,y)=E\Psi(x,y)$ for the second spinor $\Psi_{2}(x)$ component can
be reduced to

$\omega^{2}_{c}\hat{a}^{\dagger 3}\hat{a}^{3}\Psi_{2}(x)=E^{2}\Psi_{2}(x),$
(8)

which can be solved by harmonic oscillator states
$\Psi_{2}(x)=\big{|}n\big{>}$. Associated energies are given as

$E_{n,s}=s\omega_{c}\sqrt{n(n-1)(n-2)},$ (9)

where $s=\pm 1$. Using as an ansatz $\Psi_{2}(x)=\big{|}n\big{>}$ in this
equation, we can obtain full normalized spinors as

$\displaystyle\Psi_{n}^{(h)}(x,y)=\dfrac{1}{\sqrt{2}}\left(\begin{array}[]{c}s\big{|}n-3\big{>}\\\
\big{|}n\big{>}\\\ 0\\\ 0\end{array}\right)e^{ik_{y}y},\qquad n\geq 3,$ (14)

and using a similar relation for the fourth component, we find

$\displaystyle\Psi_{n}^{(-h)}(x,y)=\dfrac{1}{\sqrt{2}}\left(\begin{array}[]{c}0\\\
0\\\ -s\big{|}n-3\big{>}\\\
\big{|}n\big{>}\end{array}\right)e^{ik_{y}y},\qquad n\geq 3,$ (19)

Here, the label ${(h)}$ stands for spinors with contributions that come solely
from the upper block and ${(-h)}$ for the lower block in the Hamiltonian.

A special case must be considered for describing the $n<3$ states since they
require different normalization and do not come as pairs with opposing signs.
For the case $m=0,1,2$, we find the spinors

$\displaystyle\Psi_{l}^{(h)}(x,y)=\left(\begin{array}[]{c}0\\\
\big{|}l\big{>}\\\ 0\\\ 0\end{array}\right)e^{ik_{y}y},\quad l<3,$ (24)
$\displaystyle\Psi_{l}^{(-h)}(x,y)=\left(\begin{array}[]{c}0\\\ 0\\\ 0\\\
\big{|}l\big{>}\end{array}\right)e^{ik_{y}y},\quad l<3,$ (29)

where we find that the zero energy ground state is sixfold degenerate. This
observation is interesting because it contrasts that of single-layer graphene,
which has a zero energy state with no degeneracy. This ground-state degeneracy
has interesting consequences and impacts the Hall conductivity. We begin our
discussion of this with the typical formula for the Hall conductivity at
finite temperature

$\sigma_{xy}=i\sum_{n,n^{\prime}}\left(f\left(E^{n}\right)-f\left(E^{n^{\prime}}\right)\right)\frac{\left\langle\psi_{n}\left|j_{x}\right|\psi_{n^{\prime}}\right\rangle\left\langle\psi_{n^{\prime}}\left|j_{y}\right|\psi_{n}\right\rangle}{\left(E^{n}-E^{n^{\prime}}\right)^{2}},$
(30)

where $|\psi_{n}\rangle$ is any eigenstate of the Hamiltonian,
$f(E)=(\exp(\beta(E-\mu))+1)^{-1}$ is the Fermi-Dirac distribution, $\mu$ the
chemical potential, and $\beta$ the inverse temperture. Current operators are
found in the usual way by taking the derivative of the Hamiltonian with
respect to a vector potential. For our eigenstates and eigenvalues, we find a
unitless reduced Hall conductivity $\sigma_{xy}^{(r)}$, which is given as

$\displaystyle\sigma_{xy}^{(r)}$
$\displaystyle=\frac{h\sigma_{xy}}{4e^{2}}=\frac{9}{2}\omega_{c}^{2}\sum_{n,\sigma,\sigma^{\prime}}n(n-1)\frac{f({E_{n+1,\sigma}})-f({E_{n,\sigma^{\prime}}})}{(E_{n,\sigma^{\prime}}-E_{n+1,\sigma})^{2}}.$
(31)

Here, the sum is over $n\in\mathbb{N}$ and
$\sigma,\sigma^{\prime}\in\\{-1,1\\}$. The result of the computation can be
seen in Fig. 1 below.

Figure 1: Plot of the dimensionless Hall conductivity $\sigma_{xy}^{(r)}$ as
function of dimensionless chemical potential $\mu^{(r)}=\mu/\omega_{c}$.
Subfigure a) shows the Hall conductivity for our problem of a cubic Dirac
semi-metal ($\beta=50/\omega_{c}$), and Subfigure b) the case of graphene
($\beta=100/\omega_{c}$).

It shows that the step in Hall conductivity at zero chemical potential is
$3/2$ times the step size of other steps in the quantized Hall conductivity.
This result differs from the case of graphene, which we show for comparison
where the first step is half the usual step size - for a typical Dirac semi-
metal, the result has the same properties as graphene. The reason for this
behavior can be found in the six-fold degenerate lowest Landau level of the
cubic Dirac semi-metal, which has 3/2 the degeneracy of other states. This
property contrasts graphene, where the lowest Landau level has half the
degeneracy of other levels. This result also differs from a typical metal
where all steps are equally spaced. The second observation one may make is
that in graphene, steps occur at a dimensionless chemical potential
$\mu^{(r)}=\sqrt{n}$. In contrast, in the cubic Dirac semi-metal case, they
occur at values that asymptotically go to $\mu^{(r)}=n^{3/2}$. Both
observations combined with the step height (it is twice that of a cubic Weyl
semi-metal [22]) - which is due to the level degeneracy of the cubic Dirac
semi-metal’s Landau levels - can be used to identify such materials via the
quantized Hall conductivity experimentally.

###  III.2 Quasi-2D isotropic case $k_{z}\neq 0$

After considering the case of semi-metal with $k_{z}=0$ (in addition to being
isotropic in the $x$-$y$ plane), we now consider the case where the dimension
of the material in $z$-direction is sufficiently large such that considering
non-zero momentum $k_{z}\neq 0$ is not only relevant for high energies. As
previously, we will restrict our discussion to the isotropic case. To keep our
description consistent in units of the cyclotron frequency $\omega_{c}$, we
define $m=\frac{v_{z}k_{z}}{2v\left(\sqrt{2B}\right)^{3}}$, which can be
interpreted as a mass-like term in the original Hamiltonian model (it causes a
gap in the excitation spectrum). The Hamiltonian then takes the form

$H=\omega_{c}\left(\begin{array}[]{cccc}m&\hat{a}^{3}&0&0\\\ \hat{a}^{\dagger
3}&-m&0&0\\\ 0&0&-m&-\hat{a}^{3}\\\ 0&0&-\hat{a}^{\dagger
3}&m\end{array}\right).$ (32)

From here, finding the Landau levels is almost immediate and follows almost
the same steps as in the previous section. We find

$E_{n,s}(m)=s\omega_{c}\sqrt{n(n-1)(n-2)+m^{2}},\quad n\geq 3,$ (33)

with $s=\pm 1$, where normalized eigenspinors are given as

$\displaystyle\Psi_{n,s}^{(h)}(x,y)=N\left(\begin{array}[]{c}sc_{n}\big{|}n-3\big{>}\\\
\big{|}n\big{>}\\\ 0\\\ 0\end{array}\right)e^{ik_{y}y},$ (38)
$\displaystyle\Psi_{n}^{(-h)}(x,y)=N\left(\begin{array}[]{c}0\\\ 0\\\
-sc_{n}\big{|}n-3\big{>}\\\ \big{|}n\big{>}\end{array}\right)e^{ik_{y}y},$
(43)

with the constants

$\displaystyle c_{n}=\dfrac{E_{n,s}(0)}{E_{n,s}(m)-m},\ \
N=\dfrac{1}{\sqrt{c_{n}^{2}+1}}.$ (44)

It is important to note that the case of $n<3$ again requires special care and
has the same eigenspinors as given in Eq. (29) with corresponding eigenvalues

$E_{n}^{(\pm h)}=\mp\omega_{c}m,$ (45)

which are each thrice degenerate. It is important to stress that the ”$-$”
sign corresponds to the upper block of the Hamiltonian and the ”$+$” sign to
the lower block, which is unlike the $n>3$ result that has both signs in both
cases. Much like the truly one-dimensional case, we can plot a Hall
conductivity and find the result in Fig. 2 below

Figure 2: Plot of the dimensionless Hall conductivity $\sigma_{xy}^{(r)}$ as
function of dimensionless chemical potential $\mu^{(r)}=\mu/\omega_{c}$ at
inverse temperature $\beta=50/\omega_{c}$ and with ”mass” term $m=-3$.

We observe that the effect of the mass term $m$ is that some plateaus have
been broadened in width while others have been shrunk. This behavior would
allow experimentalists to recognize if electrons carry momentum in the
z-direction.

##  IV Perturbative energy levels for the case $v_{y}\neq v_{x}$

While it was easy to find exact energy levels in the isotropic case of
$v_{y}=v_{x}$, this task is not simple using analytical means for the
anisotropic case of $v_{y}\neq\ v_{x}$. Introducing an anisotropy parameter
$\lambda=\frac{1-v_{y}/v_{x}}{2}$ is useful. This approach allows us to write
the Hamiltonian in a convenient form

$H=\omega_{c}(h_{0}+\lambda h_{1})\otimes\sigma_{z},$ (46)

with $\lambda$ as perturbative parameter and definitions

$\displaystyle h_{0}=\left(\begin{array}[]{cc}0&\hat{a}^{3}\\\
\hat{a}^{\dagger 3}&0\end{array}\right),\
h_{1}=\left(\begin{array}[]{cc}0&\hat{a}^{\dagger 3}-\hat{a}^{3}\\\
\hat{a}^{3}-\hat{a}^{\dagger 3}&0\end{array}\right).$ (51)

With this, $h_{0}$ corresponds to an exactly solvable isotropic Hamiltonian,
and $\lambda h_{1}$ is considered a perturbation. Eigenvalues $E^{0}_{n}$ are
non-degenerate in the case $n\geq 3$ and corrections to the energy to first
order can be computed straightforwardly by the standard expressions for
Reighleigh-Schrödinger perturbation theory [23] as

$\displaystyle E_{n}=E_{n}^{0}\left(1-\lambda\right).$ (52)

Like in previous sections, the degenerate zero eigenvalues case requires
special care and must be treated using degenerate perturbation theory.
Interestingly, first-order degenerate perturbation theory yields no correction
such that the expression (52) that includes the first-order corrections to
energies above remains valid in all cases. Therefore, it is clear that a small
anisotropy is not expected to have much of an impact on the qualitative
picture of the Hall conductivity (to leading order
$\omega\to(1-\lambda)\omega_{c}$ is the only modification).

##  V Semi-classical treatment

Since a completely quantum mechanical treatment for the anisotropic case is
challenging, we will resort to a semi-classical treatment in this section.
This approach allows us to gain further insights, such as the shapes of
particle trajectories under the influence of a magnetic field, and further
insights into the anisotropic regime, which we stress is difficult to access
in a complete quantum treatment.

###  V.1 Review of Bohr-Sommerfeld-type semiclassics for matrix Hamiltonians

A general and powerful method to determine energy levels by studying classical
trajectories is the so-called Gutzwiller approach to the semi-classical
density of states [1]. This approach has been generalized to arbitrary matrix-
valued Hamiltonians in [3]. While the method developed in [3] is quite general
and allows for treatments of various special cases, our current problem of
determining Landau levels for a cubic Dirac semi-metal is much simpler.
Notably, we have three significant simplifications:

  1. 1.

It is enough to restrict ourselves to one-dimensional Hamiltonians. We may do
so because, in our case of Landau levels, $k_{y}$ and $k_{z}$ can be treated
as parameters entering a 1D Hamiltonian - that is, as complex numbers rather
than bonafide operators.

  2. 2.

The eigenvalues of the classical limit Hamiltonian $H(p,x)$, where the
momentum operator $\hat{p}$ was mapped to a complex number $p$, are non-
degenerate. This simplification is possible because the $h(k)$ blocks in
equation (1) are not coupled, and it is, therefore, enough to consider them
separately, and their eigenvalues can be found to be non-degenerate.

  3. 3.

We may ignore degeneracy factors that arise from different Landau orbits over
the plane. This simplification is possible because they only enter the density
of states and not in expressions for the energy levels that interest us.

These three simplifications mean that in what follows, we may directly employ
the modified Bohr-Sommerfeld quantization condition that was derived in [3]
and is given below

$2S_{\alpha}-\hbar\left[\pi(\nu_{\alpha}+4n)-2\int_{0}^{T_{\alpha}}dtM_{\alpha}\right]=0.$
(53)

In what follows, we will choose units with $\hbar=1$. We now move on to
explain the different terms that enter this expression. Hereby, $S_{\alpha}$
is the classical action for a primitive periodic orbit - an orbit in phase
space that has only been traversed once. The index $\alpha$ accounts for
different periodic orbits being possible at a fixed energy. There are two ways
this can happen.

  1. 1.

Potentials can allow for different periodic orbits at a given energy $E$ like
in Fig. 3 below.

Figure 3: Plot of a potential landscape $V(x)$ (thick black line) and energy
as a dashed line. Two distinct potential pots are visible. Each will have its
own associated classical periodic orbits with an action $S_{\alpha}$.

  2. 2.

We deal with matrix Hamiltonians. That is, to the lowest order, the
Schrödinger equation is replaced by a simple matrix equation

$H(x,p_{x}=\partial_{x}S)V=EV.$ (54)

For an $n\times n$ Hamiltonian this can be diagonalized to determine $n$
eigenvectors $V_{i}$ and $n$ eigenvalue equations that determine actions
$S_{\alpha}$

$H_{i}(x,p_{x}=\partial_{x}S_{\alpha})=E,$ (55)

which are ordinary Hamilton-Jacobi equations and $H_{i}$ are eigenvalues of
$H(x,\partial_{x}S)$. These Hamilton-Jacobi equations can have one or more
periodic orbits associated with a given energy $E$.

Next, in our expression, Eq. 53 enters $\nu_{\alpha}$, which is the so-called
Maslov index which for 1D motion counts the number of reflections by a
potential barrier (or more generally directional changes) that a particle
experiences while it traverses a periodic orbit. The last term that enters is
the so-called semi-classical phase factor

$M_{\alpha}=\mathrm{Im}(V_{\alpha}^{\dagger}[\partial_{p_{x}}H(x,p_{x})]\partial_{x}V_{\alpha}),$
(56)

which bears some resemblance to a Berry phase as discussed in [24, 25, 26, 27,
3]. Hereby, $V_{\alpha}$ is the eigenvector of the classical Hamiltonian
matrix $H(x,\partial_{x}S_{\alpha})$ corresponding to the orbit $\alpha$.

###  V.2 The Hamiltonian system

For our computations of semi-classical energies, we will consider the upper
block of (1) only, that is, Eq. (2) with the magnetic field introduced as
before via minimal substitution. This simplification is possible because the
second block in equation (1) will yield (up to a sign) the same energies as
the other block. The semi-classical limit of the Hamiltonian (2) is then
obtained by replacing momentum operators with canonical momenta
$\hat{p}_{i}\rightarrow\partial_{x_{i}}S=p_{i}$. We obtain the classical
expression

$h(\pi_{i})=2v_{x}\begin{pmatrix}0&\pi_{+}^{3}+\lambda(\pi_{-}^{3}-\pi_{+}^{3})\\\
\pi_{-}^{3}+\lambda(\pi_{+}^{3}-\pi_{-}^{3})&0\end{pmatrix},$ (57)

and we introduced the shorthand notations

$\displaystyle\pi_{\pm}=p_{x}\pm i(p_{y}+eBx),$ (58)

where $p_{i}$ are canonical momenta and $\pi_{i}$ are kinetic momenta. The
parameter $\lambda=\frac{1-v_{y}/v_{x}}{2}$ measures the degree of anisotropy
as discussed previously.

###  V.3 Isotropic case $k_{z}=0$

In the isotropic case, we may set $\lambda=0$ and $v_{x}=v_{y}=v$ and for a
thin sheet, low energies correspond to $k_{z}=0$ to find

$h(\pi)=2v\left(\begin{array}[]{cc}0&\pi^{3}_{+}\\\ \pi^{3}_{-}&0\\\
\end{array}\right).$ (59)

The eigenvalues are obtained from the determinant as $\det[h(\pi_{i})-E]=0$

$\begin{split}E_{s}=&2sv\sqrt{(\pi_{+}\pi_{-})^{3}}\end{split},$ (60)

where $s=\pm 1$. The corresponding normalized eigenvectors are given by

$V=\dfrac{1}{\sqrt{2}}\begin{pmatrix}\frac{2v}{E_{s}}\pi_{+}^{3}\\\
1\end{pmatrix}.$ (61)

At this level, we see that there is two classical Hamiltonians of the form

$H(p_{i},x)=2sv\left(p_{x}^{2}+\left(p_{y}+eBx\right)^{2}\right)^{3/2}.$ (62)

####  V.3.1 Time for an orbit

Another ingredient that enters our expression for the generalized Bohr-
Sommerfeld quantization condition is the time for an orbit. Here, we start by
investigating the equations of motion, which are given by Hamilton’s equations

$\dot{q}=\dfrac{\partial H}{\partial p},\quad\dot{p}=-\dfrac{\partial
H}{\partial q},$ (63)

with $q$ is the position and $p$ is momentum coordinates. In our case, we find

$\displaystyle\dot{x}=3s\left(4|E|v^{2}\right)^{\frac{1}{3}}p_{x},$ (64)
$\displaystyle\dot{y}=3s\left(4|E|v^{2}\right)^{\frac{1}{3}}(p_{y}+eBx),$
$\displaystyle\dot{p}_{x}=3seB\left(4|E|v^{2}\right)^{\frac{1}{3}}(p_{y}+eBx),$
$\displaystyle\dot{p}_{y}=0,$

From (64) one finds (taking a time derivative of $\dot{x}$ and using the
equation for $\dot{p}_{x}$) that the equation of motion for the $x$-coordinate
is given as

$\ddot{x}=-\omega^{2}(x+x_{s}),$ (65)

which constitutes a shifted harmonic oscillator with position shift
$x_{s}=p_{y}/(eB)$ (note that $p_{y}$ is constant as shown above). It is
important to note that this equation is valid for both classical Hamiltonians
in Eq. (62). It also means that it is enough to consider only one of the
eigenvalues to find semi-classical energy levels. Indeed, in the remainder of
this section, we will focus on only one of the Hamiltonians and, therefore,
drop all indices related to orbits. The oscillator frequency is found to be

$\omega=3(4|E|v^{2})^{\frac{1}{3}}eB.$ (66)

The time for an orbit as per definition $\omega=2\pi/T$ then is then found to
be

$T=\dfrac{2\pi}{3(4|E|v^{2})^{\frac{1}{3}}eB}.$ (67)

####  V.3.2 Classical orbit, action, and Maslov index

Now, from Eq. (65), we find the solutions

$\displaystyle x=x_{s}+x_{0}\sin(\omega t),$ $\displaystyle
p_{x}=seBx_{0}\cos\left(\omega t\right),$ (68)

with amplitude

$x_{0}=\frac{1}{eB}\left(\frac{|E|}{2v}\right)^{\frac{1}{3}},$ (69)

which can be found for a given energy by placing the solutions (68) in Eq.
(60). Plotted, we get the typical phase space curve of a harmonic oscillator
that can be seen below

Figure 4: Phase space curve of an electron in cubic semi-metal subjected to a
constant magnetic field.

which is just an elliptic orbit. This curve directly tells us that the
particle along each orbit changes direction twice, and therefore, the Maslov
index is

$\nu_{\alpha}=2.$ (70)

We may next use $p_{x}=\partial_{x}S$ in Eq. (60) to find the action for an
orbit as

$\begin{split}S&=\int\sqrt{\left(\frac{|E|}{2v}\right)^{\frac{2}{3}}-(p_{y}+eBx)^{2}}dx\\\
\end{split}.$ (71)

This integral for a single period of a primitive orbit gives

$S=\frac{\pi}{eB}\left(\frac{|E|}{2v}\right)^{\frac{2}{3}}.$ (72)

####  V.3.3 Semi-classical Berry-like phase and energy levels

For the semi-classical phase term using (56) we find

$\displaystyle M_{\pm}=9eB\left(\frac{|E|v^{2}}{2}\right)^{\frac{1}{3}},$ (73)

which is a constant, and therefore the integral $\int_{0}^{T}dtM=MT=3\pi$ is
trivial.

We may use our results (72), (70), (73) and put them into (53) to find that
our semiclassical energies are then given by

$E=\pm\omega_{c}(n-1)^{3/2},$ (74)

with $\omega_{c}=2v\left(2eB\right)^{3/2}$.

####  V.3.4 Performance of semi-classical results

Understanding under what conditions the semi-classical results perform well
compared to the exact results obtained earlier is helpful. Of course, we see
directly that the results are different. This observation poses the question
under which conditions the results agree.

Typically, semi-classical results are reliable for large quantum numbers.
Therefore, we do an expansion for large $n$ of both the exact
$E_{\mathrm{exact}}$ and semi-classical $E_{\mathrm{semcl}}$ results to find

$\displaystyle
E_{\mathrm{exact}}=\pm\omega_{c}\left(n^{3/2}-\frac{3}{2}n^{\frac{1}{2}}+\frac{3}{8}n^{-\frac{1}{2}}\right)+\mathcal{O}\left(n^{-\frac{3}{3}}\right),$
(75) $\displaystyle
E_{\mathrm{semcl}}=\pm\omega_{c}\left(n^{3/2}-\frac{3}{2}n^{\frac{1}{2}}-\frac{1}{8}n^{-\frac{1}{2}}\right)+\mathcal{O}\left(n^{-\frac{3}{3}}\right),$

where we can see that results agree well in the limit of large quantum numbers
- they only disagree at order $\mathcal{O}(n^{-1/2})$ as it is typical in the
semi-classical approach - because typically larger quantum numbers correspond
to a larger classical action. To better understand how accurate the
approximation is, we plot the density of states below and energy levels as a
function of $n$. The approximation’s effectiveness is also visualized in the
figure below.

Figure 5: (color online) The left column of the figure shows a plot of the
density of states, and the right column shows the relative error between exact
and approximate energies.

The semi-classical energy levels agree very well with exact results except for
the lowest energy. Furthermore, as expected, the approximation error decreases
rapidly as a function of the quantum number $n$.

###  V.4 Quasi-2D isotropic limit

Next, we consider the case $p_{z}\neq 0$, which can describe a sheet of finite
thickness. To obtain a simpler notation that we can compare to the exact
quantum case, we set $p_{z}=m\omega_{c}/(2v_{z})$. The notation is suggestive
and clarifies that $p_{z}$ takes the role of a mass term $m$. The Hamiltonian
then takes the form

$h(\pi)=2v\left(\begin{array}[]{cc}\frac{m\omega_{c}}{2v}&\pi^{3}_{+}\\\
\pi^{3}_{-}&-\frac{m\omega_{c}}{2v}\\\ \end{array}\right).$ (76)

The eigenvalues are readily found as

$E=2vs\sqrt{(\pi_{+}\pi_{-})^{3}+\left(\frac{m\omega_{c}}{2v}\right)^{2}},$
(77)

with $s$ the sign of the energy. The corresponding eigenvectors can be written
as

$\displaystyle
V=\sqrt{2}\sqrt{1-\frac{m\omega_{c}}{E}}\begin{pmatrix}\frac{\pi_{+}^{3}v}{E-2m\omega_{c}}\\\
1\end{pmatrix}.$ (78)

The classical system will reduce to a harmonic oscillator after using the same
methods as in the previous section. However, this time with a slightly
modified driving frequency

$\omega=3\left(\frac{2v\left(\left|E|^{2}-m^{2}\omega_{c}^{2}\right.\right)}{|E|^{3/2}}\right)^{\frac{2}{3}}eB.$
(79)

Considering the limit of small mass and the relation to the frequency in the
massless case, Eq. (66) is interesting. We will use the term $\omega_{0}$ to
refer to the frequency of the massless case below. In the limit of small
masses, the oscillation frequency becomes

$\omega\approx\omega_{0}-\frac{2}{3}\frac{m^{2}\omega_{c}^{2}}{|E|^{2}}\omega_{0},$
(80)

which tells us that a mass term reduces the oscillation frequency just as
intuition would dictate. Because we have harmonic motion, we find a Maslov
index

$\nu=2.$ (81)

The action of a primitive can be computed in analogy to the previous case and
is given as

$S=\frac{\pi(|E|^{2}-m^{2}\omega_{c}^{2})^{\frac{1}{3}}}{2^{2/3}eBv^{2/3}}.$
(82)

Lastly, we find that the semi-classical phase factor has the form

$M=\frac{9Be\left(v\left(|E|^{2}-4m^{2}\omega^{2}/v^{2}\right)\right)^{2/3}}{2^{\frac{1}{3}}E},$
(83)

which is a constant. Therefore, the integral $\int_{0}^{T}dtM=2\pi
M/\omega=3\pi$ is trivial.

Applying the Bohr-Sommerfeld quantization condition Eq. (53) as before, we
find that energy levels are given as

$E_{\text{semcl}}=\omega_{c}\sqrt{m^{2}+(n-1)^{3}}.$ (84)

We may now compare the semi-classical result $E_{\text{semcl}}$ to the exact
result Eq. (33), which we title $E_{\text{exact}}$ to find

$E_{\text{semcl}}-E_{\text{exact}}=(n-1)\frac{\omega_{c}}{2m}+\mathcal{O}(m^{-3}),$
(85)

which tells us that, similar to expectations, a mass term turns the energy
levels more classical (the approximation error shrinks for large mass).

###  V.5 Maximally anisotropic case $v_{x}=0$

In the last part of the section, we consider a case that is difficult to solve
by analytical means in an exact quantum treatment - the case of a maximally
anisotropic cubic Dirac semi-metal. Maximum anisotropy is achieved for
$v_{x}=0$ and $v_{y}=v$, where the classical Hamiltonian takes the form

$H=v\begin{pmatrix}0&\pi_{-}^{3}-\pi_{+}^{3}\\\
\pi_{+}^{3}-\pi_{-}^{3}&0\end{pmatrix}.$ (86)

For this case, it is advantageous to combine Hamilton-Jacobi equations for
different eigenvalues into a single Hamilton-Jacobi equation $\det(E-H)=0$,
and one finds

$E^{2}+v^{2}\left(\pi_{-}^{3}-\pi_{+}^{3}\right)^{2}=0.$ (87)

Eigenvectors that correspond to different solutions for $E$ are given as

$V_{\pm}=\frac{1}{\sqrt{2}}\begin{pmatrix}\pm 1\\\ 1\end{pmatrix},$ (88)

and are independent of momenta and positions, which directly lets one realize
that the semi-classical phase $M$ will have no contributions. It is,
therefore, important to recognize that the main contribution to orbits will be
from the classical phase space orbit. The combined Hamilton-Jacobi equation
(87) can be interpreted as an implicit relation for the orbit in phase space

$E^{2}-4v^{2}(Bex+p_{y})^{2}\left((Bex+p_{y})^{2}-3p_{x}^{2}\right)^{2}=0,$
(89)

which leads to an interesting phase space orbit shown in Fig. 6.

Figure 6: Plot of phase space trajectories in terms of unitless momentum
$\tilde{p}_{i}=p_{i}(eB)^{-1/2}$, unitless energy
$\epsilon=Ev^{-1}(eB)^{-3/2}$ and position $\tilde{x}=x(eB)^{1/2}$ (recall $v$
is not a velocity and we set $\hbar=1$). In both cases, we set $p_{y}=0$
because this term only leads to a shift of the trajectory center along the
$x$-axis. The red curve is for a value $\epsilon=1$ and the blue one for a
value $\epsilon=0.02$

We observe that semi-classical particles in this problem traverse star-shaped
orbits in phase space. These orbits do not close except at infinity (for a
non-zero value of $v_{x}$, the closure happens at a finite value, though). It
is clear from the figure that despite orbits that only close at infinity, the
phase space volume enclosed by the trajectory is finite. This observation is
interesting and peculiar because the finite phase space volume implies
discretized energy levels corresponding to localized states. At the same time,
the classical particles are not localized - they move to infinity along the
jags of the star. This observation demonstrates how a particle can be
localized for a quantum mechanical problem. The classical dynamics do not
reflect this in the real space trajectory but only in phase space through a
finite volume. The Maslov index can also easily be found from the figure if we
note that barrier reflections happen when the momentum derivative becomes
infinite - this happens at four spots along the orbit (the two spikes at
$\tilde{x}=0$ and the indents at $\tilde{p}_{x}=0$). Therefore, the Maslov
index in equation (53) is $\nu=4$.

The last ingredient in equation (53) to find semi-classical energies is the
action, which is the volume enclosed by the orbit. The simplest way to do this
is to realize that it is enough to compute the area of one of the four star
sectors, which is marked in blue Fig. 7 below.

Figure 7: Plot of one sector of the star orbit. Marked in blue is the area one
needs to compute. The blue curve is given by
$\tilde{p}_{x}=\sqrt{x^{2}/3+\epsilon/(6x)}$ and the orange curve by
$\tilde{p}_{x}=\sqrt{x^{2}/3-\epsilon/(6x)}$

We find that the area can be separated into different pieces, giving a
unitless action $s$ with

$s/4=\int_{0}^{\infty}dx\sqrt{\frac{2x^{3}+\epsilon}{6x}}-\int_{\epsilon^{1/3}2^{-1/3}}^{\infty}dx\sqrt{\frac{2x^{3}-\epsilon}{6x}}.$
(90)

The integral simplifies considerably to give

$s=\frac{2^{4/3}\sqrt{\pi}\epsilon^{2/3}\Gamma\left(\frac{7}{6}\right)}{\Gamma\left(\frac{5}{3}\right)},$
(91)

where $\Gamma$ is the Gamma function. After units are reintroduced, the result
is inserted into the Bohr-Sommerfeld quantization condition, and the resulting
relation solved for the energy. We find the semi-classical energy levels to be
given by

$E_{n}=\hbar\omega_{c}\frac{\pi^{3/4}\Gamma\left(\frac{5}{3}\right)^{3/2}}{8\Gamma\left(\frac{7}{6}\right)^{3/2}}n^{3/2},$
(92)

which much like the isotropic case behave as $E_{n}\propto n^{3/2}$ just with
a different proportionality constant. Of course, it is interesting that we
obtained relatively simple closed-form results semiclassically, where a full
quantum treatment is difficult. Even a proper numerical treatment is
numerically expensive because particles are only weakly localized (for
details, see appendix A). We compare our closed-form semi-classical
approximation to the exact numeric results below to find out if our
approximation is reliable.

Figure 8: (color online) The left column of the figure shows a plot of the
density of states, and the right column shows the relative error between exact
and approximate energies.

The results are in excellent agreement with a complete numerical treatment.
Therefore, we can see that the classically non-localized particles are weakly
quantum localized. Moreover, the functional form of the energies tells us that
even in an anisotropic case, one can expect similar behavior for the Hall
conductivity as we have predicted for the isotropic case.

##  VI Conclusion

In summary, we have shown that the Landau levels of a cubic Dirac semi-metal
exhibit intriguing features both quantum mechanically and from a classical
perspective. We have shown in a fully quantum mechanical treatment that one
would expect experimentally detectable signatures of a cubic Dirac semi-metal
in the Hall conductivity, and we have found explicit analytical expressions.
While the anisotropic case could not be understood by employing a full quantum
treatment, we made considerably much more progress on the semiclassical side.

Interestingly, we found that in the anisotropic case, the electrons do not
localize classically. Rather they move infinitely far from the origin.
However, they were found to be quantum mechanically weakly localized, giving
rise to discrete energy levels. Our work also provides an exciting example of
quantum corrections playing a significant role in localization, enriching the
phenomenological literature in this respect. Nevertheless, weak localization
from the classical side has left its imprint in that the numerics for finding
exact quantum energy levels are relatively expensive when using a local set of
basis functions.

##  VII Acknowledgements

The authors gratefully acknowledge the support provided by the Deanship of
Research Oversight and Coordination (DROC) at King Fahd University of
Petroleum & Minerals (KFUPM) for funding this work through exploratory
research grant No. ER221002.

## References

  * Gutzwiller [1971] M. C. Gutzwiller, Periodic Orbits and Classical Quantization Conditions, Journal of Mathematical Physics 12, 343 (1971), https://pubs.aip.org/aip/jmp/article-pdf/12/3/343/19155865/343_1_online.pdf .
  * Kramers [1926] H. A. Kramers, Wellenmechanik und halbzahlige quantisierung, Zeitschrift für Physik 39, 828 (1926).
  * Vogl  _et al._ [2017] M. Vogl, O. Pankratov, and S. Shallcross, Semiclassics for matrix hamiltonians: The gutzwiller trace formula with applications to graphene-type systems, Physical Review B 96, 035442 (2017).
  * Loewe Yanez [2017] B. A. Loewe Yanez,  _Bridges between quantum and classical mechanics: Directed polymers, flocking and transitionless quanum driving_ , Ph.D. thesis, Georgia Institute of Technology (2017).
  * Cvitanović  _et al._ [2016] P. Cvitanović, R. Artuso, R. Mainieri, G. Tanner, and G. Vattay,  _Chaos: Classical and Quantum_ (Niels Bohr Inst., Copenhagen, 2016).
  * Heller [1984] E. J. Heller, Bound-state eigenfunctions of classically chaotic hamiltonian systems: Scars of periodic orbits, Phys. Rev. Lett. 53, 1515 (1984).
  * [7] Semiclassical time evolution amplitude, in  _Path Integrals in Quantum Mechanics, Statistics, Polymer Physics, and Financial Markets_ (World Scientific) pp. 369–457.
  * Amiot and Marleau [2006] P. Amiot and L. Marleau, Mécanique classique ii, Département de physique, de génie physique et d’optique. Université de Laval, Québec, Canada (2006).
  * Novoselov  _et al._ [2004] K. S. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang, Y. Zhang, S. V. Dubonos, I. V. Grigorieva, and A. A. Firsov, Electric field effect in atomically thin carbon films, Science 306, 666 (2004), https://www.science.org/doi/pdf/10.1126/science.1102896 .
  * Katsnelson  _et al._ [2006] M. I. Katsnelson, K. S. Novoselov, and A. K. Geim, Chiral tunnelling and the klein paradox in graphene, Nature Physics 2, 620 (2006).
  * Liu  _et al._ [2014a] Z. Liu, B. Zhou, Y. Zhang, Z. Wang, H. Weng, D. Prabhakaran, S.-K. Mo, Z. Shen, Z. Fang, X. Dai, _et al._ , Discovery of a three-dimensional topological dirac semimetal, $\text{Na}_{3}\text{Bi}$, Science 343, 864 (2014a).
  * Liu and Zunger [2017] Q. Liu and A. Zunger, Predicted realization of cubic dirac fermion in quasi-one-dimensional transition-metal monochalcogenides, Physical Review X 7, 021019 (2017).
  * Liu  _et al._ [2014b] Z. K. Liu, J. Jiang, B. Zhou, Z. J. Wang, Y. Zhang, H. M. Weng, D. Prabhakaran, S.-K. Mo, H. Peng, P. Dudin, T. Kim, M. Hoesch, Z. Fang, X. Dai, Z. X. Shen, D. L. Feng, Z. Hussain, and Y. L. Chen, A stable three-dimensional topological dirac semimetal $\text{Cd}_{3}\text{As}_{2}$, Nature materials 13, 677 (2014b).
  * Kushwaha  _et al._ [2015] S. K. Kushwaha, J. W. Krizan, B. E. Feldman, A. Gyenis, M. T. Randeria, J. Xiong, S.-Y. Xu, N. Alidoust, I. Belopolski, T. Liang, _et al._ , Bulk crystal growth and electronic characterization of the 3d dirac semimetal $\text{Na}_{3}\text{Bi}$, APL materials 3, 041504 (2015).
  * Han  _et al._ [2020] J. Han, S.-W. Kim, W. Kyung, C. Kim, G. Cao, X. Chen, S. Wilson, S. Cheon, and J. Lee, Nonsymmorphic dirac semimetal and carrier dynamics in the doped spin-orbit-coupled mott insulator $\text{Sr}_{2}\text{IrO}_{4}$, Physical Review B 102, 041108 (2020).
  * Fang  _et al._ [2016] C. Fang, H. Weng, X. Dai, and Z. Fang, Topological nodal line semimetals, Chinese Physics B 25, 117106 (2016).
  * Yu  _et al._ [2017] R. Yu, Z. Fang, X. Dai, and H. Weng, Topological nodal line semimetals predicted from first-principles calculations, Frontiers of Physics 12, 127202 (2017).
  * Fang  _et al._ [2012] C. Fang, M. J. Gilbert, X. Dai, and B. A. Bernevig, Multi-weyl topological semimetals stabilized by point group symmetry, Physical Review Letters 108, 10.1103/physrevlett.108.266802 (2012).
  * Bouhlal  _et al._ [2021] A. Bouhlal, A. Jellal, H. Bahlouli, and M. Vogl, Tunneling in an anisotropic cubic dirac semi-metal, Annals of Physics 432, 168563 (2021).
  * Bouhlal  _et al._ [2022] A. Bouhlal, A. Abbout, A. Jellal, H. Bahlouli, and M. Vogl, Tunneling phase diagrams in anisotropic multi‐weyl semimetals, Annalen der Physik 534, 10.1002/andp.202200267 (2022).
  * Jin  _et al._ [2022] Y. J. Jin, Y. Xu, Z. J. Chen, and H. Xu, Type-ii quadratic and cubic weyl fermions, Phys. Rev. B 105, 035141 (2022).
  * Barlas  _et al._ [2012] Y. Barlas, K. Yang, and A. H. MacDonald, Quantum hall effects in graphene-based two-dimensional electron systems, Nanotechnology 23, 052001 (2012).
  * Schrödinger [1926] E. Schrödinger, Quantisierung als eigenwertproblem, Annalen der Physik 385, 437 (1926), https://onlinelibrary.wiley.com/doi/pdf/10.1002/ andp.19263851302 .
  * Carmier and Ullmo [2008] P. Carmier and D. Ullmo, Berry phase in graphene: Semiclassical perspective, Phys. Rev. B 77, 245413 (2008).
  * Keppeler [2003] S. Keppeler, Semiclassical quantisation rules for the dirac and pauli equations, Annals of Physics 304, 40 (2003).
  * Bolte and Keppeler [1999] J. Bolte and S. Keppeler, A semiclassical approach to the dirac equation, Annals of Physics 274, 125 (1999).
  * Bolte and Keppeler [1998] J. Bolte and S. Keppeler, Semiclassical time evolution and trace formula for relativistic spin-1/2 particles, Phys. Rev. Lett. 81, 1987 (1998).

##  Appendix A Numerical solution of the maximally anisotropic case

For a numerical treatment, we have expressed the Hamiltonian (86) in terms of
the usual Landau basis as

$h=\frac{\omega_{c}}{2}\begin{pmatrix}0&{a^{{\dagger}}}^{3}-a^{3}\\\
{a}^{3}-{a^{{\dagger}}}^{3}&0\end{pmatrix},$ (93)

and we used that we may write

$(a^{3})_{ij}=\sqrt{2i+3i^{2}+i^{3}}\delta_{i+3,j}.$ (94)

A naive expansion letting $i$ and $j$ run from zero to a cut-off does not
yield valid results. Indeed, as one can already check in the isotropic case,
such an approach leads to issues with unphysical states at zero energy. One
finds a sixfold fold degeneracy instead of the actual physical threefold
degeneracy. Luckily, there is a simple way to fix this. If one replaces

$a^{3}\to(a^{3})_{\text{chop,x,y}},$ (95)

where $(a^{3})_{\text{c,x}}$ has the last three columns chopped from the
matrix and $(a^{3})_{\text{c,y}}$ has the last three rows chopped. In the
isotropic case, this replacement leads to the correct number of zero modes and
the correct eigenvectors for the zero modes, which justifies this replacement.

For the anisotropic case, we may then write the numerical Hamiltonian as

$h_{\text{chop}}=\frac{\omega_{c}}{2}\begin{pmatrix}0&({a^{{\dagger}}}^{3})_{\text{c,x}}-(a^{3})_{\text{c,x}}\\\
({a}^{3})_{\text{c,x}}-({a^{{\dagger}}}^{3})_{\text{c,y}}&0\end{pmatrix},$
(96)

which leads to the results employed in the main text. We should also note that
this approach is slowly converging, so a basis of $\sim 2^{22}$ states was
needed (close to state-of-the-art system sizes in exact diagonalization) to
generate results from the main text. The resulting Python code employing
sparse matrix diagonalization is available from the authors upon reasonable
request.