\section{Introduction}
Knowledge graphs (KGs) represent factual information in the form of triplets $(s,r,o)$, e.g., (\textit{Joe Biden}, \textit{is president of}, \textit{USA}), where $s$, $o$ are the subject and the object of a fact, and $r$ is the relation between $s$ and $o$. KGs have been extensively used to aid the downstream tasks in the field of artificial intelligence, e.g., recommender systems \cite{DBLP:conf/aaai/WangWX00C19} and question answering \cite{DBLP:conf/aaai/ZhangDKSS18}. By incorporating time information into KGs, temporal knowledge graphs (TKGs) represent every fact with a quadruple $(s,r,o,t)$, where $t$ denotes the temporal constraint specifying the time validity of the fact. With the introduction of temporal constraints, TKGs are able to describe the ever-changing knowledge of the world. For example, due to the evolution of world knowledge, the fact (\textit{Angela Merkel}, \textit{is chancellor of}, \textit{Germany}) is valid only before (\textit{Olaf Scholz}, \textit{is chancellor of}, \textit{Germany}). TKGs naturally capture the evolution of relational facts in a time-varying context.

Though KGs are constructed with large-scale data, they still suffer from the problem of incompleteness \cite{DBLP:conf/naacl/MinGW0G13}. Hence, there has been extensive work aiming to propose KG reasoning models to infer the missing facts, i.e., the missing links, in KGs \cite{DBLP:conf/nips/BordesUGWY13,DBLP:conf/icml/TrouillonWRGB16,DBLP:conf/iclr/SunDNT19}. Similar to KGs, TKGs are also known to be highly incomplete. This draws huge attention to developing methods for link inference in TKGs \cite{DBLP:conf/www/LeblayC18,Ma2019EmbeddingMF,DBLP:conf/iclr/LacroixOU20}. Most of these methods require a huge amount of data for all the relations to learn expressive representations, however, Xiong et al. \cite{DBLP:conf/emnlp/XiongYCGW18} and Mirtaheri et al. \cite{DBLP:conf/akbc/MirtaheriR0MG21} find that a large portion of KG and TKG relations are long-tail (i.e., these relations only occur for a handful of times), and this leads to the degenerated link inference performance of the existing KG and TKG reasoning methods. To tackle this problem, a line of few-shot learning (FSL) methods \cite{DBLP:conf/emnlp/XiongYCGW18,DBLP:conf/emnlp/ChenZZCC19,DBLP:conf/aaai/ZhangYHJLC20,DBLP:conf/emnlp/ShengGCYWLX20} employ the meta-learning framework and generalize the relational information from few-shot examples to all the link prediction (LP) queries. Based on these methods, Mirtaheri et al. \cite{DBLP:conf/akbc/MirtaheriR0MG21} develop a method aiming to alleviate these problems for TKGs by considering temporal dependencies between facts. They formulate the one-shot extrapolated LP task for TKGs and propose new datasets based on benchmark TKG databases ICEWS \cite{DVN/28075_2015} and GDELT \cite{leetaru2013gdelt}.

In this paper, we extend the existing TKG LP tasks, i.e., interpolated and extrapolated LP, to the one-shot setting, and we propose a model learning \textbf{m}eta representations of \textbf{o}ne-\textbf{s}ho\textbf{t} relations for solving both tasks in TKGs (MOST).
The main contribution of our work is three folded:
\begin{itemize}
\item We propose four new large-scale datasets for one-shot relational learning in TKGs. For every long-tail relation, we have a substantial number of associated TKG facts, which greatly alleviates instability in model training and evaluation.
\item We fix the unfair evaluation settings employed by previous KG FSL methods, and report comprehensive evaluation results on our datasets.
\item We propose a model solving both interpolated and extrapolated LP tasks for one-shot relations in TKGs, and evaluate our model on all four newly-proposed datasets. Our model achieves state-of-the-art performance on all datasets in both LP tasks while keeping a low time cost.
\end{itemize}

\section{Related Work}
KG representation learning has brought success to KG models in recent years. A line of methods are translational models \cite{DBLP:conf/nips/BordesUGWY13,DBLP:conf/aaai/LinLSLZ15,DBLP:conf/iclr/SunDNT19,DBLP:conf/nips/AbboudCLS20}, while another series of methods are based on tensor factorization \cite{DBLP:conf/icml/NickelTK11,DBLP:journals/corr/YangYHGD14a,DBLP:conf/emnlp/BalazevicAH19}. On top of them, neural-based models \cite{DBLP:conf/esws/SchlichtkrullKB18,DBLP:conf/iclr/VashishthSNT20} incorporates deep learning into these two lines of models and achieve strong performance.
To further model temporal dynamics of TKGs, a number of studies have been conducted \cite{DBLP:conf/www/LeblayC18,Ma2019EmbeddingMF,DBLP:conf/iclr/LacroixOU20,jin-etal-2020-recurrent,DBLP:conf/emnlp/WuCCH20,DBLP:conf/emnlp/HanDMGT21,DBLP:conf/iclr/HanCMT21,DBLP:conf/kdd/JungJK21,DBLP:journals/corr/abs-2112-07791}, and have shown great effectiveness.

Xiong et al. \cite{DBLP:conf/emnlp/XiongYCGW18} find that data scarcity problem exists in KGs and traditional KG representation learning methods fail to model sparse KG relations. To solve this problem, several researches \cite{DBLP:conf/emnlp/XiongYCGW18,DBLP:conf/emnlp/ChenZZCC19,DBLP:conf/aaai/ZhangYHJLC20,DBLP:conf/emnlp/ShengGCYWLX20,DBLP:conf/sigir/NiuLTGDLWSHS21} employ FSL paradigm \cite{DBLP:conf/nips/VinyalsBLKW16,DBLP:conf/cvpr/SungYZXTH18} and make improvement in learning sparse relations. Mirtaheri et al. \cite{DBLP:conf/akbc/MirtaheriR0MG21} find the same problem in TKGs and develop a one-shot TKG reasoning model aiming to better model sparse relations in TKGs. Please refer to Appendix \ref{appendix: related work} for more details.

\section{Task Formulation}
Let $\mathcal{E}$, $\mathcal{R}$ and $\mathcal{T}$ represent a finite set of entities, relations and timestamps, respectively. A temporal knowledge graph $\mathcal{G}$ is a relational graph consisting of a finite set of facts denoted with quadruples in the form of $(s, r, o, t)$, i.e., $\mathcal{G} = \{(s,r,o,t)|s, o\in \mathcal{E}, r \in \mathcal{R}, t \in \mathcal{T}\}$.
We define the TKG link prediction (LP) tasks as follows. A complete TKG contains both the observed facts and the unobserved true facts, i.e., $\mathcal{G} = \{\mathcal{G}_{obs}, \mathcal{G}_{un}\}$. Given $\mathcal{G}_{obs}$, TKG LP aims to predict the ground truth object (or subject) entities of queries $(s, r, ?, t)$ (or $(?, r, o, t)$), where $(s, r, o, t)\in \mathcal{G}_{un}$. For any link prediction query $(s, r, ?, t_q)$ (or $(?, r, o, t_q)$), while in the interpolation task the prediction can be based on all the observed facts $\{(s,r,o,t_i)|t_i\in \mathcal{T}\}$ from any timestamp, the extrapolated LP task further regulates that the prediction can only be based on the observed past facts $\{(s,r,o,t_i)|t_i<t_q\}$. In our work, we only consider object prediction since we add reciprocal relations
for every quadruple , i.e., adding $(o,r^{-1},s,t)$ for
every $(s,r,o,t)$. The restriction to only predict
object entities does not lead to a loss of generality.
\subsection{One-shot Temporal Knowledge Graph Link Prediction Setup}
TKG LP can be generalized to the one-shot setting as follows. Given a TKG $\mathcal{G}$, all its relations $\mathcal{R}$ can be classified into two groups, i.e., frequent relations $\mathcal{R}_{freq}$ and sparse relations $\mathcal{R}_{sp}$. Following the standard KG one-shot learning setting \cite{DBLP:conf/emnlp/XiongYCGW18,DBLP:conf/akbc/MirtaheriR0MG21}, we assume that we have access to a set of training tasks. Each training task $T_r$ corresponds to a sparse KG relation $r\in \mathcal{R}_{sp}^{train}$ ($\mathcal{R}_{sp}^{train} \subset \mathcal{R}_{sp}$). $T_r = \{\mathcal{S}_r, \mathcal{Q}_r\}$, where $\mathcal{S}_r = \{(s_0,r,o_0,t_0)\}$ is a set containing only one support quadruple $(s_0,r,o_0,t_0)$, and $\mathcal{Q}_r = \{(s_q, r, o_q, t_q)|s_q,o_q\in \mathcal{E}, t_q \in \mathcal{T}\} \setminus \{(s_0,r,o_0,t_0)\}$ is a set containing the rest of quadruples concerning the sparse relation $r$. We name $\mathcal{S}_r$ and $\mathcal{Q}_r$ as the support set and the query set of the task $T_r$, respectively. And the set of all training tasks is also denoted as the meta-training set $\mathbb{T}_{meta-train}$ ($\mathbb{T}_{meta-train} = \{T_r|r \in \mathcal{R}_{sp}^{train}\}$). For every sparse relation $r$, the goal of one-shot TKG LP is to accurately predict the missing entities of all the LP queries $(s_q, r, ?, t_q)$ (or $(?, r, o_q, t_q)$) derived from query quadruples $(s_q, r, o_q, t_q)\in \mathcal{Q}_r$, with only one observed $r$-specific support quadruple $(s_0, r, o_0, t_0)$ from $\mathcal{S}_r$. After the meta-training process, TKG reasoning models will be evaluated on meta-test tasks $\mathbb{T}_{meta-test}$ corresponding to unseen sparse relations $\mathcal{R}_{sp}^{test}$, where $\mathcal{R}_{sp}^{test} \subset \mathcal{R}_{sp}$ and $\mathcal{R}_{sp}^{train} \cap \mathcal{R}_{sp}^{test} = \phi$. We also validate the model performance with a meta-validation set $\mathbb{T}_{meta-valid}$ ($\mathcal{R}_{sp}^{valid} \subset \mathcal{R}_{sp}$, $\mathcal{R}_{sp}^{train} \cap \mathcal{R}_{sp}^{valid} = \phi, \mathcal{R}_{sp}^{valid} \cap \mathcal{R}_{sp}^{test} = \phi$). Similar to meta-training, for each sparse relation in meta-validation and meta-test, only one associated quadruple is observed in its support set, and all the links in its query set are to be predicted. Besides, same as \cite{DBLP:conf/emnlp/XiongYCGW18,DBLP:conf/akbc/MirtaheriR0MG21}, a background graph $\mathcal{G}' = \{(s,r,o,t)|s, o\in \mathcal{E}, r \in \mathcal{R}_{freq}, t \in \mathcal{T}\}$ is constructed by including all the quadruples concerning frequent relations, and it is also observable to the TKG reasoning model.

For each sparse relation $r$, in the interpolated LP, there is no constraint for the timestamp $t_0$ of its support quadruple $(s_0, r, o_0, t_0)$, while in the extrapolated LP, temporal constraint is imposed that $t_0 < \text{min}(\{t_q|(s_q,r,o_q,t_q) \in \mathcal{Q}_r\})$. Moreover, in the interpolated LP, we assume that the whole background graph is always observable, while in the extrapolated LP, only the background graph before the one-shot support quadruple is observable. Following the extrapolation setting in \cite{DBLP:conf/akbc/MirtaheriR0MG21}, we keep the time span of meta-learning sets ($\mathbb{T}_{meta-train}$, $\mathbb{T}_{meta-valid}$, $\mathbb{T}_{meta-test}$) in a non-overlapped sequential order (Figure \ref{fig: ext_setting3}).

\section{Our Method}
\begin{figure*}[htbp]
\centering
\includegraphics[scale=0.35]{pics/model3.png}
\caption{\label{fig: model} Overview of MOST. Assume we have a sparse relation \textit{Appeal for change in leadership} and the associated support quadruple is (\textit{Iraq}, \textit{Appeal for change in leadership}, \textit{Nuri al-Maliki}, 2013-01-15) $\in \mathcal{S}_r$. MOST first finds the nearest observable temporal neighbors of support entities, i.e., \textit{Iraq} and \textit{Nuri al-Maliki}, and computes their contextualized time-aware representations. Concatenated support entity representations are then through a multilayer perceptron, which generates a meta-representation for the sparse relation. Assume we have a query quadruple (\textit{Protestor (Egypt)}, \textit{Appeal for change in leadership}, \textit{Head of Government (Egypt)}, 2013-01-23) $\in \mathcal{Q}_r$ and it derives a link prediction query (\textit{Protestor (Egypt)}, \textit{Appeal for change in leadership}, $?$, 2013-01-23). Candidate entity representations are updated with a time encoder and the scoring function will compute the plausibility scores of the quadruples with different candidates.}
\end{figure*}
Figure \ref{fig: model} shows the overview of MOST. MOST consists of two components: (1) Time-aware relational graph encoder for learning contextualized time-aware entity representations; (2) Meta-relational decoder employed to generate meta representations for sparse relations and compute the plausibility scores of the quadruples.
\subsection{Time-aware Relational Graph Encoder}
MOST employs a time-aware relational graph encoder to learn the contextualized entity representations of support entities. For every support entity $e$, MOST first finds $e$'s temporal neighbors. It searches for the background facts whose object entity corresponds to $e$, and constructs a temporal neighborhood, i.e., $\mathcal{N}_{e} = \{(e',r', t')|r' \in \mathcal{R}_{freq}, (e', r', e, t') \in \mathcal{G}'\}$.
To avoid including excessive noisy information, MOST filters these background facts and only keeps a fixed number of temporal neighbors nearest to the timestamp of the support quadruple $t_0$ (MOST ensures the kept neighbors are prior to $t_0$ for extrapolated LP). The number of sampled neighbors is a hyperparameter and can be tuned. We denote the filtered neighborhood as $\Tilde{\mathcal{N}}_e$.

MOST then computes $e$'s time-aware representation by aggregating the information provided by its temporal neighbors. Inspired by \cite{DBLP:conf/iclr/VashishthSNT20}, our relational graph encoder derives the entity representations as follows:
\begin{equation}
\label{eq: agg}
\begin{aligned}
    &\Tilde{\mathbf{h}}_{e}^{l+1} = \frac{1}{|\Tilde{\mathcal{N}}_{e}|} \sum_{(e',r', t') \in \Tilde{\mathcal{N}}_{e}}\mathbf{W}^l (f(\mathbf{h}_{e'}^{l}\|\boldsymbol{\Phi}(t'))\circ\mathbf{h}_{r'}),\\
    &\mathbf{h}_{e}^{l+1} = \mathbf{h}_{e}^{l}+ \delta_1 \sigma(\Tilde{\mathbf{h}}_{e}^{l+1}).
\end{aligned}
\end{equation}
$\Tilde{\mathbf{h}}_e^{l+1}$ denotes $e$'s entity representation in the ${(l+1)}^{th}$ layer. $\mathbf{h}_{r'}$ denotes the relation representation of the frequent relation $r'$. Moreover, $\circ$, $\|$ represent Hadamard product and concatenation operation, respectively. $\mathbf{W}^l$ is a weight matrix that processes the information in the $l^{th}$ layer. $f$ is a layer of feed-forward neural network. $\delta_1$ is a trainable parameter deciding how much temporal information is included in updating entity representation. $\sigma(\cdot)$ is an activation function. $\boldsymbol{\Phi}(t')$ denotes the time encoding function that encodes timestamp $t'$ as: $\sqrt{\frac{1}{d_t}}[cos(\omega_1t'+\phi_1), \dots, cos(\omega_{d_t}t'+ \phi_{d_t}))]$, where $d_t$ is the dimension of the time representation, $\omega_1 \dots \omega_{d_t}$ and $\phi_1 \dots \phi_{d_t}$ are trainable parameters. We name our model with this timestamp encoder as MOST-TA. Due to the success of \cite{DBLP:conf/iclr/XuRKKA20} who encodes time differences for dynamic graphs, we input $t_0-t'$ instead of $t'$ into the time encoder $\boldsymbol{\Phi}(\cdot)$ and derive another model variant MOAT-TD. We show in experiments (Section \ref{exp results}) that both MOST-TA and MOST-TD can achieve state-of-the-art performance in one-shot TKG LP tasks.

\subsection{Meta-relational Decoder}
Assume we have a meta-learning task $T_r$, and we have a support quadruple $(s_0, r, o_0, t_0) \in \mathcal{S}_r$, with several query quadruples $(s_q, r, o_q, t_q) \in \mathcal{Q}_r$.
After obtaining the time-aware representations for support entities, i.e., $\mathbf{h}_{s_0}$ and $\mathbf{h}_{o_0}$, we learn a meta representation for the sparse relation $r$:
\begin{equation}
\label{eq: relproj}
    \mathbf{h}_{r} = \text{MLP}(\mathbf{h}_{s_0}\|\mathbf{h}_{o_0}),
\end{equation}
where \text{MLP} is a multilayer perceptron consisting of three layers of feed-forward neural network.

For a link prediction query $(s_q, r, ?, t_q)$ (generated from a query quadruple $(s_q, r, o_q, t_q)$), we compute the time-aware representations of $s_q$ as well as every candidate entity as follows:
\begin{equation}
\label{eq: queryemb}
    \mathbf{h}_{e}^{l+1} = \mathbf{h}_{e}^l+ \delta_2 f(\mathbf{h}_{e}^l\|\boldsymbol{\Phi}(t_q)).
\end{equation}
$\delta_2$ is a trainable parameter controlling the amount of the injected temporal information. We do not search temporal neighbors for every candidate entity to avoid huge time consumption. Similarly, MOST-TD adapts Equation \ref{eq: queryemb} to the following form to enable time difference learning:
\begin{equation}
\label{eq: queryemb2}
    \mathbf{h}_{e}^{l+1} = \mathbf{h}_{e}^l+ \delta_2 f(\mathbf{h}_{e}^l\|\boldsymbol{\Phi}(t_q-t_0)).
\end{equation}

Inspired by \cite{DBLP:conf/iclr/SunDNT19}, we map the entity representations to the complex space, and treat the meta representations of sparse relations as element-wise rotation.
We use the following scoring function to compute the scores regarding different candidates $e_{c}$:
\begin{equation}
\label{equa: scorefunc}

    \Psi_{e_c}^r = \frac{1}{1+exp\left( - \left(\mathbf{h}_{s_q} \circ \Tilde{\mathbf{h}}_r \right)^T \mathbf{h}_{e_{c}} \right) }

\end{equation}
$\Tilde{\mathbf{h}}_r = \frac{1}{\lVert\mathbf{h}_r\rVert_\infty} \mathbf{h}_r$, and $\lVert\mathbf{h}_r\rVert_\infty$ denotes the infinity norm of the vector $\mathbf{h}_r$.
In our framework, the meta representations of sparse relations are computed from the
one-shot support entity representations, thus naturally coupled with noise. We adaptively regularize the lengths of meta representations by dividing themselves by their infinity norm. We do not divide by L2 norm since we argue that both the directions and the lengths of meta representations contribute to prediction. Diving by L2 norm will lose length information. We name this process as norm regularization.

\subsection{Parameter Learning}
We use the binary cross entropy loss to train our model.
The loss function in our work is in the following form:
\begin{equation}
\label{eq: loss}
\begin{aligned}
    l_{q}^{e_c} = & \ y_{q}^{e_c}\textrm{log}\Psi_{e_c}^r + (1-y_{q}^{e_c})(1-\textrm{log}\Psi_{e_c}^r),\\
    \mathcal{L} = & -\sum \limits_{r \in \mathcal{R}_{sp}^{train}} \frac{1}{2|\mathcal{Q}_r|} \sum \limits_{q} \frac{1}{|\mathcal{E}|} \sum \limits_{e_c \in \mathcal{E}} l_{q}^{e_c}.

\end{aligned}
\end{equation}
$q$ denotes a LP query derived from the query quadruple $(s_q, r, o_q, t_q) \in \mathcal{Q}_r$ corresponding to a sparse relation $r \in \mathcal{R}_{sp}^{train}$. $y_{q}^{e_c}$ is the binary label indicating whether a candidate entity $e_c$ is the ground truth missing entity from $q$ or not. Note that in Equation \ref{eq: loss}, for every query quadruple $q$, we average the loss over all the entities $e_{c} \in \mathcal{E}$ other than $o_q$ (or $s_q$), rather than sampling a fixed number of negative samples as in previous KG FSL methods \cite{DBLP:conf/emnlp/XiongYCGW18,DBLP:conf/emnlp/ChenZZCC19,DBLP:conf/aaai/ZhangYHJLC20,DBLP:conf/sigir/NiuLTGDLWSHS21,DBLP:conf/akbc/MirtaheriR0MG21}.

\section{Experiments}
We evaluate MOST on both interpolated and extrapolated TKG LP in the one-shot setting. We propose four new large-scale datasets and evaluate previous KG FSL methods on them.

\subsection{Datasets}
\begin{table}[htbp]
    \centering
    \resizebox{\columnwidth}{!}{
\begin{tabular}{c c c c c} \hline
Dataset&$|\mathcal{E}|$&$|\mathcal{R}|$&$|\mathcal{T}|$&$|\mathbb{T}|$\\ \hline
ICEWS05-15-one\_ext  & $7,934$ & $109$ & $4,017$ & $53/6/11$  \\
ICEWS05-15-one\_int  & $10,356$ & $155$ & $4,017$ & $74/9/10$  \\
GDELT-one\_ext  & $6,647$ & $155$ & $2,751$ & $55/7/11$  \\
GDELT-one\_int  & $7,677$ & $181$ & $2,751$ & $64/8/8$ \\
\hline
\end{tabular}}
\caption{Dataset statistics. $|\mathbb{T}|$ denotes the number of meta-learning tasks in $\mathbb{T}_{meta-training}$, $\mathbb{T}_{meta-valid}$, $\mathbb{T}_{meta-test}$.}
\label{tab: data}
\end{table}
By taking subsets of two benchmark TKG databases, i.e., ICEWS \cite{DVN/28075_2015} and GDELT \cite{leetaru2013gdelt}, Mirtaheri et al. \cite{DBLP:conf/akbc/MirtaheriR0MG21} propose two one-shot extrapolation LP datasets, i.e., ICEWS17 and GDELT. They first set upper and lower thresholds, and then select the relations with frequency between them as sparse relations (frequency 50 to 500 for ICEWS17, 50 to 700 for GDELT). To prevent time overlaps among meta-learning sets ($\mathbb{T}_{meta-train}$, $\mathbb{T}_{meta-valid}$, $\mathbb{T}_{meta-test}$), they further remove a significant number of quadruples regarding sparse relations. For example, a relation $r$ is selected as a sparse relation and $T_r \in \mathbb{T}_{meta-train}$. The ending timestamp of the meta-training set is $t'$. Then all the quadruples in $\{(s,r,o,t)|s,o \in \mathcal{E}, t> t'\}$ are removed from the dataset. This leads to a considerably smaller query set $\mathcal{Q}_r$ when a large number of $r$-related events take place after $t'$. If $r$'s frequency is close to the lower threshold before removal, it is very likely that after removal, the number of associated quadruples left in $\{(s,r,o,t)|s,o \in \mathcal{E}, t\leq t'\}$ becomes extremely small, leading to a tiny query set $\mathcal{Q}_r$ that causes instability during training. Similarly, if $T_r \in \mathbb{T}_{meta-valid}$ or $T_r \in \mathbb{T}_{meta-test}$, evaluation on a tiny query set $\mathcal{Q}_r$ makes it hard to accurately determine the performance of the model since the test data is not comprehensive. As shown in Figure \ref{fig: data comparison}, a large portion of sparse relations in ICEWS17 and GDELT have very few associated quadruples, which introduces instability in model training and evaluation.

To overcome this problem, we construct two new large-scale extrapolation LP datasets, i.e., ICEWS-one\_ext and GDELT-one\_ext, also by taking subsets from ICEWS \cite{DVN/28075_2015} and GDELT \cite{leetaru2013gdelt}. ICEWS-one\_ext contains timestamped political facts happening from 2005 to 2015, while GDELT-one\_ext contains global social facts from Jan. 1, 2018 to Jan. 31, 2018. For sparse relation selection, we set the upper and lower thresholds of frequency to 100 and 1000 for ICEWS-one\_ext, 200 and 2000 for GDELT-one\_ext, and then split these relations into train/valid/test groups. We take the relations with higher frequency as frequent relations $\mathcal{R}_{freq}$ and build background graphs $\mathcal{G}'$ with all the quadruples containing them. Following \cite{DBLP:conf/akbc/MirtaheriR0MG21}, we then remove a part of quadruples associated with sparse relations to prevent time overlaps among meta-learning sets. After removal, we further discard the relations with too few associated quadruples (less than 50 for ICEWS-one\_ext, 100 for GDELT-one\_ext). In this way, we prevent including meta-tasks $T_r$ with extremely small query set $\mathcal{Q}_r$.
From Figure \ref{fig: data comparison}, we observe that ICEWS-one\_ext and GDELT-one\_ext have a substantial number of associated quadruples for each sparse relation, which greatly alleviates instability in model training and evaluation.
\begin{figure}[htbp]
  \begin{subfigure}[htbp]{\columnwidth}
    \includegraphics[width=\columnwidth]{pics/icews_comparison.png}

    \label{fig: icews_comparison}
  \end{subfigure}
  \begin{subfigure}[htbp]{\columnwidth}
    \includegraphics[width=\columnwidth]{pics/gdelt_comparison.png}

    \label{fig: gdelt_comparison}
  \end{subfigure}
  \caption{Sparse Relation frequency comparison between ICEWS-one\_ext and ICEWS17; GDELT-one\_ext and GDELT.}
  \label{fig: data comparison}
\end{figure}

Similarly, we construct two more datasets, i.e., ICEWS-one\_int and GDELT-one\_int, for interpolated LP in the one-shot setting. Since in interpolated LP, there exists no constraint on the support timestamp $t_0$, we do not have to remove quadruples to eliminate time overlaps among meta-learning sets. To this end, we set the upper and lower thresholds of sparse relations' frequency to 50 and 500 for ICEWS-one\_int, 100 and 1000 for GDELT-one\_int, and then split these relations into train/valid/test groups. Statistics of our datasets are presented in Table \ref{tab: data} and we provide more details about data construction in Appendix \ref{appendix: data construct}.

\subsection{Baseline Methods}
We compare our model with two groups of baseline methods on both interpolated and extrapolated LP in the one-shot setting.
\subsubsection*{Few-shot Relational Learning Methods}
We consider five static KG FSL methods, i.e., Gmatching \cite{DBLP:conf/emnlp/XiongYCGW18}, MetaR \cite{DBLP:conf/emnlp/ChenZZCC19}, FSRL \cite{DBLP:conf/aaai/ZhangYHJLC20}, FAAN \cite{DBLP:conf/emnlp/ShengGCYWLX20}, GANA \cite{DBLP:conf/sigir/NiuLTGDLWSHS21}, and a TKG FSL method, i.e., OAT \cite{DBLP:conf/akbc/MirtaheriR0MG21}. In \cite{DBLP:conf/akbc/MirtaheriR0MG21}, static KG FSL methods are trained and evaluated on an unweighted static KG derived from collapsing the original TKG, which greatly decreases the inductive bias brougt by the original TKG and causes poor performance of these methods. In our work, we provide static KG FSL methods with all the facts in the original datasets, and neglect time information, i.e., neglecting $t$ in $(s,r,o,t)$. We provide a detailed explanation and prove our assertion empirically with further experiments in Appendix \ref{appendix: unweighted}.
\subsubsection*{Temporal Knowledge Graph Embedding Methods}
Three TKG interpolation methods, i.e., TNTComplEx \cite{DBLP:conf/iclr/LacroixOU20}, ATiSE \cite{DBLP:journals/corr/abs-1911-07893}, TeLM \cite{DBLP:conf/naacl/XuCNL21}, and three TKG extrapolation methods, i.e., TANGO \cite{DBLP:conf/emnlp/HanDMGT21}, CyGNet \cite{DBLP:conf/aaai/ZhuCFCZ21}, xERTE \cite{DBLP:conf/iclr/HanCMT21} are selected as our baselines. For each interpolation dataset, we build a training set for these methods by adding all the quadruples of the background graph $\mathcal{G}'$ and the quadruples associated with all the meta-training relations $\mathcal{R}_{sp}^{train}$. We further add the support quadruple associated with each sparse relation $r \in \{\mathcal{R}_{sp}^{valid}, \mathcal{R}_{sp}^{test}\}$ into the training set. For each extrapolation dataset, we build a training set by adding all the background quadruples during meta-training time $\mathcal{G}'_{train}$, as well as all the quadruples concerning every $r \in \mathcal{R}_{sp}^{train}$. We do not include any quadruple regarding $r \in \{\mathcal{R}_{sp}^{valid}, \mathcal{R}_{sp}^{test}\}$ into the training set due to the time constraint in the extrapolation setting. But we allow the models to have access to the support quadruples ($\mathcal{S}_r$, $r \in \{\mathcal{R}_{sp}^{valid}, \mathcal{R}_{sp}^{test}\}$) during inference. We test TKG embedding baselines with the same quadruples tested by FSL methods to ensure fair comparison.
\begin{table*}[t]
    \centering
    \resizebox{\textwidth}{!}{
    \large\begin{tabular}{@{}|l|ccccc|ccccc|@{}}
\toprule
        Datasets & \multicolumn{5}{|c}{\textbf{ICEWS-one\_int}} & \multicolumn{5}{|c|}{\textbf{GDELT-one\_int}}\\
\midrule
        Model & MRR & Hits@1 & Hits@3 & Hits@5 & Hits@10 & MRR & Hits@1 & Hits@3 & Hits@5 & Hits@10 \\
\midrule

        TNTComplEx & 23.34 & 14.57 & 27.21 & 31.54 & 36.88
        & 11.95 & 6.76 & 11.98 & 15.58 & 21.78
         \\
        ATiSE & 34.40 & 22.03 & 39.51 & 49.25 & 60.57
        & 7.77 & 5.10 & 6.72 & 8.13 & 12.13
         \\
        TeLM  & 35.38 & 24.42 & 39.21 & 47.74 & 59.12
        & 10.41 & 5.97 & 10.37 & 13.28 & 18.87
         \\

\midrule
        GANA  & 13.83 & 6.07 & 19.21 & 24.00 & 27.23
       & 5.89 & 2.53 & 6.54  & 8.35 & 12.20
       \\
       MetaR  & 27.69 & 7.88 & 41.02 & 52.58 & 61.78
       & 9.91 & 0.18 & 14.61 & 19.62 & 26.79
       \\
         GMatching & 30.59 & 15.46 & 39.33 & 48.80 & 58.62
         & 12.53 & 6.55 & 12.80 & 17.14 & 24.15
         \\

    FSRL & 33.98 & 18.94 & 44.48 & 52.61 & 59.82
       & 14.11 & 7.61 & 14.67 & 19.56 & 27.54
       \\

       FAAN & 35.48 & 23.27 & 43.36 & 49.45 & 57.73
       & 14.77 & 7.67 & 16.19 & 21.35 & 27.11
        \\
\midrule
      OAT  & 11.55 & 5.47 & 10.09 & 14.81 & 23.42
      & 12.28 & 7.70 & 12.59 & 15.18 & 21.47
        \\
\midrule
        MOST-TA
        & \textbf{47.79} & \textbf{39.91} & 51.79 & \textbf{57.01} & 62.25
        & \textbf{17.71} & 11.56 & \textbf{19.07} & \textbf{23.25} & \textbf{29.76}
        \\
        MOST-TD
        & 47.60 & 39.43 & \textbf{51.98} & 56.83 & \textbf{62.38}
        & 17.36 & \textbf{11.67} & 18.18 & 22.74 & 28.63
        \\
\bottomrule
    \end{tabular}
    }
    \caption{Interpolated LP results for one-shot relational learning on ICEWS-one\_int and GDELT-one\_int. Evaluation metrics are filtered MRR and Hits@1/3/5/10.
    The best results are marked in bold. More discussions in Appendix \ref{app: model analysis}.}\label{tab: tid int results}
\end{table*}
\begin{table*}[t]
    \centering
    \resizebox{\textwidth}{!}{
    \large\begin{tabular}{@{}|l|ccccc|ccccc|@{}}
\toprule
        Datasets & \multicolumn{5}{|c}{\textbf{ICEWS-one\_ext}} & \multicolumn{5}{|c|}{\textbf{GDELT-one\_ext}}\\
\midrule
        Model & MRR & Hits@1 & Hits@3 & Hits@5 & Hits@10 & MRR & Hits@1 & Hits@3 & Hits@5 & Hits@10 \\
\midrule
        TANGO & 10.23 & 3.94 & 11.40 & 15.88 & 25.78
        & 13.88 & 9.61 & 13.17 & 16.93 & 22.29
         \\
        CyGNet & 22.30 & 12.61 & 25.51 & 30.46 & 39.13
        & 9.42 & 4.87 & 9.74 & 13.13 & 16.81
         \\
        xERTE  & 30.02 & 19.79 & 36.63 & 42.13 & 51.16
        & 16.38 & 10.88 & 18.23 & 22.19 & 27.76
         \\
\midrule
        GANA  & 11.34 & 3.70 & 15.52 & 19.25 & 25.67
       & 7.12 & 4.85 & 7.02  & 8.89 & 11.13
       \\
       MetaR  & 23.50 & 9.01 & 32.95 & 40.18 & 48.73
       & 9.66 & 0.03 & 13.79  & 19.52 & 26.30
       \\
         GMatching & 20.30 & 12.35 & 21.06 & 28.80 & 38.02
         & 12.26 & 8.41 & 11.44 & 13.76 & 19.01
         \\

    FSRL & 18.06 & 12.09 & 17.68 & 21.06 & 32.23
       & 6.96 & 2.52 & 8.81 & 11.58 & 14.13
       \\
       FAAN & 25.73 & 15.86 & 29.14 & 35.95 & 43.73
       & 14.36 & 8.71 & 15.31 & 18.46 & 23.71
        \\
\midrule
      OAT  & 9.24 & 7.02 & 7.31 & 7.40 & 11.88
      & 14.06 & 6.71 & 13.43 & 18.59 & 28.11
        \\
\midrule
        MOST-TA
        & 32.94 & 26.35 & 34.64 & 39.97 & 47.19
        & 15.69 & 10.14 & 16.49 & 20.54 & 26.38
        \\
        MOST-TD
        & \textbf{38.46} & \textbf{31.51} & \textbf{40.73} & \textbf{46.02} & \textbf{52.32}
        & \textbf{17.36} & \textbf{11.64} & \textbf{18.37} & \textbf{22.46} & \textbf{28.15}
        \\
\bottomrule
    \end{tabular}
    }
    \caption{Extrapolated LP results for one-shot relational learning on ICEWS-one\_ext and GDELT-one\_ext. Evaluation metrics are filtered MRR and Hits@1/3/5/10.
    The best results are marked in bold. More discussions in Appendix \ref{app: model analysis}.}\label{tab: tid ext results}
\end{table*}

\subsection{Experimental Results}
\label{exp results}
Previous KG FSL methods only report object prediction results. To achieve comprehensive results, for each test quadruple $(s_q, r_q, o_q, t_q)$, we derive two LP queries, i.e., $(s_q, r_q, ?, t_q)$ and $(?, r_q, o_q, t_q)$, and perform prediction on both of them. We employ two evaluation metrics, i.e., Hits@1/3/5/10 and mean reciprocal rank (MRR), to evaluate model performance in all experiments. We also follow \cite{DBLP:conf/nips/BordesUGWY13} and use filtered results for fairer evaluation. More details of evaluation protocol are provided in Appendix \ref{appendix: eval protocol}.

Table \ref{tab: tid int results} and Table \ref{tab: tid ext results} report the experimental results of one-shot interpolated and extrapolated LP, respectively. We find that static KG FSL methods can achieve competitive or even better performance compared with traditional TKG embedding methods, implying the effectiveness of FSL in modeling sparse relations in KGs. MOST outperforms baseline methods on all datasets in both LP tasks. While MOST-TA performs better than MOST-TD in the interpolation task, MOST-TD outperforms MOST-TA in predicting future facts. This can be explained as follows.
In the interpolated LP, part of the TKG at every timestamp is observable during training, thus enabling the time encoder to learn from all the timestamps. However, in the extrapolated LP, meta-training set does not span across the whole timeline, which leads to the degenerated model performance during inference when we sample the temporal neighbors from the timestamps unseen in the meta-training set. For extrapolation, modeling time differences achieves better results since almost all time differences we encounter during inference are already seen and learned during meta-training.
\begin{table}[t]
    \centering
    \resizebox{\columnwidth}{!}{
    \large\begin{tabular}{@{}|l|ccc|ccc|@{}}
\toprule
        Datasets & \multicolumn{3}{|c}{\textbf{ICEWS-one\_int}} &  \multicolumn{3}{|c|}{\textbf{ICEWS-one\_ext}}\\
\midrule
        Variants & MRR & Hits@1 & Hits@10 & MRR & Hits@1 & Hits@10  \\
\midrule
        A1  & 46.22 & 37.04  & 62.20 & 31.75  & 23.85 & 45.55
         \\
        A2  & 45.85 & 37.22  & 61.88 & 32.03 & 24.36 & 47.16
         \\
\midrule
        B1 & 46.15 & 37.86  & 62.23 & 31.66 & 23.05 & 47.07
         \\
        B2  & 11.75 & 6.14  & 23.62 & 11.40 & 4.86 & 25.80
         \\
         B3  & 16.27 & 7.36  & 32.47 & 26.13 & 17.97 & 43.65
         \\
\midrule
         C1  & 43.15 & 32.81  & 62.22 & 31.77 & 24.74 & 44.16
         \\
         C2  & 42.81 & 31.52  & 62.10 & 32.34 & 26.01 & 45.60
         \\
\midrule
        MOST-TA &  \textbf{47.79} & \textbf{39.91} & \textbf{62.25} & \textbf{32.94}
        & \textbf{26.35} & \textbf{47.19}
        \\
\bottomrule
    \end{tabular}
    }
    \caption{Ablation studies of MOST-TA variants on ICEWS-one\_int and ICEWS-one\_ext. The best results are marked in bold.}\label{tab: ablation most-ta}
\end{table}
\begin{table}[t]
    \centering
    \resizebox{\columnwidth}{!}{
    \large\begin{tabular}{@{}|l|ccc|ccc|@{}}
\toprule
        Datasets & \multicolumn{3}{|c}{\textbf{ICEWS-one\_int}} &  \multicolumn{3}{|c|}{\textbf{ICEWS-one\_ext}}\\
\midrule
        Variants & MRR & Hits@1 & Hits@10 & MRR & Hits@1 & Hits@10  \\
\midrule
        A1  & 44.01 & 35.40 & 60.89 & 35.96  & 28.63 & 49.83
         \\
        A2  & 43.77 & 33.81  & 62.25 & 35.89 & 29.01 & 50.25
         \\
\midrule
        B1 & 42.94 & 32.02  & 62.28 & 36.30 & 29.52 & 51.22
         \\
        B2  & 11.78 & 6.21  & 23.44 & 11.67 & 4.78 & 24.78
         \\
         B3  & 23.68 & 13.87  & 44.40 & 27.35 & 18.82 & 43.78
         \\
\midrule
         C1  & 41.06 & 29.04  & 62.25 & 36.81 & 30.24 & 50.16
         \\
         C2  & 43.33 & 32.96  & 62.36 & 36.84 & 29.52 & 50.67
         \\
\midrule
        MOST-TD &  \textbf{47.60} & \textbf{39.43} & \textbf{62.38} & \textbf{38.46}
        & \textbf{31.51} & \textbf{52.32}
        \\
\bottomrule
    \end{tabular}
    }
    \caption{Ablation studies of MOST-TD variants on ICEWS-one\_int and ICEWS-one\_ext. The best results are marked in bold.}\label{tab: ablation most-td}
\end{table}
\subsection{Ablation Study}
To validate the effectiveness of model components, we derive model variants for both MOST-TA and MOST-TD, and conduct several ablation studies on ICEWS-one\_int and ICEWS-one\_ext. We present the experimental results in Table \ref{tab: ablation most-ta} and Table \ref{tab: ablation most-td}. We devise model variants from the following angles:

\textbf{(A) Temporal Neighbor Sampling Strategy:}
In A1, we keep all the temporal neighbors, without limiting the size of the sampled neighborhood. In A2, we still keep a fixed number of temporal neighbors, but we randomly sample them instead of keeping the nearest ones. We observe that by sampling nearest neighbors, we exclude excessive noise from farther neighbors and achieve better results.

\textbf{(B) Decoder Variants:} In B1, we switch the infinity norm to the L2 norm during norm regularization. In B2, we take the raw entity representations of query entities as the input into our decoder, without injecting temporal information through Equation \ref{eq: queryemb} and Equation \ref{eq: queryemb2}. In B3, we employ the LSTM-based matcher \cite{DBLP:journals/neco/HochreiterS97} proposed in \cite{DBLP:conf/emnlp/XiongYCGW18} instead of the scoring function (Equation \ref{equa: scorefunc}) in our meta-relational decoder. We observe that dividing by infinity norms helps to adaptively regularize the lengths of meta representations and improve model performance. Compared with \cite{DBLP:conf/emnlp/XiongYCGW18}, our meta-relational decoder employs a stronger scoring function. During score computation, introducing temporal information into the query entities greatly increases performance.

\textbf{(C) Graph Encoder Variants:}
In C1, we use RGCN \cite{DBLP:conf/esws/SchlichtkrullKB18} instead of our graph encoder. In C2, we remove time encoder $\boldsymbol{\Phi}(\cdot)$ during aggregation in our graph encoder. We observe that incorporating temporal information into the graph encoder helps to improve model performance.

\subsection{Performance over Different Relations}
We report the model performance over every sparse relation in ICEWS-one\_int (Table \ref{tab: diff rel int}) and ICEWS-one\_ext (Table \ref {tab: diff rel ext}). We compare MOST with the strongest FSL baseline FAAN. Experimental results show that MOST achieves performance gain in all sparse relations, which implies its robustness.
\begin{table}[htbp]
    \centering
    \resizebox{\columnwidth}{!}{
    \large\begin{tabular}{@{}|cc|cc|cc|@{}}
\toprule
        \multicolumn{2}{|c}{\textbf{ICEWS-one\_int}} & \multicolumn{2}{|c|}{MRR} &  \multicolumn{2}{c|}{Hits@10}\\
\midrule
        CAMEO & Frequency & MOST-TA & FAAN & MOST-TA & FAAN  \\
\midrule
        1313 & 65 & \textbf{35.25} & 28.13 & \textbf{57.03} & 46.09 \\
        1411 & 89 &  \textbf{53.65} & 40.86 & \textbf{71.02} &  67.61 \\
        035 & 92 & \textbf{50.50} & 36.38 & \textbf{71.42} & 63.73 \\
        151 & 118 & \textbf{75.75} & 53.06 & \textbf{85.04} & 82.47 \\
        023 & 146 & \textbf{52.69} & 49.37 & \textbf{64.48} & 63.10 \\
        191 & 175 & \textbf{63.83} & 49.41 & \textbf{75.57} & 75.28 \\
        1721 & 282 & \textbf{38.45} & 34.33 & 51.24 & \textbf{53.20} \\
        015 & 269 & \textbf{39.54} & 25.66 & \textbf{53.17} & 48.50 \\
        064 & 349 & \textbf{29.75} & 15.18 & \textbf{47.12} & 34.48 \\
        128 & 436 & \textbf{57.84} & 43.37 & 71.95 & \textbf{72.52} \\
\bottomrule
    \end{tabular}
    }
    \caption{Performance over each sparse relation in ICEWS-one\_int. CAMEO denotes the CAMEO code for each relation. The best results are marked in bold.}\label{tab: diff rel int}
\end{table}
\begin{table}[htbp]
    \centering
    \resizebox{\columnwidth}{!}{
    \large\begin{tabular}{@{}|cc|cc|cc|@{}}
\toprule
        \multicolumn{2}{|c}{\textbf{ICEWS-one\_ext}} & \multicolumn{2}{|c|}{MRR} &  \multicolumn{2}{c|}{Hits@10}\\
\midrule
        CAMEO & Frequency & MOST-TD & FAAN & MOST-TD & FAAN  \\
\midrule
        074 & 55 & \textbf{32.85} & 17.63 & \textbf{57.40} & 24.07 \\
        1122 & 57 & \textbf{21.39} & 4.99 & \textbf{29.46} & 13.39 \\
        0241 & 65 & \textbf{22.05} & 21.60 & \textbf{35.15} & 32.81 \\
        015 & 67 & \textbf{26.31} & 23.30 & \textbf{42.42} & 38.63 \\
        064 & 93 & \textbf{30.41} & 18.52 & \textbf{50.00} & 44.02 \\
        113 & 107 & \textbf{21.67} & 15.40 & \textbf{33.96} & 31.13 \\
        033 & 127 & \textbf{37.13} & 28.26 & \textbf{53.96} & 53.17 \\
        0214 & 128 & \textbf{42.91} & 43.88 & \textbf{55.11} & 52.36 \\
        072 & 130 & \textbf{33.51} & 9.74 & \textbf{47.28} & 32.17 \\
        154 & 180 & \textbf{59.61} & 32.96 & \textbf{72.06} & 65.64 \\
        115 & 184 & \textbf{47.47} & 36.28 & \textbf{61.20} & 54.09 \\

\bottomrule
    \end{tabular}
    }
    \caption{Performance over each sparse relation in ICEWS-one\_ext. CAMEO denotes the CAMEO code for each relation. The best results are marked in bold.}\label{tab: diff rel ext}
\end{table}

\subsection{Time Cost Analysis}
We report in Figure \ref{fig: time cost} and Table \ref{tab: eval time main body} the time cost of MOST and several strong baselines on both ICEWS-based datasets. We observe that MOST achieves the best performance on LP tasks while keeping a low time cost. Though MOST-TD and MOST-TA achieve weaker performance than their counterpart in the interpolated and extrapolated LP, respectively, they require much shorter training time and can still achieve superior performance as reported in Table \ref{tab: tid int results} and Table \ref{tab: tid ext results}. We provide more time cost details in Appendix \ref{appendix: time cost}.
\begin{figure}[htbp]
\centering
\includegraphics[scale=0.25]{pics/Runtime.png}
\caption{Training time comparison among MOST and the strongest baselines on ICEWS-based datasets.}
\label{fig: time cost}
\end{figure}
\begin{table}[htbp]
    \centering
    \resizebox{\columnwidth}{!}{
    \large\begin{tabular}{@{}lccccc@{}}
\toprule
        Model & MOST-TA & MOST-TD & FAAN & TeLM & xERTE\\
\midrule
        ICEWS-one\_int  & 0.10 & 0.11 & 35.93 & 0.02 & -
       \\
       ICEWS-one\_ext  & 0.20 & 0.23 & 27.20 & - & 5.23
       \\
\bottomrule
    \end{tabular}
    }
    \caption{Test time (min) comparison among MOST and the strongest baselines on ICEWS-based datasets.}\label{tab: eval time main body}
\end{table}
\section{Conclusion}
We extend both TKG interpolated and extrapolated LP tasks to the one-shot setting, and propose a model learning \textbf{m}eta representations of \textbf{o}ne-\textbf{s}ho\textbf{t} relations for solving both tasks (MOST). MOST learns meta representations from time-aware entity representations of the entities in the one-shot examples. It further employs a scoring function together with a norm regularizer for predicting missing entities. We propose four large-scale datasets, fix the unfair evaluation settings employed by previous KG FSL methods, and compare MOST with a large group of baselines on both TKG LP tasks. Experimental results show that MOST achieves stat-of-the-art performance on both one-shot LP tasks while keeping a low time cost.

\bibliographystyle{acl_natbib}


\section{Introduction}
Heat is often considered an energy waste, yet it can be converted into useful electrical energy using thermoelectric (TE) materials. However, compared to other energy conversion schemes, TE devices often have lower efficiency. One can assess the TE performance through the so-called figure of merit: $ZT = \frac{S^2\sigma}{\kappa}T$, where $S$ is the Seebeck coefficient, $\sigma$ is the electrical conductivity, $T$ is the operating temperature, and $\kappa$ is the total thermal conductivity. The total thermal conductivity, $\kappa$, is a sum of the electronic thermal conductivity ($\kappa_e$) and the lattice thermal conductivity ($\kappa_{ph}$), i.e., $\kappa = \kappa_e+\kappa_{ph}$.  The figure of merit is also proportional to the power factor (PF) by the following relation: $ZT = \mathrm{PF}. T / \kappa$, where $\mathrm{PF} = S^2 \sigma$.  From the expression of ZT, a good TE material should possess good electrical conduction and thermal isolation.  In other words, we should maximize the Seebeck coefficient and electrical conductivity, but, at the same time, minimize the thermal conductivity to obtain a good TE material. Unfortunately, the interrelation between those parameters makes it difficult to find materials with high ZT~\cite{majumdar2004thermoelectricity,vining2009inconvenient,snyder2011complex,sootsman2009new,zeier2016thinking,zebarjadi2012perspectives,zhang2015thermoelectric,zhu2017compromise,mao2018advances,he2017advances,yang2016tuning,gorai2017computationally,jain2016computational}. For example, based on the Wiedemann-Franz law, the ratio between electrical conductivity $\sigma$ and electrical thermal conductivity $\kappa_e$ is constant. Therefore, it is tricky to get material with high $\sigma$ simultaneously with low thermal conductivity $\kappa_e$~\cite{zhao2014ultralow,takahashi2012thermoelectric}.

To solve the TE trade-off problem, there have been some promising proposals, such as optimizing carrier concentration~\cite{CHEN2012535,Zhou2016}, using low-dimensional materials~\cite{hicks1993effect,hung2016quantum,dresselhaus2007new,urban2007synergism,lee2016thermal,eslamian2017inorganic,bottner2006aspects}, utilizing band energy convergence~\cite{pei2011convergence,zeier2014band,zhang2014high,gayner2016recent}, and considering hierarchical architecture~\cite{biswas2012high}. Two-dimensional (2D) materials such as 2D Dirac materials have been shown to feature good electronic and thermal properties which are potential for TE materials~\cite{hippalgaonkar2017high,duan2016high}. Moreover, a solely parabolic band in a 2D material has been found to give more desirable TE properties than a three-dimensional (3D) material due to the quantum confinement that can lead to an improvement in the figure of merit~\cite{hicks1993effect}.

The band engineering method such as tuning the band gap~\cite{pei2011alloying} or the effective mass~\cite{bilc2004resonant,heremans2012resonant} in terms of the curvature of the band also becomes an effective way to scan potential TE materials because this method uses a relatively cheap computational method by considering only the energy dispersion relation $E(k)$ and the scattering lifetime $\tau(E)$.  One can tune the bandgap or change the combination of the band structures to know the optimized structure that will give better TE properties.  Several works related to the band engineering have been performed for many types of band structures, such as pudding-mold bands~\cite{usui2017enhanced,usui2013large,kuroki2007pudding,isaacs2019remarkable,wei2020strain}, parabolic bands~\cite{xia2019leveraging,xia2019high,adhidewata2021thermoelectric}, and linear Dirac bands~\cite{hasdeo2019optimal}.  In particular, the existence of a heavy band alongside the Dirac bands gives a high power factor (PF) in a semimetal like \ch{CoSi}.  This result is contrary to the mainstream impression that semimetals are generally poor for thermoelectricity~\cite{xia2019high}.  Furthermore, the existence of a heavy parabolic band that intersects Dirac bands at the Fermi level is found to enhance PF due to the improved electron-phonon scattering in the form of a sharp spike of density of states (DOS)~\cite{xia2019leveraging}. Specifically, the existing heavy band acts as a filter for the low-energy carriers to be excited. This unbalance condition will boost the Seebeck coefficient.  However, rather than the heavy parabolic band, a single pudding-mold band can increase the electrical conductivity resulting in a higher power factor and figure of merit~\cite{usui2017enhanced}.  Regarding this fact, we are intrigued to ask, can the pudding-mold band alongside the Dirac bands lead to a better TE performance than the parabolic one?

To answer the question, in this paper, we will discuss the TE properties of 2D materials that possess Dirac bands with a heavy band. As we will see later, the existence of a heavy band alongside Dirac bands will improve the thermoelectric performance even with a narrow bandgap. We achieved the maximum value of $ZT$ when the heavy band intersects the Dirac bands at the Fermi level for the parabolic band but not for the other heavy bands~\cite{xia2019leveraging}. However, removing one linear band will give a better figure of merit $ZT$ for every type of heavy band.
\section{Model and methods}
\label{sec:th}
We assume the electronic energy dispersion $E(\mathbf{k})$ for the linear/Dirac and nonlinear/heavy bands by a common formula:
\begin{equation}
    E^{(n)} = A_n|\mathbf{k}|^n+\Delta_n,
    \label{eq:disp}
\end{equation}
where $n$ is the order of dispersion, $\mathbf{k}$ is the electronic wavevector, $\Delta$ is the energy shift from Fermi level, and $A_n$ is the normalization constant. The order $n=1$ and $n=2$ correspond to the Dirac and parabolic bands, respectively, while $n=4$ and $n=6$ both correspond to the pudding-mold band. The three-band system with parabolic, pudding-mold 4th order, and pudding-mold 6th order are depicted in Figs.~\ref{bandDOS3}(a), \ref{bandDOS3}(c), and \ref{bandDOS3}(e) respectively with the corresponding DOS in Figs.~\ref{bandDOS3}(b), \ref{bandDOS3}(d), and \ref{bandDOS3}(f). By removing the linear conduction band, the two-band system is modeled as shown in Figs.~\ref{bandDOS2}(a), \ref{bandDOS2}(c), and \ref{bandDOS2}(e) with the corresponding DOS in Figs.~\ref{bandDOS2}(b), \ref{bandDOS2}(d), and \ref{bandDOS2}(f).

\begin{figure}[tb]
    \centering
    \includegraphics[width=8cm]{FigTE-3bandsDOS.pdf}
    \caption{Schematics of the energy dispersion relations and density of states for three-band systems that consist of a parabolic band in (a) and (b), 4th-order pudding-mold band in (c) and (d), and 6th-order pudding-mold band in (e) and (f).}
    \label{bandDOS3}
\end{figure}

\begin{figure}[h]
    \centering
    \includegraphics[width=8cm]{FigTE-2bandsDOS.pdf}
    \caption{The energy dispersion and density of states for two-band system hat consist parabolic band in (a) and (b), pudding-mold 4th order band in (c) and (d), and pudding mold 6th order band in (e) and (f).}
    \label{bandDOS2}
\end{figure}

We use the Boltzmann transport theory within linear response regime and apply relaxation time approximation (RTA). Using this approximation, thermoelectric properties can be calculated from transport coefficients with the following TE kernel integration:
\begin{equation}
    \mathcal{L}_i = \int_{-\infty}^{\infty} \mathcal{T}(E)(E-\mu)^i\left(-\frac{\partial f}{\partial E}\right) dE,
\end{equation}
where $\mathcal{T}(E)$ is the transport distribution function (TDF), $\mu$ is the chemical potential,$f(E)$ is the Fermi-Dirac distribution function.
The TDF can be expressed with the following formula:
\begin{equation}
    \mathcal{T}(E) = v^2(E)\tau(E)\mathcal{D}(E),
\end{equation}
where $v$ is the velocity, $\tau$(E) is the relaxation time, and $\mathcal{D}(E)$ is the density of states (DOS).

Under RTA, the relaxation time is inversely proportional to the DOS
\begin{equation}\label{taudos}
    \tau(E)=\frac{\mathcal{C}}{\mathcal{D}(E)}
\end{equation}
where $\mathcal{C}$ is scattering coefficient. For 2D material, the longitudinal velocity $v(E)$ is related to the group velocity as written below
\begin{equation}
    v^2(E) = \frac{v^2_g}{2},
\end{equation}
where $v_g = \frac{1}{\hbar}\frac{dE}{dk}$ and $d$ is the dimension of the material. Then we get,

\begin{equation}\label{velocity}
    v^2(E) = \frac{n^2A^2_n}{2\hbar^2} \left(\frac{E-\Delta}{A_n}\right)^{\frac{2(n-1)}{n}}.
\end{equation}

For each electronic band, the transport coefficient will be calculated within proper integration boundaries. By using Boltzmann transport theory within relaxation time approximation the transport coefficients of the thermoelectric properties can be calculated ~\cite{hasdeo2019optimal}
\begin{equation}
    \sigma = q^2 \mathcal{L}_0,
\end{equation}
\begin{equation}\label{Seebeck}
    S = \frac{1}{qT} \frac{\mathcal{L}_1}{\mathcal{L}_0},
\end{equation}
and
\begin{equation}
    \kappa_e = \frac{1}{T} \left(\mathcal{L}_2 - \frac{(\mathcal{L}_1)^2}{\mathcal{L}_0}\right).
\end{equation}
In the case of the multiband system, the TE properties can be calculated by summing up the contribution from all the bands involved:
\begin{equation}
    \sigma = q^2 \sum_{j=1}^{N_B}\mathcal{L}_{0j},
\end{equation}
\begin{equation}
    S = \frac{1}{qT} \frac{\sum_{j=1}^{N_B}\mathcal{L}_{1j}}{\sum_{j=1}^n\mathcal{L}_{0j}},
\end{equation}
and
\begin{equation}
    \kappa_e = \frac{1}{T} \left(\sum_{j=1}^{N_B}\mathcal{L}_{2j} - \frac{(\sum_{j=1}^{N_B}\mathcal{L}_{1j})^2}{\sum_{j=1}^{N_B}\mathcal{L}_{0j}}\right)
\end{equation}
where $i$ is the index of the band and $n$ is the number of bands which are considered.
From those parameters, the figure of merit can be calculated as
\begin{equation}
    ZT = \frac{S^2\sigma}{\kappa_e+\kappa_{ph}} T,
\end{equation}
where $\kappa_{ph}$ is the lattice thermal conductivity and the power factor, PF
\begin{equation}
    PF = S^2\sigma.
\end{equation}

The TDF of each band using approximation in Eqs.~\eqref{taudos} and~\eqref{velocity} can be written as
\begin{equation}
    \mathcal{T}^{(n)}(E) = \mathcal{C} \frac{n^2A^2_n}{2\hbar^2} \left(\frac{E-\Delta}{A_n}\right)^{\frac{2(n-1)}{n}}
\end{equation}
By using dimensionless variables $\varepsilon = E/k_BT$, $\eta = \mu/k_bT$,  and $\tilde{\Delta} = \Delta/k_BT$ then the TE integral for Dirac conduction band
\begin{equation}
    \mathcal{L}^{(1)}_{i,c}=\frac{A^2_1\mathcal{C}}{2\hbar^2}(k_BT)^i \mathcal{F}_{i,c}(\eta-\tilde{\Delta})
\end{equation}
with
\begin{equation}
    \mathcal{F}_{i,c}(\eta) = \int_{-\eta}^\infty \frac{x^i}{(e^x+1)^2} dx.
\end{equation}

For the Dirac valence band we get
\begin{equation}
        \mathcal{L}^{(1)}_{i,v}=\frac{A^2_1\mathcal{C}}{2\hbar^2}(k_BT)^i \mathcal{F}_{i,v}(\eta+\tilde{\Delta})
\end{equation}
with
\begin{equation}
    \mathcal{F}_{i,v}(\eta) = \int_{-\infty}^{-\eta} \frac{x^i}{(e^x+1)^2} dx.
\end{equation}

For the parabolic, pudding mold 4th and 6th order, the TE integral can be written as
\begin{equation}
    \begin{split}
    \mathcal{L}^{(2)}_{i}=&\frac{4A_2\mathcal{C}}{2\hbar^2}(k_BT)^{i+1} \left[\mathcal{F}_{i+1}(\eta-\tilde{\Delta})+ (\eta-\tilde{\Delta})\right. \\ &\left. \times\mathcal{F}_i(\eta-\tilde{\Delta})-\mathcal{F}_{i+1}(\eta-\tilde{E}_0) \right.\\
    &\left.-(\eta-\tilde{\Delta})\mathcal{F}_i(\eta-\tilde{E}_0)\right],
\end{split}
\end{equation}

\begin{align}
    \mathcal{L}^{(4)}_{i}=\frac{16(A_4)^{\frac{1}{2}}\mathcal{C}}{2\hbar^2}(k_BT)^{i+\frac{3}{2}}\mathcal{H}_{4,i}(\eta,\tilde{\Delta},\tilde{E}_0)
\end{align}
with
\begin{align}
    \mathcal{H}_{4,i}(\eta,\tilde{\Delta},\tilde{E}_0) = \int_{\tilde{\Delta}}^{\tilde{E}_0} (\varepsilon-\tilde{\Delta})^{\frac{3}{2}}(\varepsilon-\eta)^i\frac{e^{\varepsilon-\eta}}{(e^{\varepsilon-\eta}+1)^2} d\varepsilon,
\end{align}
\begin{align}
    \mathcal{L}^{(6)}_{i}=\frac{36(A_4)^{\frac{1}{3}}\mathcal{C}}{2\hbar^2}(k_BT)^{i+\frac{5}{3}} \mathcal{H}_{6,i}(\eta,\tilde{\Delta},\tilde{E}_0)
\end{align}
with
\begin{align}
    \mathcal{H}_{6,i}(\eta,\tilde{\Delta},\tilde{E}_0) = \int_{\tilde{\Delta}}^{\tilde{E}_0} (\varepsilon-\tilde{\Delta})^{\frac{5}{3}}(\varepsilon-\eta)^i\frac{e^{\varepsilon-\eta}}{(e^{\varepsilon-\eta}+1)^2} d\varepsilon,
\end{align}
respectively.

The units of thermoelectric properties $S_0$, $\sigma_0$, and $\kappa_0$ are given by
\begin{equation}
    S_0 = \frac{k_B}{e} = 86.17\, \mu \text{V/K},
\end{equation}
\begin{equation}
    \sigma_0 = \frac{e^2\mathcal{C}}{2 \hbar},
\end{equation}
and
\begin{equation}
    \kappa_0 = \frac{1}{2 \hbar T}.
\end{equation}

\section{Results and Discussion}
\label{sec:res}

\subsection{Three-band systems}
\label{subsec:3band}
We calculate the thermoelectric properties of band structures which consist of three bands. Here we use $T = 300$K which corresponds to $k_BT = 0.026$ eV. The gap between Dirac bands is set to $2k_BT$ while the bandwidth of the nonlinear band is $5$ eV. We assume the ideal system which the lattice thermal conductivity $\kappa_{ph} =0$.

First, we calculate the thermoelectric properties of three-band model with these set of constants $A_1^c = 2.5\, \text{eVm}$ , $A_1^v = -2.5\, \text{eVm}$, $A_2 = 2.5\, \text{eVm}^2$, $A_4 = 2.5\, \text{eVm}^4$ and $A_6 = 2.5 \, \text{eVm}^6$. We plot the band structure model for a three-band system with parabolic inside in Fig.~\ref{3combined-symmetry}(a) and its thermoelectric properties such as Seebeck coefficient, electrical conductivity, and electronic thermal conductivity in Figs.~\ref{3combined-symmetry}(b-d). We plot the thermoelectric properties as a function of a dimensionless variable from its chemical potential $\mu$. For each thermoelectric property, we plot each type of three-band system with a different heavy band. In addition to the three-band system, we also plot the thermoelectric properties of the pure Dirac material system.

In Fig.~\ref{3combined-symmetry}(b), we have plotted the Seebeck coefficient for those three-band systems. First, we can see that the value of the Seebeck coefficient is similar to all of the systems in negative chemical potential as shown in that the curves coincide. Because in this negative region, only the Dirac band exists for all of the systems. At the Fermi level ($\mu/k_BT = 0$), it starts to have a different value for each heavy band system. The difference between each system becomes more evident in the positive chemical potential. The value of the Seebeck coefficient of the three-band systems is just slightly larger than the Dirac cones system. It is due to the existence of a heavy band in the system that generates energy filtering for impeding the low carrier to get excited. This process will imbalance the number of careers available in the band. As we know that the more imbalance carrier between the band, the higher the Seebeck coefficient will be.
\begin{figure}[ht!]
    \centering
    \includegraphics[width=8cm]{FigTE-3bands-combined-symmetry.pdf}
    \caption{Thermoelectric properties of a three-band system with symmetric Dirac bands. The band structure in (a), the Seebeck coefficient ($S$) in (b), electrical conductivity in (c), and electronic thermal conductivity in (d) for each combination. The results for $S,\sigma$ and $\kappa_e$ are expressed in the units of $S_0$, $\sigma_0$, and $\kappa_0$ respectively.}
    \label{3combined-symmetry}
\end{figure}

In the case of electrical conductivity in Fig.~\ref{3combined-symmetry}(c), the electrical conductivity in the vicinity of the Fermi level for the three-band system is slightly higher than the Dirac material. It becomes much different in the positive chemical potential region. For the Dirac material, the electrical conductivity is symmetric with respect to the Fermi level (between p-type and n-type) while the three-band system is not symmetric anymore. The highest electrical conductivity is achieved by the pudding-mold 6th order system in positive chemical potential. This large value comes from the dependence of $n$ in the thermoelectric integral $\mathcal{L}_0$ wherein the Seebeck coefficient cancels out because of the division of $\mathcal{L}_1$ and $\mathcal{L}_0$ in Eq.~\eqref{Seebeck}. On the other hand, the more evident difference among the nonlinear band appears also in the electronic thermal conductivity. The electronic thermal conductivity is somewhat more related to the electrical conductivity than the other TE properties. Near the Fermi level, the electrical conductivity of parabolic and pudding-mold 6th order is similar which is also higher than the pudding-mold 4th order. In the positive chemical potential, the electrical conductivity of the pudding-mold 6th order system becomes larger. This suggests that the three-band system will provide better TE performance on n-doping than p-doping.

In general, the Seebeck coefficient will be higher if the career concentration is small enough and there is a difference between the conduction band career and valence band career that are electron and hole respectively. The Seebeck coefficient for each nonlinear three-band system is rather similar. We see many differences between this nonlinear band system in the electrical conductivity and electronic thermal conductivity. From these three TE properties, we plot the power factor and the figure of merit in Figs.~\ref{PFZT3-symmetry}. We did some changes to the position of the nonlinear band. We shifted up $4k_BT$ and $8k_BT$ in Figs.\ref{PFZT3-symmetry}(b and e) and (c and f) respectively. We could see a higher PF when the nonlinear band intersects the Dirac band at $\Delta = 1k_BT$. As the nonlinear band shifted up, the value of PF of the three-band system will reduce to the PF of Dirac material. We could see that the combination with pudding-mold 6th order features the highest PF among the others in the positive region of chemical potential.

However, from Figs.~\ref{PFZT3-symmetry}(c-f) we can see that the figure of merit of the three-band system is much similar to the Dirac material. The most evident difference can be seen when the nonlinear band intersects the same point with the Dirac band compared to others ($\Delta = 5k_BT$ and $\Delta = 9k_BT$). Despite higher PF for the three-band system where the nonlinear band intersects the Dirac band, the large electronic thermal conductivity $\kappa_e$ of this system cancels the promising value of PF as the figure of merit is inversely proportional to $\kappa_e$. We also noted that the further we shift up the nonlinear band, it will reduce to a two-band system as if the nonlinear band does not exist and only Dirac cones exist in the system.

\begin{figure*}[ht!]
    \centering
    \includegraphics[width=17cm]{PFZT-3bands.pdf}
    \caption{Power factor (PF) in (a-c) and electronic figure of merit ($ZT_e$) in (d-f) within three-band model. The position of the nonlinear band is varied corresponding to the Fermi level. The nonlinear band is shifted up $1k_BT$ in (a) and (d), $5k_BT$ in (b) and (e), and $9k_BT$ in (c) and (f).}
    \label{PFZT3-symmetry}
\end{figure*}

One example of real material that represents the most similar to this three-band model was found on the square-octagonal \ch{MoS2}. From the band structure in Fig. \ref{so-MoS2BandDOS}, we can see that there are Dirac cone and a parabolic band near the Fermi level. By doing numerical fitting to this DFT band structure using the dispersion relation in Eq.~\eqref{eq:disp}, we found $A_1^c = 2.472 \, \text{eVm}$, $A_1^v = -2.556 \, \text{eVm}$, and $A_2 = 2.563 \, \text{eVm}^2$. Unfortunately, the thermoelectric properties in Fig.~\ref{3combined-symmetry} are not in the same trend as shown from the first principle result in Fig.~\ref{TEBoltzTrap2} though we have used a similar value of the coefficient $A_n$ for each respective band.

From the previous result, we do some changes deliberately to these sets of coefficients to get a similar result to the first principle calculation. We found these set of coefficients $A_1^c = 1.4 \, \text{eVm}$, $A_1^v = -3.3 \, \text{eVm}$, and $A_2 = 2.5 \, \text{eVm}^2$ give the similar result to the first principle result as shown in Fig~\ref{3combined-assymmetry}. Compared to the numerical fitting result, this set of coefficients is much different from the linear band coefficient. This difference can be understood as an effective model to capture the effect of the other bands. In the first principle calculation, there is more than three bands are considered to generate the thermoelectric properties of the material.

\begin{figure}[ht!]
    \centering
    \includegraphics[width=8cm]{FigTE-3bands-combined-asymmetry.pdf}
    \caption{Thermoelectric properties of a three-band system with asymmetric Dirac bands. The band structure in (a), the Seebeck coefficient ($S$) in (b), electrical conductivity in (c), and electronic thermal conductivity in (d) for each combination. The results for $S,\sigma$ and $\kappa_e$ are expressed in the units of $S_0$, $\sigma_0$, and $\kappa_0$ respectively.}
    \label{3combined-assymmetry}
\end{figure}

\subsection{Two-band systems}
\label{subsec:2band}
Now, we remove the linear conduction band to make a two-band system. The two-band system is similar to the Dirac material system unless the conduction linear band is substituted by a nonlinear band as we can see in Fig.~\ref{2combined}(a). Doing the same calculation for the thermoelectric properties, we can see from Fig.~\ref{2combined}(b) that the Seebeck coefficient is generally higher than the three-band system or even the pristine Dirac material system. In this two-band system, we can see from Fig. \ref{2combined}(b) that the difference for each combination is slightly more evident than in the three-band systems. The two-band model that consists of pudding-mold 4th order features a higher Seebeck coefficient though it is just a tiny difference. This large Seebeck coefficient comes from the asymmetry condition in this two-band system. In the case of a pristine Dirac material system, we plot the same Dirac or linear band in both the valence and conduction band. This will generate a small Seebeck coefficient.

Moreover, we can see from Fig. \ref{2combined}(c) that the electrical conductivity is lower than in the three-band system for all combinations. Lower electrical conductivity will also compensate for the higher Seebeck coefficient that we got previously. Consequently, the power factor will not generally be higher than in the three-band system. Figure \ref{2combined}(f) depicts the electron thermal conductivity for the two-band system. The electron thermal conductivity is lower than in the three-band systems. It will result in delivering a higher power factor and $ZT_e$.

Now we shift up the nonlinear band from the original position $1k_BT$ to some higher values such as $5k_BT$ and $9k_BT$. From Figs.~\ref{PFZT2}(a-c), the power factor will increase as the nonlinear band is shifted up. However, the power factor is similar between $\Delta 5k_BT$ and $\Delta 9k_BT$. The large difference between the two-band and the three-band system comes in the value of the electronic figure of merit $ZT_e$ in Figs.~\ref{PFZT2}(d-f). The more we shifted up the nonlinear band, means higher the bandgap, the higher $ZT_e$ will be. The smallest $ZT_e$ when $\Delta = 1k_BT$ for a two-band system is 10 times larger than the same $\delta$ for a three-band system. For the larger bandgap in Figs.\ref{PFZT2}(e and f), the $ZT_e$ becomes extremely larger than the previous one. Another fact that we found in this two-band system is the clear difference between parabolic, pudding-mold 4th, and pudding-mold 6th band in the $ZT_e$. We can see from Figs.~\ref{PFZT2}(d-f) that the two-band system with a pudding-mold band is higher than the parabolic system as it is presumed from the theory that the pudding-mold band has a drastic rise in the DOS as shown in Fig.~\ref{bandDOS2}(d and f).

\begin{figure}[ht!]
    \centering
    \includegraphics[width=8cm]{FigTE-2bands-combined.pdf}
    \caption{Thermoelectric properties of the two-band system. Seebeck coefficient ($S$) in (a), electrical conductivity in (b), and electronic thermal conductivity in (c) for each combination. The results for $S,\sigma$ and $\kappa_e$ are expressed in the units of $S_0$, $\sigma_0$, and $\kappa_0$ respectively.}
    \label{2combined}
\end{figure}

\begin{figure*}[ht!]
    \centering
    \includegraphics[width=17cm]{PFZT-2bands.pdf}
    \caption{Power factor (PF) in (a-c) and electronic figure of merit ($ZT_e$) in (d-f) within two-band model. The position of the nonlinear band is varied corresponding to the Fermi level. The nonlinear band is shifted up $1k_BT$ in (a) and (d), $5k_BT$ in (b) and (e), and $9k_BT$ in (c) and (f).}
    \label{PFZT2}
\end{figure*}

\subsection{Effect of bandwidth}
Here, we will present the effect of the bandwidth of the nonlinear band on the thermoelectric properties. A previous study point out that there was no optimum bandwidth for $\tau = \frac{\mathcal{C}}{DOS}$ scattering approximation\cite{PhysRevLett.107.226601}. In a three-band system, we found that there exists an optimum value of the bandwidth to bring the optimum thermoelectric properties. The Seebeck coefficient in Fig. ~\ref{3band-bandwidth} (a,f,k), there is no large difference for each bandwidth in every type of nonlinear band. In contrast with the Seebeck coefficient, the electrical conductivity (Fig.~\ref{3band-bandwidth} (b,g,l) and electronic thermal conductivity (Fig.~\ref{3band-bandwidth} (c,h,m)) do have apparent difference for each bandwidth, specifically in the positive regime of the Fermi level. Less bandwidth utilized features less electrical conductivity and electronic thermal conductivity. This results in a higher power factor for the higher bandwidth. From Fig.~\ref{3band-bandwidth}(d, i,n), we can see that the largest PF is achieved when the bandwidth is set to $10k_BT$. Increasing the bandwidth above $10k_BT$ such as $15k_BT$ does not give us a higher PF than the former one. Lastly, the figure of merit in Fig.~\ref{3band-bandwidth} (e,j,o) show a similar result for each bandwidth. The large power factor is compensated by the large electronic thermal conductivity. It results in a lower figure of merit (less than $0.1$).

We apply the same calculation for the two-band system. In this case, we found that the Seebeck coefficient for small bandwidth ($5k_BT$) in Fig.~\ref{2band-bandwidth}(a,f,k) is higher than the others at the chemical potential above the Fermi level. A similar trend to the three-band system, the electrical conductivity (Fig.~\ref{2band-bandwidth}(b,g,l)) and the electron thermal conductivity (Fig.`\ref{PFZT2}(c,h,m) is less in lower bandwidth. However, the power factor (Fig.~\ref{2band-bandwidth}) for each bandwidth does not have any large difference among them. The clear difference between the bandwidth happens in the figure of merit (Fig.~\ref{2band-bandwidth}(e,j,o)) where the small bandwidth features larger $ZT_e$ than the higher bandwidth.  It is physically reasonable because, in a two-band system, less bandwidth of the nonlinear band will make this band looks flatter.
\begin{figure*}[ht!]
    \centering
    \includegraphics[width=17cm]{FigTE-3bands-bandwidth.pdf}
    \caption{Seebeck coefficient, electrical conductivity, electronic thermal conductivity, power factor, and $ZT_e$ within three-band models with parabolic in (a-e), pudding-mold 4th in (f-j), and pudding mold 6th band in (k-o) varied by different bandwidth of the nonlinear band.}
    \label{3band-bandwidth}
\end{figure*}

\begin{figure*}[ht!]
    \centering
    \includegraphics[width=17cm]{FigTE-2bands-bandwidth.pdf}
    \caption{Seebeck coefficient, electrical conductivity, electronic thermal conductivity, power factor, and $ZT_e$ within two-band models with parabolic in (a-e), pudding-mold 4th in (f-j), and pudding mold 6th band in (k-o) varied by different bandwidth of the nonlinear band.}
    \label{2band-bandwidth}
\end{figure*}

\section{Conclusions}
\label{sec:con}
We find that in the three-band systems, any nonlinear band accompanying the Dirac band will generate similar TE properties near the Fermi level. The difference will emerge mostly in the positive chemical potential of the pudding-mold band system, which features a higher power factor and electronic figure of merit. The optimum thermoelectric power factor is found if the nonlinear band intersects the Dirac conduction band. The three-band system will reduce to a two-band Dirac material system when the nonlinear band is more than $k_BT$ above the Fermi level. The two-band system exhibits better TE performance in general than the three-band system. Compared to the three-band system, the difference given by the various nonlinear bands becomes more obvious in the two-band system, where the pudding-mold system offers the higher electronic figure of merit ($ZT_E$). We can conclude that fewer bands will give better TE performance than multiband systems but a single-band thermoelectrics in reality is not feasible. Therefore, materials that exhibit two bands near the Fermi level will give the best thermoelectric performance.

\begin{acknowledgements}
A.D. performed first principles calculations using King Abdulaziz University (KAU)’s High Performance Computing Center (Aziz Supercomputer). He is also supported by the KAU scholarship. E.S. and A.R.T.N acknowledge computational facilities from BRIN HPC.
\end{acknowledgements}