# Network Graph Based Neural Architecture Search

Zhenhan Huang1 Chunheng Jiang1 Pin-Yu Chen2&Jianxi Gao1111Contact Author
1Department of Computer Science, Rensselaer Polytechnic Institute, Troy, NY
12180
2IBM Thomas J. Watson Research Center, Yorktown Heights, NY 10598
{huangz12<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Neural architecture search enables automation of architecture design. Despite
its success, it is computationally costly and does not provide insight on how
to design a desirable architecture. Here we propose a new way of searching
neural network where we search neural architecture by rewiring the
corresponding graph and predict the architecture performance by graph
properties. Because we do not perform machine learning over the entire graph
space and use predicted architecture performance to search architecture, the
searching process is remarkably efficient. We find graph based search can give
a reasonably good prediction of desirable architecture. In addition, we find
graph properties that are effective to predict architecture performance. Our
work proposes a new way of searching neural architecture and provides insights
on neural architecture design.

##  1 Introduction

Neural networks architecture achieves a great success over last few years in
various challenging applications, such as image classification Krizhevsky et
al. (2012), speech recognition Hinton et al. (2012) and machine translation Wu
et al. (2016). The success, largely attributed to the feature engineering, is
accompanied by the arise of architecture engineering, e.g., AlexNet Krizhevsky
et al. (2012), VGGNet Simonyan and Zisserman (2014), ResNet He et al. (2016).
Increasingly complicated neural network architectures make it progressively
more difficult for manual design of architecture. Neural architecture search
(NAS), enabling automation of architecture engineering, proves its potential
in boosting the performance of neural network architecture. For example, NAS
methods have outperformed some manually designed architecture in some
applications, such as image classification Real et al. (2019), object
detection Zoph et al. (2018).

In NAS method, there are three major components Elsken et al. (2019): search
space $\mathcal{S}$, search strategy ${\mathcal{A}}_{s}$ and performance
estimation strategy ${\mathcal{A}}_{e}$. Predefined $\mathcal{S}$ confines the
total number of possible architectures if it is not unbounded. Thus it will
affect the search efficiency and optimal architecture. ${\mathcal{A}}_{s}$
determines the search efficiency and should avoid potential pitfall of local
minimum. ${\mathcal{A}}_{e}$ provides a way to evaluate architecture candidate
and feedback. The simplest way for ${\mathcal{A}}_{e}$ is to perform a
standard training and validation for targeting optimal architecture.

Powerful as NAS method is, searching an architecture in $\mathcal{S}$ and
testing its performance will take ample time. Although decomposing hand-
crafted architectures into motifs and searching motifs, blocks or cells can
boost searching speed, this process inevitably introduces bias in the search
space. Furthermore, NAS method does not throw light on why specific
architectures outperform others and general principles of designing an
architecture. Network graphs, on the other hand, have good metrics for
evaluation. For example, clustering coefficient measures the degree to which
nodes a in graph cluster. If we can bridge graph and neural network, we will
be able to relate graph properties to architecture performance. Instead of
searching in architecture space, we can perform search in graph space. By
targeting the optimal graph structure, we will able to locate the optimal
neural network architecture.

Figure 1 shows the graph-based neural architecture search. The relationship
between graph and neural network is bridged by relational graph. We predict
machine learning error based on graph properties. A set of graph features is
selected for the prediction. Based on predicted error, we apply rewiring
strategy to search a better or worse neural architecture.

Figure 1: Schematic illustration of graph-based neural architecture search

##  2 Bridge of Neural Network and Graph

###  2.1 Graph Representation of Neural Network

The relation between graph and neural network can be bridged by relational
graph You et al. (2020). For a fixed-width multilayer perceptron (MLP) where
each layer contains the same number $R$ of computation units (neurons),
suppose the input and output of the $r$-th layer ($1\leq r\leq R$) is
$\mathbf{X}^{(r)}$ and $\mathbf{X}^{(r+1)}$, respectively. A neuron in the
$r$-th layer computes You et al. (2020):

$x^{(r+1)}_{i}=\sigma(\sum_{j\in N(i)}\omega^{(r)}_{ij}x^{(r)}_{j}),$ (1)

where $w^{(r)}_{ij}$ is the $i$-th row and $j$-th column of the trainable
weight $\mathbf{W}^{(r)}$, $x^{(r)}_{j}$ is the $j$-th dimension of the input
$\mathbf{X}_{j}$, $x^{(r)}_{i}$ is the $i$-th dimension of the output
$\mathbf{X}_{i}$. $\sigma(\cdot)$ is the activation function introducing non-
linearity. $N(i)$ is the set of neurons of the $(r+1)$-th layer connecting to
$r$-th layer and defined by relational graph $G$.

###  2.2 Undirected Graph Generator

Graphs can be categorized into two groups: directed graph or undirected graph.
We consider relational graph as undirected, so its adjacency matrix is
symmetric. The classic graph generation algorithms include (1) Watts-Strogatz
(WS) model that can generate graphs with small-world properties Watts and
Strogatz (1998). The ”small world effect” is generally referred to describe
graphs whose average path length is comparable with a homogeneous random
graphs Prettejohn et al. (2011); (2) Erdős-Rényi (ER) model that can generate
random graphs Erdos et al. (1960); (3) Barabási-Albert (BA) model that
constructs scale-free graphs Barabási and Albert (1999). A scale-free graph
has a degree distribution following power law, i.e., the probability of a node
having a degree $k$ has a scale-invariant decay $P(k)\sim k^{-\gamma}$, where
$\gamma$ is a constant and confined by $\gamma>1$; (4) Harary model that
generates graphs with maximum connectivity Harary (1962). WS-flex graph
generator, proposed in You et al. (2020), is a more general form of WS model
by relaxing the constraint that all the nodes have the same degree before
random rewiring. After fixing the number of nodes to be 64, WS-flex graph
generators are capable of generating graphs encompassing almost all graphs
constructed by classic graph generation algorithms mentioned above You et al.
(2020).

##  3 Experimental Setup

We use the CIFAR-10 dataset Krizhevsky et al. (2009) for training MLPs.
CIFAR-10 dataset has 50K training images and 10K validation images. Each image
is a matrix with a dimension of $32\times 32\times 3$. 60K images are
classified into 10 categories for supervised machine learning.

###  3.1 Graph Space

We use WS-flex graph generator to generate undirected graphs. The number of
nodes is fixed to be $64$ and the average degree range is $[2,63]$. Total
number of generated graphs is $5983$. We calculate $26$ Graph properties
including average degree, clustering coefficient, heterogeneity, average path
length, bimodularity, greedy modularity, resilience parameter, degree entropy,
wedge count, gini index, average node connectivity, edge connectivity, average
closeness centrality, average closeness centrality (WF improvement), average
eccentricity, diameter, radius, average edge betweenness centrality, average
node betweenness centrality, central point of dominance, core number, minimum
Laplacian spectrum, maximum Laplacian spectrum, transitivity, local
efficiency, global efficiency. Detailed description regarding to those
properties can be referred in Appendix A.

We find that heterogeneity range for $5983$ graphs constructed by ws-flex
graph generator is $[0,0.82]$. Most graphs are not heterogeneous in the
structure. To introduce more heterogeneous graphs, we do random rewiring based
on those graphs. Total number of graphs becomes $19724$.

###  3.2 Feature Selection

Average path length and clustering coefficient are considered as good
indicator for predicting learning error of MLP You et al. (2020). We calculate
28 features for each graph and take into consideration the potential of these
features for predicting machine learning performance. To filter features, we
use sequential forward selection (SFS) algorithm Marcano-Cedeño et al. (2010);
Ververidis and Kotropoulos (2005); Cotter et al. (1999), a commonly used
method for reducing the data dimension. SFS algorithm is a bottom-up search
strategy in which an empty set $S$, the starting point, continuously add
features until all features are added. At each iteration, $S$ greedily
searches for best feature from remaining feature pool and include it.

Selecting MLP learning error as output variable and graph properties as input
variables, We use linear regression to fit results and choose mean squared
error (MSE) as metric to select features. We use random splitting strategy to
obtain training set and test set. The ratio of training set to test set is
$9:1$. At each iteration of SFS algorithm, we use training set to determine
linear regression parameters and test set to calculate evaluation parameters
MSE and Pearson correlation coefficient.

###  3.3 Neural Network Architecture

We use a 5-layer MLPs as the neural network architecture. Each MLP layer has
512 hidden units as baseline architecture. The relational graph determines the
connection of hidden units between two neighboring layer. To ensure all
networks have the approximately same complexity, we use FLOPS (# of
multiplication and addition) as metric to adjust baseline architecture such
that FLOPS for different relational graph represented networks is roughly the
same.

Fig. 1 a) shows the illustration of MLP architecture. Each MLP layer contains
Batchnorm ReLU layer to introduce non-linearity and BatchNorm layer Ioffe and
Szegedy (2015). Batch size in the training process is 128 and total number of
epochs for training a model is 200. We use a decaying learning rate with a
cosine annealing Loshchilov and Hutter (2016). Our Momentum Optimizer has an
initial learning rate of 0.1, 5e-4 weight decay, 0.9 momentum and uses
Nesterov Momentum Sutskever et al. (2013); Nesterov (1983).

###  3.4 Rewiring Algorithm

To search the best neural architecture, we adopt a greedy strategy which is
different from the conventional approaches in NAS. We start from a neural
architecture with known performance, improving the architecture by iteratively
rewiring it. Let $G_{t}$ be the current architecture, $P(G_{t})$ be the post-
training performance of $G_{t}$, and $\mathcal{A}$ be our rewiring strategy,
we can formulate the rewiring procedure as

$G_{t+1}=\operatorname*{arg\,max}_{\mathcal{A}}P(\mathcal{A}(G_{t})).$ (2)

Nevertheless, training each $\mathcal{A}(G_{t})$ for the performance
$P(\mathcal{A}(G_{t}))$ will make the search inevitably expensive. To address
this issue, we refer to a surrogate parametric model $F(G;\bm{\theta})$ for
$P(G)$, which is designed to capture the relationship between a set of simple
topological properties (see Appendix A) and the performance of $G_{t}$, i.e.
$\hat{P}(G)=F(G;\bm{\theta})$. To some extend, the “goodness” of $F$
determines our choices of the resultant architectures via rewiring. We search
the architecture by selecting a good rewired candidate or discarding a bad
one. If $F(G_{t};\bm{\theta})$ is consistent with the real performance
$P(G_{t})$, there will be very few mistakes in our choices and our search will
guarantee a constantly improved architecture. Also, some architectures may be
a local optimal, which can not be rewired to reach a better one.

As shown in Algorithm 1, our approach starts from a randomly selected relation
graph whose associated neural network post-trained performance may be
arbitrarily low. We perform a greedy search of a sequence of rewiring
operations to constantly improve the performance of the associated neural
network of the initial relational graph. Currently, our rewiring strategy
allows removal of edges $\mathcal{A}_{\rm rmv}$, building new edges
$\mathcal{A}_{\rm new}$, random rewiring $\mathcal{A}_{\rm rnd}$ and double
edge swaps $\mathcal{A}_{\rm swap}$.

Algorithm 1 Rewiring and Searching

Input: A randomly selected relation graph $G_{0}\in\mathcal{G}$, four rewiring
operations $\\{\mathcal{A}_{\rm new},\mathcal{A}_{\rm rmv},\mathcal{A}_{\rm
swap},\mathcal{A}_{\rm rnd}\\}$, accepted relative improvement
$\varepsilon\in(0,1)$, maximum number $K$ of rewiring operations, the
performance predictor $F(G;\bm{\theta})$ of the associated neural network
$\phi(G)$ with the topological properties of a relational graph $G$
Output: Rewired relation graph $G_{K}$

1: Let $k=0$, $\hat{P}_{0}=F(G_{0};\bm{\theta})$.

2: while $k<K$ do

3: Randomly select a rewiring operation $\mathcal{A}$

4: Perform $\mathcal{A}$ over $G_{k-1}$ and obtain
$\hat{G}_{k-1}=\mathcal{A}(G_{k-1})$

5: Predict the post-trained performance of $\phi(\hat{G}_{k-1})$ with
$\hat{P}_{k}=F(\hat{G}_{k-1};\bm{\theta})$

6: if $|\hat{P}_{k}/\hat{P}_{k-1}-1|\geq\varepsilon$ then

7: Accept $\mathcal{A}$ with $G_{k}=\hat{G}_{k-1}$ and set $k=k+1$

8: else

9: Reject $\mathcal{A}$ and repeat the above procedure

10: end if

11: end while

##  4 Experimental Result

###  4.1 Features for prediction

Calculating all features of a graph can be computational expensive. Hence, we
are interested in how many features are needed to predict MLP performance and
what are important features. We use SFS algorithm to quantitatively show the
relation between number of features and prediction quality, as well as to
determine significant features.

Figure 2 shows MSE and Pearson correlation coefficient variation with
including features. Instead of typical oscillating trend of SFS algorithm, MSE
shows a monotone decreasing with adding features. Even though more features
give a better prediction quality, the gain in prediction boost becomes
negligible while the computational cost increases dramatically when the number
of features is more than $10$. At the same time, we find linear regression can
give a good prediction. If we only choose one feature, MSE is $0.086$. Adding
just $4$ more features, MSE drops to $55\%$ that of one feature. Considering
the trade-off between prediction performance and computation cost, we choose
the first 10 features to predict MLP learning error. The corresponding MSE is
$0.4$ and we have a high $0.93$ Pearson correlation coefficient. The finding
suggests that a small number of features can still give a good prediction of
MLP performance.

Figure 2: Sequential forward selection for feature selection

SFS algorithm determines the eccentricity as the starting point. We are
interested, within a certain tolerance on prediction error of MLP performance,
can we choose other features as starting point that give a similar prediction
quality? There are several benefits for moderately sacrificing prediction
quality: 1) the computational cost for different features can vary a lot. If
we replace features having a high computational cost with ones having low cost
and, at the same time, the prediction quality does not change much, the MLP
prediction model will be satisfying; 2) there are highly developed theories
about some of features. If we use those features for prediction, it is
beneficial to understand why and how a neural network has a good performance.
For example, the eccentricity measures the maximum distance from one node to
all other node. The average path length measures the average of minimum
distance between all pairs of nodes. These two features are closely related to
message exchange efficiency. If we can replace the eccentricity with the
average path length, we will be able to analyze neural network performance
using the average path length related theories, e.g., small world network
theory; 3) in reality, it is unlikely a certain rewiring strategy just alters
one graph property while other property remains unchanged. If we can replace
less controllable feature with more controllable one, it becomes easier to
target optimal relational graph.

We replace the first feature, i.e., eccentricity, with other features and
perform SFS algorithm to determine remaining features. Figure 3 shows
representative MSE variation as a function of number of features in SFS
algorithm after fixing the first feature. Generally, a different starting
point of SFS algorithm will lead to different feature set $S$ and different
MSE. However, some features give a similar prediction quality and feature set.
For example, MSE for average path length and eccentricity at each iteration is
very close. Besides, the order of adding features is nearly the same. When we
consider the first $10$ features excluding first feature, the eccentricity and
the average path length have the same remaining $9$ features. The feature
adding order is also the same.

When a single feature is used for predicting MLP performance, calculated
features have a pronounced difference in the prediction quality. Some
features, such as heterogeneity, have a high $0.3$ MSE. Nevertheless, with
adding more features, all features show a fast converge (adding less than 10
features) to a small $0.4$ MSE. This indicates that a linear combination of
features gives a good prediction of MLP performance. More details about MSE
variation with adding features when fixing the first feature can be referred
in Appendix C.

Figure 3: Representative sequential forward selection after fixing the first
feature

Table 1 shows first $10$ features determined by SFS algorithm and first $10$
features by SFS algorithm after fixing first feature to be the average path
length. Except the first feature, two feature sets have the exactly the same
features and linear regression parameters associated with these features have
the same sign, i.e., consistent positive correlation or negative correlation.
In addition, MSE is close for these two feature sets. Hence, we think the
average eccentricity and average path length are interchangeable in terms of
predicting MLP performance.

Algorithm |  Selected Features |  Sign
---|---|---
SFS |  Average eccentricity |  $+$
Central point of dominance |  $+$
Resilience Parameter |  $+$
Global efficiency |  $-$
Edge connectivity |  $-$
Wedge count |  $-$
Clustering coefficient |  $-$
Average node connectivity |  $+$
Average closeness centrality |  $+$
Greedy modularity |  $-$
|  Average path length |  $+$
|  Central point of dominance |  $+$
|  Resilience Parameter |  $+$
|  Global efficiency |  $-$
SFS fixing |  Edge connectivity |  $-$
first feature |  Wedge count |  $-$
|  Clustering coefficient |  $-$
|  Average node connectivity |  $+$
|  Average closeness centrality |  $+$
|  Greedy modularity |  $-$
Table 1: Feature sets determined by SFS algorithms Figure 4: Pearson
correlation coefficient of every pair of graph features. Feature sets related
to every feature are determined by SFS algorithm

Except for the pair of the average eccentricity and average path length, we
are also curious about a question – are there other features interchangeable.
We use Pearson correlation coefficient to characterize the similarity between
every pair of feature sets. We fix first feature and use SFS algorithm to
determine the remaining $9$ features. Every remaining $9$ features constitute
a set $S$. Then we calculate Pearson correlation coefficient between pairs of
those sets. Figure 4 shows Pearson correlation coefficient between all pairs
of features. The result reveals the similar effect in predicting MLP
performance. The similarity seems to be correlated to physical meaning of
graph features. For example, the average path length measures the average of
shortest path length between pairs of nodes Albert and Barabási (2002); Achard
and Bullmore (2007). The average eccentricity measures the maximum distance
from one node to other nodes in a graph $G$ Hage and Harary (1995). Diameter
is the maximum eccentricity of a graph $G$. Average node betweenness
centrality calculates the average of node betweenness centrality which is the
sum of the fraction of all shortest paths passing through a node Brandes
(2001, 2008). These four features are related to message exchange efficiency
in a graph. At the same time, they have a highest correlation 1.

###  4.2 Prediction of MLP performance

We choose as the starting point a graph with a highest MLP top 1 error and a
graph with lowest MLP top 1 error in the graph pool. Then we use the rewiring
strategy described before to rewire graph for searching a better graph
structure step by step. At each iteration, if the predicted error decreases by
more than a threshold value $0.01$, the rewired graph is accepted and the
change of graph structure due to rewiring is kept. To verify the prediction of
performance, we use MLPs represented by relational graphs at each stage to
calculate the real top 1 error. Total time for rewiring is $6$ hours.

Figure 5 shows the top 1 error decrease path (superior rewiring) and Figure 6
shows the top 1 error increase path (inferior rewiring). At the end of the
rewiring process, both superior rewiring and inferior rewiring show a dramatic
increase in the computational cost. In both rewirings, real top 1 error
oscillates around the predicted top 1 error. The prediction is promising as
the real top 1 error overall follows the prediction path. During the rewiring
process, we do not observe the sudden jump in the top 1 error. We think that
limited change in the graph edge set does not change MLP performance
dramatically. This is consistent with the result from conventaional NAS: small
perturbations in the neural network architecture does not have a significant
effect on its performance Ru et al. (2020). During the superior rewiring
process, the computational cost increases drastically even though the learning
error does not hit the minimal error region. This might be attributed to
premature convergence to the region of suboptimal structure. There seems to be
a exploration-exploitation trade-off problem. The greedy searching strategy
might lead to a suboptimal graph structure.

(a) Top 1 error decrease path

(b) Computational cost

Figure 5: Rewiring path of MLP learning error decrease

(a) Top 1 error increase path

(b) Computational cost

Figure 6: Rewiring path of MLP learning error decrease

We choose the graph in inferior rewiring as the starting point and repeat our
rewiring strategy for 100 times. We use 100 random seeds to generate different
random rewiring processes. Figure 7 shows the statistical result of rewiring
result. We select continuous step periods, e.g., step 1 to step 10 as one step
range and then step 11 to step 20 as another step range. Then we calculate the
average of real MLP top 1 errors of each step range. As the number of rewiring
steps increases, the median number of average top 1 error increases. This
trend is consistent with our prediction, i.e., the statistical increase in top
1 error has a same trend as the predicted increase in top 1 error. Hence, the
performance predictor $F(G;\bm{\theta})$ can be used for predicting MLP
performance.

Figure 7: Statistical result of rewiring path

##  5 Discussion

The performance of neural network architecture can be related to properties of
corresponding relational graph. For example, the average path length is a
measurement of efficiency of information transport over a graph. At the same
time, it is a good indicator of MLP performance. In general, a lower average
path length is preferred over a higher one (refer to A). However, there is no
continuous function fitting well a single graph property and MLP performance.
It is impractical to just rely on one graph feature to predict MLP
performance. A combination of graph properties, nevertheless, can give a
reasonably good prediction of the performance. By predicting a learning error,
we are able to search an optimal neural network architecture just by rewiring
a graph and calculating graph properties. The benefit of avoiding performing
machine learning over an entire graph space is high time efficiency and low
computational cost. We adopt linear regression to calculate predicted
performance based on graph properties, this process is computationally cheap.
Because the graph space is unbounded and does not rely on typical manual
designed architecture, the search space is unbiased. Overal, the real MLP
performance is consistent with the predicted MLP performance. Based on
predicted MLP performance, we are able to target inferior or superior graph
structure.

We find different properties might have a similar effect in machine learning
accuracy. In other words, some properties are interchangeable. The
interchangeability is consistent to physical meaning of those properties. For
instance, the average path length, diameter and eccentricity are all related
to message exchange efficiency. They are interchangeable in terms of
predicting machine learning performance. A directed rewiring strategy,
targeting at changing a specific group of graph properties, will definitely
enhance the search efficiency.

## Acknowledgments

This work is supported by Center for Computational Innovation (CCI) of
Rensselaer Polytechnic Institute (RPI).

##  Appendix A Features of a Graph Network

Average degree. The average degree $\bar{k}$ calculates the average number of
edges for one node Luce and Perry (1949). For a graph with $n$ nodes, the
average degree $\bar{k}$ is given by:

$\bar{k}=\frac{\sum_{i\in N}k_{i}}{n},$ (3)

where $1\leq i\leq n$ and $k_{i}$ is the degree of node $i$. For a complete
undirected graph, the average degree is $n-1$.

Clustering coefficient. The clustering coefficient describes the likelihood of
a node $j$ in the neighborhood $N_{i}$ of node $i$, is immediately connected
to other nodes in $N_{i}$ Watts and Strogatz (1998). The clustering
coefficient $C_{i}$ for node $i$ in an undirected graph is given by:

$C_{i}=\frac{|\\{e_{jk},j\in N_{i},k\in N_{i},e_{jk}\in E,j\neq
k\\}|}{\begin{pmatrix}k_{i}\\\ 2\end{pmatrix}},$ (4)

where $e_{jk}$ is the edge connecting node $j$ and $k$, $E$ is the set of
edges for graph $G$. Notation $|\\{e_{jk},j\in N_{i},k\in N_{i},e_{jk}\in
E,i\neq j\\}|$ represents the number of edges within the neighbourhood of node
$i$. A larger clustering coefficient physically means nodes in a graph are
more likely to cluster together. We use the average clustering coefficient
$\frac{1}{n}\sum_{i\in N}C_{i}$ as clustering coefficient $C$ of a graph.

The heterogeneity $\mathcal{H}$ of an unweighted graph is the ratio of
variance of degree $k$ to expectation of $k$ and given by Gao et al. (2016):

$\mathcal{H}=\frac{\frac{1}{n}\sum_{i\in N}k_{i}^{2}-(\frac{1}{n}\sum_{i\in
N}k_{i})^{2}}{\frac{1}{n}\sum_{i\in N}k_{i}}$ (5)

For a graph with homogeneous degree distribution, heterogeneity is equal to 0.

The average path length $\bar{l}$ measures the average distance between any
two nodes in the network and is defined by sum of shortest path length between
all pairs of nodes normalized by the total number of node pairs Albert and
Barabási (2002); Achard and Bullmore (2007):

$\bar{l}=\frac{\sum_{i,j\in N,i\neq j}l_{i,j}}{n(n-1)}$ (6)

Where $l_{i,j}$ is the shortest path length between node $i$ and $j$. The
average path length measures the efficiency of message passing over a graph.

The modularity describes the quality of partition of a graph into communities.
A good partition separates nodes in such way that majority of edges is in
communities and minority lies between them Clauset et al. (2004); Newman and
Girvan (2004). The modularity $Q$ is defined by:

$Q=\frac{1}{2n}\sum_{i,j\in
N}(A_{ij}-\frac{k_{i}k_{j}}{2n})\delta(c_{i},c_{j})$ (7)

Where $A$ is the adjacency matrix of a graph $G$, the $\delta$-function
$\delta(c_{i},c_{j})$ is 1 if $i$ and $j$ are in the same community 0
otherwise. Bimodularity $Q_{b}$ measures the quality of partitioning a graph
into two blocks using the Kernighan-Lin algorithm. Kernighan-Lin algorithm
partitions nodes of a graph in the manner of minimizing the costs on cutting
edges. Greedy modularity $Q_{g}$ measures the quality of partitioning nodes
using Clauset-Newman-Moore greedy modularity maximization Clauset et al.
(2004).

The resilience parameter $\beta_{\text{eff}}$ of an undirected graph is given
by Gao et al. (2016):

$\beta_{\text{eff}}=\frac{\frac{1}{n}\sum_{i\in
N}k_{i}^{2}}{\frac{1}{n}\sum_{i\in N}k_{i}}$ (8)

The physical meaning of resilience is the ability of a system to maintain
basic functionality after external perturbation.

The degree entropy $H$, a measure of disorder, calculates the entropy of the
degree distribution and is given by Ji et al. (2021):

$H=\frac{1}{n}\sum_{i\in N}-\frac{k_{i}}{m}\log\frac{k_{i}}{m}$ (9)

Where $m$ is the total number of edges.

The wedge count $W$ counts the number of wedge that is defined as a two-hop
path in an undirected graph. $W$ can be calculated by Ji et al. (2021); Gupta
et al. (2016):

$W=\sum_{i\in N}\begin{pmatrix}k_{i}\\\ 2\end{pmatrix}$ (10)

The Gini index $\mathcal{G}$ measures sparsity of a graph and is defined by
Goswami et al. (2018); Ji et al. (2021):

$\mathcal{G}=\frac{2\sum_{i\in N}i\hat{k}_{i}}{n\sum_{i\in
N}\hat{k}_{i}}-\frac{n+1}{n}$ (11)

Where $\hat{k}_{i}$ is the degree of node i after sorting degrees.

The average node connectivity $\bar{\kappa}$ is the average of local node
connectivity over all pairs of nodes of a graph $G$ and defined as Beineke et
al. (2002):

$\bar{\kappa}=\frac{\sum_{i,j\in N,i\neq j}\kappa(i,j)}{\begin{pmatrix}n\\\
2\end{pmatrix}}$ (12)

Where $\kappa(i,j)$ is defined as the maximum value of $\kappa$ for which node
$i$ and $j$ are $\kappa$-connected. For node $i$ and $j$ considered as
$\kappa$-connected, there are $\kappa$ or more pairwise internally disjoint
paths between them.

The edge connectivity is the minimum number of edges needed in order to
disconnect a graph $G$ Esfahanian (2013).

The closeness centrality $\mathcal{C}_{i}$ for a node $i$ is the reciprocal of
the average shortest path length $l$ of node $i$ to all other $n_{r}-1$ nodes
Freeman (1978):

$\mathcal{C}_{i}=\frac{n_{r}-1}{\sum_{j\in N,j\neq i}l_{i,j}}$ (13)

The improved closeness centrality $\mathcal{C}_{i}^{WF}$ by Wasserman and
Faust adds a scale factor to scale down closeness centrality for unconnected
graph and is given by Wasserman et al. (1994):

$\mathcal{C}_{i}^{WF}=\frac{n_{r}-1}{n-1}\frac{n_{r}-1}{\sum_{j\in N,j\neq
i}l_{i,j}}$ (14)

For a connected graph, there is no difference between the close centrality
$\mathcal{C}_{i}$ and the improved closeness centrality
$\mathcal{C}_{i}^{WF}$.

The eccentricity is defined as the maximum distance from node in a graph to
all other nodes Hage and Harary (1995). The diameter of a graph is the maximum
eccentricity while the radius of a graph is the minimum eccentricity.

The average edge betweenness centrality calculates the average of betweenness
centrality $\mathcal{C}_{B}(e_{k})$ for all edges of a graph. The betweenness
centrality for an edge $e_{k}$ is the sum of the fraction of all shortest
paths passing through $e_{k}$. The betweenness centrality can be calculated by
Brandes (2001, 2008):

$\mathcal{C}_{B}(e_{k})=\sum_{i,j\in N,i\neq
j}\frac{\sigma_{i,j}\bigr{\rvert}_{e_{k}}}{\sigma_{i,j}}$ (15)

Where $\sigma_{i,j}$ is the number of shortest paths connecting node $i$ and
$j$. $\sigma_{i,j}\bigr{\rvert}_{e_{k}}$ is the number of those paths passing
through edge $e_{k}$.

Similar to the average edge betweeness centrality, the average node
betweenness centrality calculates the average of betweenness centrality
$\mathcal{C}_{B}(n)$ for all nodes of a graph. $\mathcal{C}_{B}(n)$ is
expressed as:

$\mathcal{C}_{B}(k)=\sum_{i,j\in N,i\neq
j}\frac{\sigma_{i,j}\bigr{\rvert}_{k}}{\sigma_{i,j}}$ (16)

Where $\sigma_{i,j}\bigr{\rvert}_{k}$ is the number of shortest paths passing
through node $k$.

The central point of dominance $\mathcal{C}_{B}^{\prime}$ is expressed as
Brandes (2008):

$\mathcal{C}_{B}^{\prime}=\frac{\sum_{i\in N}(\max_{i\in
N}{C}_{B}(i)-{C}_{B}(i))}{n-1}$ (17)

The core number of a node is the largest value $k$ for a $k$-core subgraph
containing nodes of degree larger or equal to $k$.

The minimal Laplacian spectrum is the minimum and maximum eigenvalue of the
Laplacian matrix of a graph $G$.

The transitivity is defined as the ratio of the number of triangles to the
number of triads that are consisted of two edges with a shared node.

The efficiency of a pair of nodes is the inverse of the shortest path length
between these two nodes. The local efficiency $\mathcal{E}_{L}(i)$ is the
average efficiency of the neighbors $N_{i}$ of node $i$ and defined as Latora
and Marchiori (2001):

$\mathcal{E}_{L}(i)=\frac{1}{n(n-1)}\sum_{i,j\in N_{i},i\neq
j}\frac{1}{l_{i,j}}$ (18)

We use the average local efficiency of all nodes in a graph $G$ as the local
efficiency of $G$.

The global efficiency is the average efficiency of all pairs of nodes.

##  Appendix B MLP performance and Features

Figure 8 shows the projection of MLP top 1 error into each feature space. A
good fitting of single feature and top 1 error normally requires a
discontinuous function. But after a linear combination of those features can
give a good fitting quliaty.

##  Appendix C SFS after Fixing First Feature

Figure 9 shows the MSE variation as adding features. Figure 10 shows Pearson
correlation coefficient as adding features. Despite the difference in the
first feature, they all show a similar trend with adding more features.

Figure 8: Top 1 errors as a function of different features. Total number of
graphs is $19724$ Figure 9: MSE variation with adding features Figure 10:
Pearson correlation coefficient variation with adding features

## References

  * Achard and Bullmore [2007] Sophie Achard and Ed Bullmore.  Efficiency and cost of economical brain functional networks.  PLoS computational biology, 3(2):e17, 2007.
  * Albert and Barabási [2002] Réka Albert and Albert-László Barabási.  Statistical mechanics of complex networks.  Reviews of modern physics, 74(1):47, 2002.
  * Barabási and Albert [1999] Albert-László Barabási and Réka Albert.  Emergence of scaling in random networks.  science, 286(5439):509–512, 1999.
  * Beineke et al. [2002] Lowell W Beineke, Ortrud R Oellermann, and Raymond E Pippert.  The average connectivity of a graph.  Discrete mathematics, 252(1-3):31–45, 2002.
  * Brandes [2001] Ulrik Brandes.  A faster algorithm for betweenness centrality.  Journal of mathematical sociology, 25(2):163–177, 2001.
  * Brandes [2008] Ulrik Brandes.  On variants of shortest-path betweenness centrality and their generic computation.  Social Networks, 30(2):136–145, 2008.
  * Clauset et al. [2004] Aaron Clauset, Mark EJ Newman, and Cristopher Moore.  Finding community structure in very large networks.  Physical review E, 70(6):066111, 2004.
  * Cotter et al. [1999] Shane F Cotter, BD Rao, K Kreutz-Delgado, and J Adler.  Forward sequential algorithms for best basis selection.  IEE Proceedings-Vision, Image and Signal Processing, 146(5):235–244, 1999.
  * Elsken et al. [2019] Thomas Elsken, Jan Hendrik Metzen, and Frank Hutter.  Neural architecture search: A survey.  The Journal of Machine Learning Research, 20(1):1997–2017, 2019\.
  * Erdos et al. [1960] Paul Erdos, Alfréd Rényi, et al.  On the evolution of random graphs.  Publ. Math. Inst. Hung. Acad. Sci, 5(1):17–60, 1960.
  * Esfahanian [2013] Abdol-Hossein Esfahanian.  Connectivity algorithms.  In Topics in structural graph theory, pages 268–281. Cambridge University Press, 2013.
  * Freeman [1978] Linton C Freeman.  Centrality in social networks conceptual clarification.  Social networks, 1(3):215–239, 1978.
  * Gao et al. [2016] Jianxi Gao, Baruch Barzel, and Albert-László Barabási.  Universal resilience patterns in complex networks.  Nature, 530(7590):307–312, 2016.
  * Goswami et al. [2018] Swati Goswami, CA Murthy, and Asit K Das.  Sparsity measure of a network graph: Gini index.  Information Sciences, 462:16–39, 2018.
  * Gupta et al. [2016] Rishi Gupta, Tim Roughgarden, and Comandur Seshadhri.  Decompositions of triangle-dense graphs.  SIAM Journal on Computing, 45(2):197–215, 2016.
  * Hage and Harary [1995] Per Hage and Frank Harary.  Eccentricity and centrality in networks.  Social networks, 17(1):57–63, 1995.
  * Harary [1962] Frank Harary.  The maximum connectivity of a graph.  Proceedings of the National Academy of Sciences of the United States of America, 48(7):1142, 1962.
  * He et al. [2016] Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun.  Deep residual learning for image recognition.  In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 770–778, 2016.
  * Hinton et al. [2012] Geoffrey Hinton, Li Deng, Dong Yu, George E Dahl, Abdel-rahman Mohamed, Navdeep Jaitly, Andrew Senior, Vincent Vanhoucke, Patrick Nguyen, Tara N Sainath, et al.  Deep neural networks for acoustic modeling in speech recognition: The shared views of four research groups.  IEEE Signal processing magazine, 29(6):82–97, 2012.
  * Ioffe and Szegedy [2015] Sergey Ioffe and Christian Szegedy.  Batch normalization: Accelerating deep network training by reducing internal covariate shift.  In International conference on machine learning, pages 448–456. PMLR, 2015.
  * Ji et al. [2021] Yuliang Ji, Ru Huang, Jie Chen, and Yuanzhe Xi.  Generating a doppelganger graph: Resembling but distinct.  arXiv preprint arXiv:2101.09593, 2021.
  * Krizhevsky et al. [2009] Alex Krizhevsky, Geoffrey Hinton, et al.  Learning multiple layers of features from tiny images.  Citeseer, 2009.
  * Krizhevsky et al. [2012] Alex Krizhevsky, Ilya Sutskever, and Geoffrey E Hinton.  Imagenet classification with deep convolutional neural networks.  Advances in neural information processing systems, 25:1097–1105, 2012.
  * Latora and Marchiori [2001] Vito Latora and Massimo Marchiori.  Efficient behavior of small-world networks.  Physical review letters, 87(19):198701, 2001.
  * Loshchilov and Hutter [2016] Ilya Loshchilov and Frank Hutter.  Sgdr: Stochastic gradient descent with warm restarts.  arXiv preprint arXiv:1608.03983, 2016.
  * Luce and Perry [1949] R Duncan Luce and Albert D Perry.  A method of matrix analysis of group structure.  Psychometrika, 14(2):95–116, 1949.
  * Marcano-Cedeño et al. [2010] Alexis Marcano-Cedeño, J Quintanilla-Domínguez, MG Cortina-Januchs, and Diego Andina.  Feature selection using sequential forward selection and classification applying artificial metaplasticity neural network.  In IECON 2010-36th annual conference on IEEE industrial electronics society, pages 2845–2850. IEEE, 2010.
  * Nesterov [1983] Yurii E Nesterov.  A method for solving the convex programming problem with convergence rate o (1/k^ 2).  In Dokl. akad. nauk Sssr, volume 269, pages 543–547, 1983.
  * Newman and Girvan [2004] Mark EJ Newman and Michelle Girvan.  Finding and evaluating community structure in networks.  Physical review E, 69(2):026113, 2004.
  * Prettejohn et al. [2011] Brenton J Prettejohn, Matthew J Berryman, and Mark D McDonnell.  Methods for generating complex networks with selected structural properties for simulations: a review and tutorial for neuroscientists.  Frontiers in computational neuroscience, 5:11, 2011.
  * Real et al. [2019] Esteban Real, Alok Aggarwal, Yanping Huang, and Quoc V Le.  Regularized evolution for image classifier architecture search.  In Proceedings of the aaai conference on artificial intelligence, volume 33, pages 4780–4789, 2019.
  * Ru et al. [2020] Binxin Ru, Pedro Esperanca, and Fabio Carlucci.  Neural architecture generator optimization.  arXiv preprint arXiv:2004.01395, 2020.
  * Simonyan and Zisserman [2014] Karen Simonyan and Andrew Zisserman.  Very deep convolutional networks for large-scale image recognition.  arXiv preprint arXiv:1409.1556, 2014.
  * Sutskever et al. [2013] Ilya Sutskever, James Martens, George Dahl, and Geoffrey Hinton.  On the importance of initialization and momentum in deep learning.  In International conference on machine learning, pages 1139–1147. PMLR, 2013.
  * Ververidis and Kotropoulos [2005] Dimitrios Ververidis and Constantine Kotropoulos.  Sequential forward feature selection with low computational cost.  In 2005 13th European Signal Processing Conference, pages 1–4. IEEE, 2005.
  * Wasserman et al. [1994] Stanley Wasserman, Katherine Faust, et al.  Social network analysis: Methods and applications.  Cambridge university press, 1994.
  * Watts and Strogatz [1998] Duncan J Watts and Steven H Strogatz.  Collective dynamics of ‘small-world’networks.  nature, 393(6684):440–442, 1998.
  * Wu et al. [2016] Yonghui Wu, Mike Schuster, Zhifeng Chen, Quoc V Le, Mohammad Norouzi, Wolfgang Macherey, Maxim Krikun, Yuan Cao, Qin Gao, Klaus Macherey, et al.  Google’s neural machine translation system: Bridging the gap between human and machine translation.  arXiv preprint arXiv:1609.08144, 2016.
  * You et al. [2020] Jiaxuan You, Jure Leskovec, Kaiming He, and Saining Xie.  Graph structure of neural networks.  In International Conference on Machine Learning, pages 10881–10891. PMLR, 2020.
  * Zoph et al. [2018] Barret Zoph, Vijay Vasudevan, Jonathon Shlens, and Quoc V Le.  Learning transferable architectures for scalable image recognition.  In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 8697–8710, 2018.


# Contact effects on thermoelectric properties of textured graphene
nanoribbons

David M T Kuo1 and Yia-Chung Chang2,3 1Department of Electrical Engineering
and Department of Physics, National Central University, Chungli, 320 Taiwan
2Research Center for Applied Sciences, Academic Sinica, Taipei, 11529 Taiwan
3 Department of Physics, National Cheng Kung University, Tainan, 701 Taiwan

###### Abstract

Transport and thermoelectric properties of finite textured graphene
nanoribbons (t-GNRs) connected to electrodes with various coupling strengths
are theoretically studied in the framework of the tight-binding model and
Green’s function approach. Due to quantum constriction induced by the indented
edges, such t-GNRs behave like serially-coupled graphene quantum dots (SGQDs).
These types of SGQDs can be formed by tailoring zigzag GNRs (ZGNRs) or
armchair GNRs (AGNRs). Their bandwidths and gaps can be engineered by varying
the size of the quantum dot and the neck width at indented edges. Effects of
defects and contact junction on electrical conductance, Seebeck coefficient
and electron thermal conductance of t-GNRs are calculated. When a defect
occurs in the interior site of textured ZGNRs (t-ZGNRs), the maximum power
factor within the central gap or near the band edges is found to be
insensitive to the defect scattering. Furthermore, we found that SGQDs formed
by t-ZGNRs have significantly better electrical power outputs than those of
textured ANGRs due to the improved functional shape of the transmission
coefficient in t-ZGNRs. With a proper design of contact the maximum power
factor ( figure of merit) of t-ZGNRs could reach $90\%$ ($95\%$) of the
theoretical limit.

##  I Introduction

Due to global warming, the Kyoto protocol aiming to reduce $CO_{2}$ emissions
was proposed in 1997. Since then, renewable energies including solar, wind,
rain, tides, and geothermal heat are topics of tremendous scientific interest
[ChenG, 1]. Thermoelectric devices can be used as power generators and
refrigerators. The electrical power based on thermoelectric effects is one of
the most important type of green energies [ChenG, 1-Heremans, 6]. Because of
the inherent physics of thermoelectric effect, thermoelectric devices can
fully avoid $CO_{2}$, hydrofluorocarbons, and perfluorocarbons emissions. In
addition, thermoelectric devices can sustain long operation time and avoid
mechanical noises.

Recently, energy harvesting applications in nanoscale systems have attracted
considerable attention[Nika, 7-Luo, 12]. How to obtain the maximum
thermoelectric efficiency of heat engines with optimized electrical power
output has been a key issue [Whitney, 9-Luo, 12]. To improve the performance
of a heat engine it is preferable to have the electron transport in ballistic
regime and phonon transport in the diffusive scattering regime. Therefore, a
high-performance thermoelectric device needs to provide a channel length
shorter than the electron mean free path ($\lambda_{e}$) but much longer than
the phonon mean free path ($\lambda_{ph}$) in order to reduce the ratio of
phonon thermal conductance ($\kappa_{ph}$) to electron thermal conductance
($\kappa_{e}$) [Xu, 13]. It was pointed out in ref. [Whitney, 9] that a Carnot
heat engine favors the electron transport in an energy range where the
transmission coefficient has a steep change with respect to energy, e.g. with
a square-form (SF). Up to now it remains unclear how to realize a SF
transmission coefficient in realistic thermoelectric devices with a short
channel length between thermal contacts, [Darancet, 14,ChenIJ, 11].

Considerable scientific efforts paved the way to answer such an intriguing
problem. Hicks and Dresselhaus theoretically demonstrated that the
thermoelectric performance can be significantly enhanced in one-dimensional
(1D) nanowires due to the reduced phonon thermal conductance and the enhanced
Seebeck coefficient ($S$) [Hicks, 2,Dresselhaus, 3]. According to theoretical
modeling, $\kappa_{ph}$ of 1D silicon quantum dot (QD) superlattices can be
reduced by one order of magnitude in comparison with 1D nanowires [Nika,
7,HuM, 8]. Due to limitations in technology the dot-size fluctuation in the 1D
silicon quantum-dot (QD) arrays still remains a serious issue. [Kagan,
15,Lawrie, 16]

The discovery of graphene in 2004 opened the door for realizing 1D nanowires
with small cross-section to degrade $\kappa_{ph}$ [Novoselovks, 17], since one
can fabricate graphene nanoribbons (GNRs) with atomic precision via the
bottom-up approach[Cai, 18]. This approach has been successfully applied to
build more complex systems, such as armchair GNRs with periodically corrugated
edges, called as textured AGNRs (t-AGNRs) here.[ChenYC, 19-DJRizzo, 25] The
scanning tunneling microscopy (STM) spectra of serially-coupled graphene
quantum dots (SGQDs) synthesized by t-AGNRs were reported experimentally
[DJRizzo, 25]. Novel graphene-based electronic devices have also focused on
AGRNs and t-AGNRs [Llinas, 26]. Recently, the existence of edge states with
topological protection in textured zigzag GNRs (t-ZGNRs) has been demonstrated
theoretically [Chou, 27].

Several experimental studies of thermoelectric properties of graphene-related
materials have been reported in recent years [ZuevYM, 28-WangYY, 34].
Nevertheless, there is a paucity of studies to consider the contact effect on
thermoelectric properties of t-AGNRs and t-ZGNRs[HeJ, 35]. In this work, our
goal is to optimize the transmission coefficient of t-GNRs with length shorter
than $\lambda_{e}$, but larger than $\lambda_{ph}$ under different coupling
strengths with the electrical contact. We found that minibands can be formed
in these t-GNRs, while their energy gaps can be tuned by varying the size of
GQDs and the inter-dot coupling strength.Furthermore, the calculated power
factor is very robust against scattering from defects occurring inside the
interior of these t-ZGNRs. The transmission coefficients through the minibands
of t-ZGNRs under optimized tunneling rates provide desired characteristics
showing steep change with respect to energy near the central gap that mimics
the theoretical limit obtained by using a SF transmission coefficient. Due to
this feature, our calculations reveal that the thermoelectric performance of
t-ZGNRs can be significantly better than the t-AGNRs. The maximum power factor
(figure of merit) of t-ZGNRs could reach $90\%$ ($95\%$) of that obtained with
the SF transmission coefficient. Therefore, the t-ZGNR is a promising
candidate for applications in nanoscale energy harvesting.

Figure 1: Schematic diagrams of (a) ZGNR, (b) textured-ZGNR (t-ZGNR), (c)
AGNR, and (d) t-AGNR. In (a), we also show the electrodes with $\Gamma_{L}$
($\Gamma_{R}$) denoting the tunneling rate of the electrons between the left
(right) electrode and the leftmost (rightmost) atoms of the ZGNR. $T_{L}$ and
$T_{R}$ denote the temperature of the left ($L$) and the right ($R$)
electrodes, respectively. Note that t-ZGNRs can be formed by periodically
removing some carbon atoms on the zigzag edges of ZGNRs. In (b) the t-ZGNR
consists of GQDs with size characterized by $(N_{a},N_{z})=(8,7)$. In (c) the
AGNR has a length $N_{a}=36$ and width $N_{z}=5$. In (d), the width and length
of the unit cell for the t-AGNR superlattice are characterized by $N_{z}=5$
and $N_{a}=8$, respectively.

##  II Calculation method

To model the thermoelectric properties of t-GNRs connected to the electrodes,
it is a good approximation to employ a tight-binding model with one $p_{z}$
orbital per atomic site to describe the electronic states.[Nakada,
36-Sevincli, 40] The Hamiltonian of the nano-junction system depicted in Fig.
1, including four different GNR structures, can be written as
$H=H_{0}+H_{GNR}$ [Haug, 41], where

$\displaystyle H_{0}$ $\displaystyle=$
$\displaystyle\sum_{k}\epsilon_{k}a^{\dagger}_{k}a_{k}+\sum_{k}\epsilon_{k}b^{\dagger}_{k}b_{k}$
$\displaystyle+$
$\displaystyle\sum_{\ell}\sum_{k}V^{L}_{k,\ell,j}d^{\dagger}_{\ell,j}a_{k}+\sum_{\ell}\sum_{k}V^{R}_{k,\ell,j}d^{\dagger}_{\ell,j}b_{k}+h.c.$

The first two terms of Eq. (1) describe the free electrons in the left ($L$)
and right ($R$) electrodes. $a^{\dagger}_{k}$ ($b^{\dagger}_{k}$) creates an
electron with wave number $k$ and energy $\epsilon_{k}$ in the left (right)
electrode. $V^{L}_{k,\ell,j=1}$ ($V^{R}_{k,\ell,j=N_{z}(N_{a})}$) describes
the coupling between the left (right) lead with its adjacent atom in the
$\ell$-th row.

$\displaystyle H_{GNR}$ $\displaystyle=$
$\displaystyle\sum_{\ell,j}E_{\ell,j}d^{\dagger}_{\ell,j}d_{\ell,j}$
$\displaystyle-$
$\displaystyle\sum_{\ell,j}\sum_{\ell^{\prime},j^{\prime}}t_{(\ell,j),(\ell^{\prime},j^{\prime})}d^{\dagger}_{\ell,j}d_{\ell^{\prime},j^{\prime}}+h.c,$

where $E_{\ell,j}$ is the on-site energy for the $p_{z}$ orbital in the
${\ell}$-th row and $j$-th column. Here, the spin-orbit interaction is
neglected. $d^{\dagger}_{\ell,j}(d_{\ell,j})$ creates (destroys) one electron
at the atom site labeled by ($\ell$,$j$) where $\ell$ and $j$, respectively
are the row and column indices as illustrated in Fig. 1.
$t_{(\ell,j),(\ell^{\prime},j^{\prime})}$ describes the electron hopping
energy from site ($\ell$,$j$) to site ($\ell^{\prime}$,$j^{\prime}$). The
tight-binding parameters used for GNRs is $E_{\ell,j}=0$ for the on-site
energy and $t_{(\ell,j),(\ell^{\prime},j^{\prime})}=t_{pp\pi}=2.7$ eV for
nearest-neighbor hopping strength.

To study the transport properties of a GNR junction connected to electrodes,
it is convenient to use the Keldysh Green’s function technique [Haug, 41].
Electron and heat currents leaving the electrodes can be expressed as

$\displaystyle J$ $\displaystyle=$
$\displaystyle\frac{2e}{h}\int{d\varepsilon}~{}{\cal
T}_{LR}(\varepsilon)[f_{L}(\varepsilon)-f_{R}(\varepsilon)],$ (3)

and

$Q_{e,L(R)}=\frac{\pm 2}{h}\int{d\varepsilon}~{}{\cal
T}_{LR}(\varepsilon)(\varepsilon-\mu_{L(R)})[f_{L}(\varepsilon)-f_{R}(\varepsilon)]$
(4)

where
$f_{\alpha}(\varepsilon)=1/\\{\exp[(\varepsilon-\mu_{\alpha})/k_{B}T_{\alpha}]+1\\}$
denotes the Fermi distribution function for the $\alpha$-th electrode, where
$\mu_{\alpha}$ and $T_{\alpha}$ are the chemical potential and the temperature
of the $\alpha$ electrode. $e$, $h$, and $k_{B}$ denote the electron charge,
the Planck’s constant, and the Boltzmann constant, respectively. ${\cal
T}_{LR}(\varepsilon)$ denotes the transmission coefficient of a GNR connected
to electrodes, which can be solved by the formula ${\cal
T}_{LR}(\varepsilon)=4Tr[\Gamma_{L}(\varepsilon)G^{r}(\varepsilon)\Gamma_{R}(\varepsilon)G^{a}(\varepsilon)]$
[Kuo4, 42,Phung, 43], where $\Gamma_{L}(\varepsilon)$ and
$\Gamma_{R}(\varepsilon)$ denote the tunneling rate (in energy units) at the
left and right leads, and ${G}^{r}(\varepsilon)$ and ${G}^{a}(\varepsilon)$
are the retarded and advanced Green’s functions of the GNR, respectively. The
tunneling rates are described by the imaginary part of self-energy correction
on the interface atoms of the GNR due to the coupling with nearby atoms in the
adjacent electrodes, i.e.
$\Gamma_{L(R)}(\varepsilon)=-Im(\Sigma^{r}_{L(R)}(\varepsilon)$). Such a
coupling depends on the contact quality with the electrodes, which is
characterized by the interaction strength
$V^{L}_{k,\ell,j=1}(V^{R}_{k,\ell,j=N_{z}(N_{a})})$) with the left (right)
lead. Here, we have adopted energy-independent tunneling rates
$\Gamma_{L(R)}(\varepsilon)=\Gamma_{L(R)}$ which is reasonable in the wide-
band limit for the electrodes [Kuo4, 42].Note that $\Gamma_{\alpha}$ and
Green’s functions are matrices in the basis of tight-binding orbitals.
$\Gamma_{\alpha}$ for the boundary atoms have diagonal entries given by the
same constant $\Gamma_{t}$. Because of the line contact as illustrated in Fig.
1, the contact resistance can be much smaller than that of surface contact
[ChenRS, 44]. Meanwhile, the variation of tunneling rates could reveal
different contact properties such as the Schottky barrier or ohmic contact
[Mahan, 45].

In the linear response regime, the electrical conductance ($G_{e}$), Seebeck
coefficient ($S$) and electron thermal conductance ($\kappa_{e}$) are given by
$G_{e}=e^{2}{\cal L}_{0}$, $S=-{\cal L}_{1}/(eT{\cal L}_{0})$ and
$\kappa_{e}=\frac{1}{T}({\cal L}_{2}-{\cal L}^{2}_{1}/{\cal L}_{0})$ with
${\cal L}_{n}$ ($=0,1,2$) defined as

${\cal L}_{n}=\frac{2}{h}\int d\varepsilon~{}{\cal
T}_{LR}(\varepsilon)(\varepsilon-\mu)^{n}\frac{\partial
f(\varepsilon)}{\partial\mu}.$ (5)

Here $f(\varepsilon)=1/(exp^{(\varepsilon-\mu)/k_{B}T}+1)$ is the Fermi
distribution function of electrodes at equilibrium temperature $T$ and
chemical potential $\mu$. As seen in Eq. (5), the transmission coefficient
${\cal T}_{LR}(\varepsilon)$ plays a significant role for electron transport
between the left($L$) and right ($R$) electrodes. At zero temperature, the
electrical conductance is given by $G_{e}(\mu)=\frac{2e^{2}}{h}{\cal
T}_{LR}(\mu)$. In the current study, only even numbers of $N_{a}$ are
considered to avoid unwanted dangling-bond states. We have developed an
efficient computation method that allows us to calculate the Green’s functions
of large-sized quantum structures. For the current study, textured ZGNRs and
AGNRs with lengths up to $14nm$ are considered.

##  III Results and discussion

###  III.1 Graphene nanoribbons

Topological states offer promising applications in electronics and
optoelectronics, owing to their robustness in transport characteristics
against defect scattering. Many studies have confirmed that 2D and 1D
topological states (TSs) exist in certain material structures [ChangYC,
46-Xian, 51]. For instance, there exists a 2D topologically protected
interface state in HgTe/CdTe superlattices [ChangYC, 46, 47]. One-dimensional
(1D) TSs were theoretically predicted to exist in square selenene and
tellurene [Xian, 51]. Recently, zero-dimensional (0D) TSs of finite-size GNRs
have been extensively studied [ChenYC, 19-DJRizzo, 25] because the 0D TS
offers more flexibility in the design of complicated electronic circuits [Yan,
22]. Before illustrating electron coherent transport through SGQDs formed by
t-ZGNRs, we first examine the characteristics of 0D TSs of finite-size GNRs by
calculating their transmission coefficient, ${\cal T}_{LR}(\varepsilon)$.

Figure 2 shows the calculated electron conductance, $G_{e}$ at $k_{B}T=0$ as a
function of $\mu$ for various $N_{a}$ with $N_{z}=7$. We note that there are
two zigzag edge states localized at the top and bottom sides of the GQDs
depicted in Fig. 1(a). For $N_{a}=40$, their wave functions are well separated
along the armchair direction [Wakabayashi, 37]. The electrical conductance
spectrum shows that $G_{e}=2G_{0}$ at $\mu=0$ and $G_{e}=G_{0}$ for other
electronic states, where $G_{0}=\frac{2e^{2}}{h}$ is the quantum conductance.
Such zero energy modes were observed experimentally by STM [DRizzo, 23]. When
two electronic states are closely spaced (with energy separation less than the
broadening), $G_{e}$ can become large than $G_{0}$. As $N_{a}$ decreases,
$G_{e}$ for the zero-energy mode is split into two peaks corresponding to the
bonding and antibonding states of coupled zigzag edge states, as seen in Fig.
2c. When $N_{a}=12$, these two peaks are well separated. Here, $\epsilon_{HO}$
and $\epsilon_{LU}$ denote the energy levels of the highest occupied molecular
orbital (HOMO) and the lowest unoccupied molecular orbital (LUMO),
respectively. Due to the short decay lengths along the armchair edge
directions for zigzag edges states, ${\cal T}_{LR}(\varepsilon)$ of
$\Sigma_{0}$ also depends on whether the zigzag or armchair edges are coupled
to the electrodes. In Ref. [MahanG, 52] it is proposed that a single QD can be
used to realize a Carnot heat engine. However, the channel length for
$N_{z}=7$ ($L_{z}=0.738~{}nm$) is too small to avoid the serious degradation
of the figure of merit due to the phonon heat conductivity, $\kappa_{ph}$. On
the other hand, for the large $N_{z}$ case ($N_{z}\gg N_{a}$), a finite GNR
shows metallic behavior, leading to unfavorable thermoelectric properties. To
maintain a sizable gap ($~{}10k_{B}T$) around the charge neutrality point
(CNP) while keeping $\kappa_{ph}$ small enough to preserve decent figure of
merit ($ZT=\frac{S^{2}G_{e}T}{\kappa_{e}+\kappa_{ph}}$), the t-ZGNR becomes a
suitable candidate as we shall discuss below.

Figure 2: Electrical conductance, $G_{e}$ of a finite GNR with armchair edges
coupled to the electrodes at $k_{B}T=0$ as a function of the chemical
potential, $\mu$ for various widths with $N_{a}=40,32,12,$ and 8 with $N_{z}$
fixed at 7. We have adopted electron tunneling rates $\Gamma_{t}=90$ meV.

###  III.2  SGQDs formed by textured ZGNRs

For energy harvesting applications at room temperature, we need to design a
${\cal T}_{LR}(\varepsilon)$ spectrum with a square shape near the central gap
to obtain optimized $ZT$ and electrical power output [Whitney, 9]. Let’s
consider an SGQD formed by t-ZGNR. A t-ZGNR can be realized by tailoring GNR
with periodic indentation on the zigzag edges such as the structure shown in
Fig.1(b). This t-ZGNR consists of GQDs with size characterized by $N_{a}=12$
and $N_{z}=7$. A single GQD of this size has resonance energies
$\varepsilon_{LU(HO)}=\pm 0.247eV$ near the CNP, as shown in Fig. 2(c). For
two coupled GQDs (2GQDs), there is one satellite peak on the left (right) side
of the $\varepsilon_{HO(LU)}$ peak (see Fig. 3(b)). For six coupled GQDs
(6GQDs), we have five additional peaks with $\varepsilon_{e(h),1}=\pm
0.3105eV$, $\varepsilon_{e(h),2}=\pm 0.468eV$, $\varepsilon_{e(h),3}=\pm
0.6705eV$,$\varepsilon_{e(h),4}=\pm 0.882eV$ and $\varepsilon_{e(h),5}=\pm
1.0665eV$ on the right (left) side of $\varepsilon_{LU(HO)}=\pm 0.247eV$ (see
Fig. 3(f)). We note that $G_{e}(\mu)$ is fully suppressed for $\mu$ between
$\varepsilon_{HO}$ and $\varepsilon_{LU}$. Therefore, SGQDs formed by t-ZGNRs
function as filters to block electrons with energies near the CNP, which is
different from that of serially coupled antidots realized by ZGNRs with
nanopores [ChangPH, 53]. The separation between peaks is inhomogeneous. In
addition, the $\varepsilon_{HO}$ and $\varepsilon_{LU}$ peaks become sharper
with increasing GQD number. These features are attributed to a special parity
symmetry in the supercell of ZGNRs. Due to the coupling of transverse and
longitudinal wave numbers in ZGNRs, the density distributions of electronic
state near CNP are inhomogeneous [Wakabayashi2, 54]. AGNRs do not have such a
parity in the supercell. The feature provides a significant effect on the
transmission coefficient of SGQDs formed by t-ZGNRs. To further illustrate the
characteristics of the edge states we show the charge densities of states with
$\varepsilon_{LU}$, $\varepsilon_{e,1}$ and $\varepsilon_{e,2}$ for the case
of 6GQDs in Fig. A.1 of the appendix.

Figure 3: Electrical conductance, $G_{e}$ as a function of $\mu$ for SGQDs
with various of GQD numbers. A SGQD consists of GQDs with the size
characterized by $N_{a}=12$ and $N_{z}=7$. Other physical parameters
correspond to those of Fig. 2(c).

Previous theoretical studies have demonstrated that edge defects can
significantly reduce the electron quantum conductance of ZGNRs [Areshkin,
55,Martins, 56]. Here, we investigate how defects (either at the edge or in
the interior region) influence the thermoelectric characteristics of t-ZGNRs
by introduce energy shifts $\Delta_{\ell,j}$ on designated defect sites.
$\Delta_{\ell,j}$ could be positive or negative, which depends on the kind of
defects [Martins, 56]. To model a vacancy within a tight binding model, one
typically takes $\Delta_{\ell,j}\rightarrow\infty$. The larger orbital-energy
shift, the stronger effect on the electrical conductance can be seen [LiTC,
57]. Here, we consider the case of a positive and large $\Delta$ to
investigate the effects of defects on the electron transport of t-ZGNRs. We
calculate $G_{e}$, $\kappa_{e}$, $S$ and power factor $PF=S^{2}G_{e}$ for
different defect locations for the case of 15 GQDs ($L_{z}=14.5$nm or
$N_{z}=119$) with $\Gamma_{t}=0.54eV$ and show the results in Fig. 4.

As seen in Fig. 4(a) and 4(b), the size of the central gap and the shape of
$G_{e}$ and $\kappa_{e}$ spectra (green lines) are changed only slightly when
the defects are located in the interior region. However the $G_{e}$ and
$\kappa_{e}$ values near the central gap become seriously suppressed (as shown
by red lines) when a single defect is located at a zigzag edge site labeled by
$(1,4)$, where the charge density of the $\varepsilon_{LU(HO)}$ electronic
state is peaked. As seen in Fig. 4(c) the antisymmetric pattern of the Seebeck
coefficient $S(\mu)$ (with respect to the sign change of $\mu$) due to the
electron-hole symmetry is distorted in the presence of defects. The peak value
of the power factor $PF=S^{2}G_{e}$ near the central gap is seriously serious
reduced when defects occur on the edge, but much less affected by defects in
the interior region as illustrated in Fig. 4(d). Here and henceforth,
$\kappa_{e}$ is in units of $\kappa_{0}=0.62nW/K$, $S$ is expressed in units
of $k_{B}/e=86.25\mu V/K$, and the power factor ($PF$) in units of
$2k^{2}_{B}/h=0.575pW/K^{2}$. It is remarkable to see that $PF$ is very robust
against defect scattering when defects are away from zigzag edges. Such
location-dependent effects can be depicted by the charge density distribution
in Fig. A1. To reduce defect effect on electron transport, one needs to avoid
creating defects randomly located on the zigzag edges.

Figure 4:  (a) Electrical conductance $G_{e}$, (b) electron thermal
conductance $\kappa_{e}$, (c) Seebeck coefficient $S$, and (d) power
factor$PF=S^{2}G_{e}$ as functions of $\mu$ for defect locations at
$k_{B}T=27meV$. The tunneling rates used are $\Gamma_{t}=0.54eV$. The length
of SGQD is $L_{z}=14.5nm$ ($N_{z}=119$). Each GQD in the SGQD structure has
size $N_{a}=12$ and $N_{z}=7$.

We note that SGQDs are formed by periodically removing some carbon atoms on
the zigzag edges of ZGNRs. To understand the relationship between the quantum
states near the central gap in the SGQD formed by t-ZGNRs and the quantum
states near the zero energy modes of ZGNRs, we show a comparison of subband
structures of the infinite long t-ZGNR shown in Fig. 1(b) and unaltered ZGNR
shown in Fig. 1(a) (with an enlarged super cell of length $L=4a$ to match the
unit cell of t-ZGNR) in Fig. A.2(a) of the appendix. Note that the first
Brillouin zone (BZ) has one quarter the size of the BZ of ZGNR. Thus, the
subband structures shown in Fig. A.2(a) includes the zone-folding effect that
maps the zone boundary ($k_{z}=\pi/a$) of the unfolded BZ of ZGNR to $k_{z}=0$
of the folded BZ here. As expected, the unaltered ZGNR has edge states with
zero energy for $k_{z}<0.06\pi/a=0.24\pi/L$. The two edge states are strongly
coupled in the t-ZGNR with an energy splitting of $~{}2.5eV$. On the other
hand, the edge states at mini-zone boundary ($k_{z}=\pi/L$) with energies
corresponding to the HOMO (LUMO) level are nearly unaffected by the removal of
a carbon atom in the neck region of t-ZGNR, since the wave functions of these
edge states have a node at that position. Theoretical calculations of the Zak
phases of various t-ZGNRs indicate that the topological nature of edge states
of t-ZGNRs are preserved even though their energy levels are shifted away from
the CNP [Chou, 27]. Thus, these states are topologically protected and the
transport through these states should be rather insensitive to the presence of
defects inside the t-ZGNR.

Junction tunneling rates can be affected by the quality of contact between the
electrode and semiconductor,[ChenRS, 44] which is a critical issue for device
applications of two-dimensional materials.[Shen, 58] To clarify the contact
effect, we show in Fig. 5 the calculated $G_{e}$, $\kappa_{e}$, $S$, and $ZT$
at room temperature for various tunneling rates ($\Gamma_{t}$) for a t-ZGNR
with 15 GQDs, each having the dimension $N_{a}=8$ and $N_{z}=7$ as shown in
Fig. 1(b). When the tunneling rate increases, the transmission coefficient
${\cal T}_{LR}(\epsilon)$ through minibands close to HOMO (LUMO) are enhanced.
For the case of 15 GQDs, there are 15 peaks in each miniband (see Fig. A.3).
Such an enhancement leads to an increase of both $G_{e}$ and $\kappa_{e}$ as
seen in Fig. 5(a) and 5(b). However, $S$ is essentially independent to the
variation of $\Gamma_{t}$ as indicated by the collapsing of all four curves in
Fig. 5(c). As a result, the power factor $PF=S^{2}G_{e}$ also increases with
$\Gamma_{t}$. It is noted that the maximum $S$ reaches $1.51mV/K$[see Fig
5(c)], which is much larger than that observed experimentally in gapless
graphene [ZuevYM, 28,WeiP, 29]. In Fig. 5(d), the maximum $ZT$ occurs at
$\mu=\pm 0.423eV$, where $\kappa_{ph}\gg\kappa_{e}$. Therefore, the
enhancement of $ZT$ with respect to the increase of $\Gamma_{t}$ mainly arises
from the enhancement of $PF$. Here, we have included the effect of phonon
thermal conductance, $\kappa_{ph}$, and assumed
$\kappa_{ph}=F_{s}\kappa^{0}_{ph}$ in the calculation of $ZT$ in Fig. 5(d),
where $\kappa^{0}_{ph}=\frac{\pi^{2}k^{2}_{B}T}{3h}$ is the phonon thermal
conductance of an ideal ZGNR. We adopt the room-temperature value of
$\kappa^{0}_{ph}=0.285nW/K$ for ZGNR with width $N_{a}\leq 8$ obtained by a
first-principles calculation as given in ref[Zhengh, 59]. $F_{s}=0.1$ denotes
a reduction factor resulting from quantum constriction in t-ZGNRs. It has been
theoretically demonstrated that the magnitude of $\kappa_{ph}$ can be reduced
by one order magnitude for ZGNRs with ”edge vacancies” ref[Cuniberti, 60],
which is similar to the mechanism introduced in silicon nanowires with surface
roughness[Murphy, 61]. The measured $\lambda_{ph}$ can be reduced from
$300-600nm$ in a single layer graphene to $10nm$ in graphene nanostructures
(see [XuY, 30] and references therein). Recently, very short $\lambda_{ph}$
has been reported experimentally, which offers promise for enhancing the
figure of merit ($ZT$) of graphene heterostructures[WangYH, 31-WangYY, 34].

The calculated results shown in Fig. 5 imply that $ZT>3$ could be realized by
using SGQDs formed by t-ZGNRs for tunneling rates corresponding to
$\Gamma_{t}$ between $t_{pp\pi}/10=0.27$ eV and $t_{pp\pi}/5=0.54$ eV. In Fig.
A.3 of the appendix, we show the peak value of $ZT(ZT_{max})$ as a function of
tunneling rate $\Gamma_{t}$, and we found that the $ZT_{max}$ can be larger
than 3 for $\Gamma_{t}$ between $0.18$ eV and $1.45$ eV. We note that if a
semi-infinite ZGNR contact is used, it is reasonable to assume that the
coupling strength between the semi-infinite ZGNR and t-ZGNR can be comparable
to the hopping strength $t_{pp\pi}=2.7eV$. Here, we show that the maximum $ZT$
can reach 3.7 with $\Gamma_{t}=0.54eV=t_{pp\pi}/5$ that could be achievable by
a good contact.

Figure 5: (a) Electrical conductance $G_{e}$, (b) electron thermal conductance
$\kappa_{e}$, (c) Seebeck coefficient $S$, and (d) figure of merit $ZT$ as
functions of $\mu$ for various tunneling rates ($\Gamma_{t}$) at
$k_{B}T=27meV$. The t-ZGNR structure with $L_{z}=14.5nm$ considered is
illustrated in Fig. 1(b). The size of each GQD in the t-ZGNR is characterized
by $N_{a}=8$ and $N_{z}=7$.

Because $S$ is a robust physical quantity with respect to the variation of
tunneling rate and channel length (which only depend on the magnitude of
central gap and temperature in our case), the optimization of $G_{e}$ becomes
a critical issue in SGQDs when $\kappa_{ph}\gg\kappa_{e}$. To provide better
understanding of the effect of tunneling rate on thermoelectric properties, we
further investigate the relation between tunneling rate and the spectral shape
of ${\cal T}_{LR}(\varepsilon)$. In Fig. 6, we show a comparison of the
tunneling spectra, ${\cal T}_{LR}(\varepsilon)$ of a t-ZGNR with 15 GQDs
calculated for two different tunneling strengths $\Gamma_{t}=0.54eV$ and
$\Gamma_{t}=t_{pp\pi}=2.7eV$. In Fig. 6(a) the area below the ${\cal
T}_{LR}(\varepsilon)$ curve shows a right-triangle shape that has a steep
change with respect to $\varepsilon$ on the side toward the central gap. In
Fig. 6(b) (with $\Gamma_{t}=t_{pp\pi}=2.7eV$) the electron tunneling
probability through the electronic states near $\varepsilon_{LU(HO)}$ is
highly suppressed, leading to an arch-like area under the ${\cal
T}_{LR}(\varepsilon)$ curve. The corresponding maximum $ZT$ for this case is
$~{}2.3$, which is much smaller than that of Fig. 5(d). The results of Fig.
6(a) and Fig. 6(b) indicate that the shape of the area under ${\cal
T}_{LR}(\varepsilon)$ curve depends on the tunneling rate (or the quality of
contact). An arch-like area under ${\cal T}_{LR}(\varepsilon)$ will reduce the
electrical conductance $G_{e}$ near the central gap.

As seen in Fig. 6(c) the optimized $PF$ found near $\Gamma_{t}=0.54eV$ is very
close to that obtained by using an ideal SF transmission spectrum (indicated
by the red line), which exhibits the quantum limit of power factor for 1D
systems with $PF_{QB}=1.2659(\frac{2k^{2}_{B}}{h})$ [Whitney, 9,ChenIJ, 11].
We obtain the optimized $PF_{max}=0.9PF_{QB}$, which is the same as that of
the quantum interference heat engine [Samuelsson, 62]. This analysis suggests
that using SGQDs with suitable tunneling rates could achieve the performance
close to an optimum heat engine with maximum electrical power output and high
thermoelectric efficiency. With the same $\kappa_{ph}$, $ZT_{max}$ of t-ZGNRs
can reach $95\%$ of that obtained by a SF transmission coefficient (with
$ZT_{max}=3.951$). So far, we have only considered periodically indented
structures on both the top and bottom zigzag edges (see Fig. 1(b)). We note
that the same optimized result of $PF_{max}=0.9PF_{QB}$ can also be achieved
in periodically indented structures corresponding to Fig. 6 but with textured
pattern only on one zigzag side provided that the tunneling rate can be
increased to $\Gamma_{t}=0.72eV$. The results are shown in Fig. A.3 in
appendix.

Figure 6: (a) and (b) Transmission coefficients ${\cal T}_{LR}(\varepsilon)$
as functions of energy for a t-ZGNR with $\Gamma_{t}=0.54eV$ and
$\Gamma_{t}=2.7eV$. The ${\cal T}_{LR}(\varepsilon)$ with square form (SF) is
also plotted by a curve with red color. (c) and (d) Power factor as functions
of $\mu$ for two tunneling rates at $k_{B}T=27$ meV. Other physical parameters
are the same as those used in Fig. 5.

###  III.3 Armchair graphene nanoribbons

As shown in Figs. 5-6, SGQDs formed by metallic ZGNRs can become
semiconductors. Next, we investigate whether its thermoelectric performance is
better than that of AGNRs or t-AGNRs. We note that many designs have focused
on the optimization of thermoelectric performance of AGNRs [Mazzamuto,
63-Merino, 66]. Here, we first consider AGNRs with their zigzag ends coupled
to electrodes as depicted in Fig. 1(c). Figure 7 shows the calculated
transmission coefficient ${\cal T}_{LR}(\epsilon)$ for different defect
locations with $N_{z}=7$, $N_{a}=100$ and $\Delta=5.4$eV. We choose
$\Gamma_{t}=2.7$eV, which gives the optimized shape of the transmission
coefficient in defect free situation in Fig. 7(a). For finite-size AGNRs the
conduction (valence) subband states are quantized, leading to closely-spaced
peaks with staircase-like structures as revealed by the ${\cal
T}_{LR}(\varepsilon)$ spectrum. The area under ${\cal T}_{LR}(\varepsilon)$
curve for states derived from the first subband has a parabolic shape, which
does not meet the criterium for achieving optimized thermoelectric property.
For a single defect at edge location (1,3), all but the transmission
coefficient spectrum of the first conduction subband are affected
significantly. When the defect is located at site (3,3), the spectra for the
first conduction subband and the first valence subband are unaffected .
Nevertheless, when the defect occurs at (4,4) the spectra in these two
subbands show a remarkable change. In current studies we do not observe the
Anderson localization effect, which leads to vanishing $G_{e}$ in the subband
regions when a single defect occurs on the edges of ZGNRs and AGNRs [LeePA,
67]. This implies that GNRs considered are not exact 1D systems.

Figure 7: Transmission coefficient of AGNRs as functions of $\varepsilon$ for
different defect locations at$\Gamma_{t}=2.7eV$, $N_{z}=7$ and $N_{a}=100$
($L_{a}=10.5$nm). We have adopted $\Delta=5.4$eV.

To illustrate the effects of defects on thermoelectric properties of finite
AGNRs, we show the calculated electrical conductance ($G_{e}$), Seebeck
coefficient ($S$), power factor ($PF$), and figure of merit ($ZT$) as
functions of $\mu$ at $k_{B}T=27$ meV in Fig. 8. Four curves in each diagram
correspond to the effects of defects located at four sites considered in Fig.
7. As seen in Fig. 8a, the spectra of $G_{e}$ for defects at locations (1,3)
and (4,4) are degraded seriously either in the valence subband or conduction
subband. Due to the robustness of the Seebeck coefficient in the central-gap
region, the effect of defects on the power factor is solely determined by
$G_{e}$. Although the maximum $S$ is larger than that of Fig. 5(c) due to the
larger gap in AGNRs with $N_{z}=7$, the $PF$ values are smaller than the
corresponding values of the optimized t-ZGNR shown in Fig. 6 (with
$\Gamma_{t}=0.54$eV). This is mainly attributed to the different shapes in the
area under the transmission-coefficient curve. As a consequence, the $G_{e}$
resulting from thermionic tunneling effect is less favorable for AGNRs. Note
that we have adopted $\kappa_{ph}=0.0285nW/K$ in Fig. 8(d), which is the same
as $\kappa_{ph}$ in Fig. 5(d).

Figure 8: (a) Electrical conductance $G_{e}$, (b) Seebeck coefficient $S$, (c)
power factor $PF$ and (d) figure of merit $ZT$ as functions of $\mu$ for
different defect locations at $k_{B}T=27meV$. Other physical parameters are
the same as those of Fig. 7. The considered $\kappa_{ph}=0.0285nW/K$ is the
same as that of Fig. 5 ($\kappa_{ph}=0.0285nW/K$).

###  III.4  SGQDs formed by textured AGNRs

Figure 9: Transmission coefficient of SGQD with $N_{z}=5$ and various $N_{a}$
values as functions of energy $\varepsilon$. $\Gamma_{t}=2.7eV$. Here,
$N_{a}=12+8\times N_{QD}$, where $N_{QD}$ is the number of GQDs in the
interior region of the structure, and they are sandwiched between two smaller
GQDs of $N_{a}=6$ connected to electrodes. Thus, $N_{a}=100,116$, and 132
correspond to $N_{QD}=11,13$, and 15, respectively. For $N_{a}=116$ and
$N_{a}=132$, the lengths of SGQDs are $12.2$ nm and $13.9$ nm, respectively.

Finally, we study the thermoelectric effect of SGQD structures based on the
t-AGNR structure as shown in Fig. 1d. The electronic structures and density of
states (DOS) of t-AGNR superlattices have been studied theoretically by using
density functional theory (DFT) [Sevincli, 40,Rizzo, 21]. Electronic
structures of t-AGNRs can also be calculated by using a tight binding model
and results are in good agreement with the DFT calculation[Sevincli, 40].
Figure 9 shows the calculated ${\cal T}_{LR}(\varepsilon)$ of SGQD junctions
made from t-AGNRs as functions of $\varepsilon$ for three different $N_{a}$
values with $N_{z}$ fixed at 5. We adopt the optimized tunneling rate with
$\Gamma_{t}=2.7eV$. Although the AGNR with $N_{z}=5$ is metallic, the textured
AGNR can be semiconducting as illustrated by the sizable central gap
($~{}1.098eV$) in Fig. 9. In Fig. A.2(b) of the appendix, we show a comparison
of subband structures of the t-AGNR and unaltered AGNR (with an enlarged super
cell of length $L=2a^{\prime}$ to match the unit cell of t-AGNR), where
$a^{\prime}=\sqrt{3}a$ is the unit cell length of AGNR. We see that the two
subbands with linear dispersion in the metallic $N_{z}=5$ AGNR are split into
two subbands with parabolic dispersion near zero wave-number once the quantum
constriction takes effect in the t-AGNR. It is worth noting that edge states
with zero energy can exist at the left and right ends (with zigzag edges) of
the truncated t-AGNR in contact with electrodes. These edge states can
contribute to ${\cal T}_{LR}(\varepsilon)$ for t-AGNRs with short SGQD
structure (see Fig. A.4 in appendix), since their wave functions decay
exponentially along the armchair direction. The area under the ${\cal
T}_{LR}(\varepsilon)$ curve has an arch-like shape that does not change much
as we vary the tunneling rate.

Figure 10 shows the calculated $G_{e}$, $S$, $PF$ and $ZT$ of SGQD junctions
made from a t-AGNR with $N_{z}=5$ and $N_{a}=132$ ($N_{QD}=15$) as functions
of temperature for different $\mu$ values. Solid and dashed lines correspond
to $\mu=0.52$eV and $\mu=0.36$eV, respectively. The first conduction subband
edge occurs at $\mu_{edge}=0.544$eV. At a given temperature $k_{B}T=27$meV,
the maximum power factor occurs at $\mu=0.52$eV. For $\mu=0.36$eV, which is
far away from $\mu_{edge}$, the calculated $G_{e}$ in the thermionic assisted
tunneling process (TATP) is extremely small, whereas its Seebeck coefficient
is highly enhanced. The temperature-dependent $S$ for this case is
complicated. Its temperature dependence can be described by three different
functions at three temperature ranges as illustrated in Fig. 10(b). In region
one ($T<90K$), we have $S_{1}=-\frac{\pi^{2}k^{2}_{B}T}{3e}\frac{\partial
ln({\cal T}_{LR}(\varepsilon))}{\partial\varepsilon}|_{\varepsilon=\mu}$. In
region three ($T>120K$) we have $S_{3}=\frac{\mu-\mu_{edge}}{eT}$ [Kuo1, 68].
In region 2, the analytic expression of $S_{2}$ is unknown. In TATP, $G_{e}$
can be described by the expression $exp^{(\mu-\mu_{edge})/(k_{B}T)}$. As a
consequence, the power factor and figure of merit are small for $\mu=0.36$eV.
When $\mu=0.52$eV, the $PF$ is enhanced quickly in the range $50K<T<120K$.
Within this temperature region, $ZT$ shows a similar behavior of $PF$. When
$T>120K$, the temperature-dependent $ZT$ has the same trend of $S$ because the
heat current is dominated by the linear-$T$ phonon thermal conductance
$\kappa_{ph}=F_{s}\frac{\pi^{2}k^{2}_{B}T}{3h}$, which cancels the factor of
$TG_{e}$ in the numerator of $ZT$. It is worth noting that the optimized $PF$
of t-AGNRs at $T=325K$ is only one-half of the power factor of the t-ZGNR
shown in Fig. 6(c).

Overall, we found no appreciable improvement in the power factor of t-AGNRs in
comparison to AGNRs. Although electron Coulomb interactions are neglected in
this calculation, our conclusion is still valid even in the Coulomb blockade
regime. When the thermoelectric behavior of SGQDs is dominated by the TATP,
electron-electron correlation functions are usually small [Kuo3, 69].
Therefore, we can neglect electron Coulomb interactions when $\mu$ is inside
the gap between two subbands.

Figure 10: (a) Electrical conductance $G_{e}$, (b) Seebeck coefficient $S$,
(c) power factor $PF$ and (d) figure of merit $ZT$ of SGQD with $N_{z}=5$ and
$N_{a}=132$ ($N_{QD}=15$) as functions of temperature for different $\mu$ at
$\Gamma_{t}=2.7eV$. The phonon thermal conductance, $\kappa_{ph}$ adopted is
the same as in Fig. 8.

##  IV Conclusion

We have theoretically investigated the transport and thermoelectric properties
of SGQDs, which are formed by tailoring ZGNRs and AGNRs. Our calculations are
based on the Green’s function approach within a tight-binding model. Electron
coherent tunneling process is found responsible for the electrical conductance
spectra of SGQDs. Subband width and central gap of SGQDs can be modulated by
varying the size of GQDs and inter-dot coupling strength. Unlike $G_{e}$ and
$\kappa_{e}$, Seebeck coefficients are found insensitive to the contact
property between the SGQDs and electrodes. As a result, the power factor and
thermoelectric figure of merit can be improved by modulating the tunneling
rates. The maximum $ZT$ values at room temperature occur when the chemical
potential, $\mu$ is close to the HOMO (LUMO) level in typical situation where
phonon carriers dominate the heat transport. As shown in Fig. 6 and Fig. A.3,
the maximum power factor (figure of merit) of t-ZGNRs at room temperature can
reach 90 $\%$ (95 $\%$) of the ideal situation with a square-form transmission
coefficient. We also found that SGQDs based on t-ZGNRs can outperform SGQDs
based on t-AGNRs for thermoelectric application. The significantly improved
thermoelectric behavior of the textured ZGNR is attributed to the sharp change
of its transmission coefficient near the central gap. We found that defects at
interior sites will not ruin the robust behavior of ${\cal
T}_{LR}(\varepsilon)$ associated with the edge states near LUMO of t-ZGNRs. It
implies that the electron mean free path,$\lambda_{e}$ contributed from the
zigzag edge states ($\lambda_{e,edge}$) can be much larger than the
contribution from bulk states ($\lambda_{e,bulk}$)). When the channel length
($L_{z}$) of t-ZGNRs is much larger than $\lambda_{ph}$ and
$\lambda_{e,bulk}$, $\kappa_{ph}$ is seriously suppressed in such a diffusing
region. The electronic states in the first miniband (near LUMO) of t-ZGNRs
showing nonlinear dispersion could remain coherent due to their unique nature
[Chou, 27]. As a consequence, the power factor of Fig. 6 remains valid,
meanwhile the corresponding $ZT$ could be further enhanced.[Xu, 13] At room
temperature with $T=324K(k_{B}T=27meV)$, the electrical power output can reach
$0.212nW/K$ for each SGQD implemented by using t-ZGNR. For an SGQD array with
density of $5\times 10^{6}cm^{-1}$ and $ZT$ larger than 3, the estimated power
output is around $mW/K$, which can be applicable for low-power wearable
electronic devices[SuarezF, 70].

Author contributions
David M T Kuo initiated the idea and performed the calculation.; David M T Kuo
and Y. C. Chang both contributed to the development of the computation code
and carried out the data analyses.David M. T. Kuo and Y. C. Chang wrote the
manuscript. All authors read and agreed to the published version of the
manuscript

Funding
This research was funded by the Ministry of Science and Technology Grant
Number MOST 110-2112-M-001-042 and 110-2119-M-008-006-MBK in Taiwan

Data availability
The data presented in this study are available upon reasonable request.

Acknowledgments
This work was supported in part by the Ministry of Science and Technology
(MOST), Taiwan under Contract Nos. 110-2112-M-001-042 and
110-2119-M-008-006-MBK.

Conflicts of interest
The authors declare no conflict of interest.

## Appendix A

###  A.1 Charge density of t-ZGNRs

To further clarify the electronic states shown in Fig. 3, we plot the charge
densities of $|\Psi_{\ell,j}(\varepsilon_{LU})|^{2}$,
$|\Psi_{\ell,j}(\varepsilon_{e,1})|^{2}$ and
$|\Psi_{\ell,j}(\varepsilon_{e,2})|^{2}$ for the 6-GQDs case formed by t-ZGNRs
with $N_{a}=12$ and $N_{z}=47$ in Fig. A.1, which correspond to the electronic
states of $\varepsilon_{LU}=0.247eV$, $\varepsilon_{e,1}=0.3105eV$ and
$\varepsilon_{e,2}=0.468eV$ in Fig. 3(f). As seen in Fig. A.1(a), the charge
density distribution in each GQD is the same for the electronic state of
$\varepsilon_{LU}$. The maximum charge densities appearing at bottom and top
zigzag edge sites are contributed by $\ell=1$ and $\ell=12$ with the condition
$j-4=8m$, where $m$ is an integer. At interior sites of each GQD the charge
density is very low. For $\varepsilon_{e,1}=0.3105eV$, the charge density at
the central two GQDs is negligible. For $\varepsilon_{e,2}=0.468eV$, the
charge density at the $2nd$ and $5th$ GQDs are dilute. For
$\varepsilon_{e,6}=1.0665eV$, the charge density at the $1st$ and the $6th$
GQDs is smaller than those of other GQDs (not shown here). This indicates that
the electronic states near $\varepsilon_{e,6}=1.0665eV$ are weakly coupled to
the electrodes. Based on the charge density distribution shown in Fig. A.1,
the defect location-dependent behavior of $G_{e}$ shown in Fig. 4 can be
clearly illustrated. The $G_{e}$ does not be degraded seriously when defects
locate at sites with low charge density.

Figure A.1: Charge density $|\Psi_{\ell,j}(\varepsilon)|^{2}$ for the t-ZGNR
with conductance spectrum shown in Fig. 3(f). Diagrams (a), (b) and (c) are
the charge densities of $\varepsilon_{LU}=0.247eV$,
$\varepsilon_{e,1}=0.3105eV$, and $\varepsilon_{e,2}=0.468eV$, respectively.

###  A.2 Electronic band structures

In Fig. A.2, we show comparisons of subband structures of graphene nanoribbons
with and without modification of side edges. Here, the length of unit cell in
t-GNRs is denoted $L$. For the t-ZGNR with its unit cell characterized by
$N_{a}=8$ and $N_{z}=8$, we have $L=4a$, where $a$ is lattice constant of
ZGNR. For the t-AGNR characterized by $N_{z}=5$ and $N_{a}=8$, we have
$L=2a^{\prime}$, where $a^{\prime}=\sqrt{3}a$ is the lattice constant of AGNR.
We note that the first Brillouin zone (BZ) has one quarter (half) the size of
the BZ of ZGNR (AGNR). Thus, the subband structures shown in Fig. A.1 includes
the zone-folding effect that maps the zone boundary $k_{z}=\pi/a$
($k_{z}=\pi/a^{\prime}$) of the unfolded BZ of ZGNR (ANGR) to $k_{z}=0$ of the
folded BZ here. Here, we show the electronic structure of $N_{a}=8$ rather
than $N_{a}=12$ in Fig. A.2(a). For $N_{a}=8$, we observe not only the central
gap, but also the gap between the first subband and the second subband.

Figure A.2: (a) Subband structures of t-ZGNR superlattice with a unit cell of
length $L$ characterized by $N_{a}=8$ and $N_{z}=8$ (solid lines) as depicted
in Fig. 1b, and corresponding ZGNR (without texture) of the same length of
unit cell (dashed lines). (b) Subband structures of t-AGNR superlattice with a
unit cell characterized by $N_{z}=5$ and $N_{a}=8$ (solid lines) as depicted
in Fig. 1(d), and corresponding AGNR (without texture) of the same length of
unit cell (dashed lines).

###  A.3 ZT optimization of textured ZGNRs

In Fig. A3, we show the calculated ${\cal T}_{LR}(\varepsilon)$ in the region
of $\varepsilon>0$ for two kinds of t-ZGNRs: (a) with indentation on both
zigzag sides and (b) with indentation on only one zigzag side of the t-ZGNR.
Both t-ZGNRs consist of 15 GQDs (same as the one considered in Figs. 5 and 6).
The solid, dashed, dotted and dot-dashed lines correspond to
$\Gamma_{t}=0.45,0.54,0.63$ and $0.72$eV, respectively. As energy
$\varepsilon$ is tuned away from LUMO, many sharp peaks appear in the spectrum
as a consequence of the Fabry-Perot interference. The period of Fabry-Perot
interference in (b) is larger than that in (a). When $\varepsilon$ gets close
to $LUMO$, electronic states can not be clearly resolved due to their
nonlinear dispersion. Diagram (c) shows $ZT_{max}$ of a t-ZGNR consisting of
15 GQDs as a function of the tunneling rate $\Gamma_{t}$ at $k_{B}T=27$meV and
$\mu=0.423$eV. It is found that $ZT_{max}$ can be larger than 3 for
$\Gamma_{t}$ between 0.18 eV and 1.45 eV for the case of (a). When the barrier
hight gets smaller, the miniband width becomes wider. This indicates that
electron interdot hopping strengths are enhanced. That is why $ZT_{max}$
occurs at $\Gamma_{t}=0.72$ eV in (b). Since the theoretical limit for an
ideal SF transmission coefficient is $ZT_{max}=3.95167$ , the maximum $ZT$
value of t-ZGNRs can reach $95\%$ of the ideal case.

Figure A.3: Transmission coefficient in the range of $\varepsilon>0$ for
different tunneling rates. (a) periodically indented structures on both
zigzag-edge sides and (b) periodically indented structures on only one zigzag-
edge side. (c) $ZT_{max}$ of a t-ZGNR consisting of 15 GQDs as a function of
the tunneling rate $\Gamma_{t}$. Each GQD in the structure has size $N_{a}=8$
and $N_{z}=7$. The maximum $ZT$ value of SF transmission coefficient at
$k_{B}T=27$meV is 3.95167 when its $\kappa_{ph}$ is adopted the same as that
of Fig. 5.

###  A.4 Local density of states of textured AGNRs

The ${\cal T}_{LR}(\varepsilon)$ spectrum cannot reveal the existence of
zigzag edge states in a finite-size AGNR that decay exponentially along the
armchair directions when the channel length is much longer than the decay
lengths. These localized states of GNRs can be probed by STM [Rizzo, 21]. In
Fig. A.4, we show the calculated local density of states (LDOS) as a function
of $\varepsilon$ for various tunneling rates. Here, we define
$LDOS=-Im(G^{r}_{\ell=2,j=1}(\varepsilon))/\pi$. The LDOS can reveal the
existence of the zigzag edge states of an AGNR coupled to the electrodes. The
zero energy modes resulting from the outer zigzag edge structures of t-AGNRs
are seen in Fig. A.4. Due to the localized nature of the zigzag edge states,
the spectra in Fig. 9 can not reveal them. Meanwhile, these edge state
localized at the interface are also strongly influenced by the variation of
coupling strength between the electrodes and the t-AGNR. With increasing
$\Gamma_{t}$, these zero energy modes become broadened in (b) and (c).

Figure A.4: Local density of states of a t-AGNR with $(N_{z},N_{a})=(5,132)$
as a function of $\varepsilon$ for various tunneling rates at site (2,1).

## References

  * (1) G. Chen, M. S. Dresselhaus, G. Dresselhaus, J. P. Fleurial, and T. Caillat, International Materials Reviews, 48, 45 (2003).
  * (2) L. D. Hicks and M. S. Dresselhaus, Phys. Rev. B 47, 16631 (1993).
  * (3) L. D. Hicks and M. S. Dresselhaus, Phys. Rev. B 47, 12727 (1993).
  * (4) A. Boukai, Y. Bunimovich, J. Tahir-Kheli, J.-K. Yu, W. A. Goddard III, and J. R. Heath, Nature 451  168 (2008).
  * (5) A. Hochbaum, A. I. Hochbaum, R. Chen, R. D. Delgado, W. Liang, E. C. Garmett, M. Najarian, A. Majumdar, and P. D. Yang, Nature 451 163 (2008).
  * (6) J. P. Heremans, C. M. Thrush, D. T. Morelli, and M. C. Wu, Phys. Rev. Lett. 88 216801 (2002).
  * (7) D. L. Nika, E. P. Pokatilov, A. A. Balandin, V. M. Fomin, A. Rastelli, and O. G. Schmidt, Phys. Rev. B 84 165415 (2011).
  * (8) M. Hu and D. Poulikakos, Nano Lett.12 5487 (2012).
  * (9) R. S. Whitney, Phys. Rev. Lett. 112, 130601 (2014).
  * (10) P. Pietzonka and U. Seifert, Phys. Rev. Lett. 190602 (2018).
  * (11) I. J Chen, A. Burke, A. Svilans, H. Linke, and C. Thelander, Phys. Rev. Lett. 120, 177703 (2018).
  * (12) R. Luo, G. Benenti, G. Casati, and J. Wang, Phys. Rev. Lett. 121, 080602 (2018).
  * (13) Y. Xu, Z. Gan, and S. C. Zhang, Phys. Rev. Lett. 112, 226801 (2014).
  * (14) P. Darancet, V.Olevano, and D. Mayou, Phys. Rev. Lett. 102, 136803 (2009).
  * (15) C. R. Kagan and C. B. Murry, Nature Nanotechnology 10, 1013 (2015).
  * (16) W. I. L. Lawrie et al, Appl. Phys. Lett. 116, 080501 (2020).
  * (17) K. S. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang, Y. Zhang, S. V. Dubonos, I. V. Grigorieva, and A. A. Firsov, Science 306 666 (2004).
  * (18) J. Cai, P. Ruffieux, R. Jaafar, M. Bieri, T. Braun, S. Blankenburg, M. Muoth, A. P. Seitsonen, M. Saleh, X. Feng, K. Mullen, and R. Fasel, Nature 466, 470 (2010).
  * (19) Y. C. Chen, T. Cao, C. Chen, Z. Pedramraz, D. Haberer, D. G. de Oteyza, F. R. Fischer, S. G. Louie and M. F. Crommie, Nature Nanotechnology 10, 156 (2015).
  * (20) O. Groning, S. Wang, X. Yao, C. A. Pignedoli, G. B. Barin, C. Daniels, A. Cupo, V. Meunier, X. Feng, A. Narita, et al., Nature 560, 209 (2018).
  * (21) D. J. Rizzo, G. Veber, T. Cao, C. Bronner, T. Chen, F. Zhao, H. Rodriguez, S. G. Louie, M. F. Crommie, and F. R. Fischer, Nature 560, 204 (2018).
  * (22) L. H. Yan and P. Liljeroth, ADVANCES IN PHYSICS: X 4, 1651672 (2019).
  * (23) D. J. Rizzo, G. Veber3, J. W. Jiang, R. McCurdy, T. Cao C. Bronner, T. Chen, Steven G. Louie1, F. R. Fischer, M. F. Crommie, Science 369, 1597 (2020).
  * (24) Q. Sun, Y. Yan, X. L. Yao, K. Mullen, A. Narita, R. Fasel, and P. Ruffieux, J. Phys. Chem. Lett. 12, 8679 (2021).
  * (25) D. J. Rizzo, J. W. Jiang, D. Joshi, G. Veber, C. Bronner, R. A. Durr, P. H. Jacobse, T. Cao, A. Kalayjian, H. Rodriguez, P. Butler, T. Chen, Steven G. Louie, F. R. Fischer, and M. F. Crommie, ACS Nano 15, 20633 (2021).
  * (26) J. P. Llinas et al, Nature communications  8, 633 (2017).
  * (27) K. S. Lin and M. Y. Chou, Nano Lett. 18, 7254 (2018).
  * (28) Y. M. Zuev, W. Chang, and P. Kim, Phys. Rev. Lett. 102, 096807 (2009).
  * (29) P. Wei, W. Z. Bao, Y. Pu, C. N. Lau, and J. Shi, Phys. Rev. Lett. 102, 166808 (2009).
  * (30) Y. Xu, Z.Y. Li, and W. H. Duan, Small 10, 2182 (2014).
  * (31) Y. H. Wang, J. Yang, L. Y. Wang, K. Du, Q. Yin, and Q. J. Yin, ACS Appl. Mater. Interfaces 9, 20124 (2017).
  * (32) M. Li, D. L. Cortie, J. X. Liu, D. H. Yu, SMKN Islam, L. L. Zhao, DRG Mitchell, R. A. Mole, M. B. Cortie, and S. X. Dou, Nano energy 53, 993 (2018).
  * (33) S. Ghosh, S. Harish, M. Ohtaki, and B. B. Saha, Material today energy,  18, 100492 (2020).
  * (34) Y.Y. Wang, D. R. Chen,J. K. Wu, T. H. Wang, C. S. Chuang , S. Y. Huang, W. P. Hsieh, M. Hofmann, Y. H. Chang, and Y. P. Hsieh, Nano Lett. 21, 6990 (2021).
  * (35) Jian He and Terry M. Tritt, Science 357, 1369 (2017).
  * (36) K. Nakada, M. Fujita, G. Dresselhaus ,and M. S. Dresselhaus Phys. Rev. B 54, 17954 (1996).
  * (37) K. Wakabayashi, M. Fujita, H. Ajiki, and M. Sigrist, Phys. Rev. B 59, 8271 (1999).
  * (38) Y. W. Son, M. L. Cohen, and Steven G. Louie, Nature 444, 347 (2006).
  * (39) M. Topsakal, H. Sevincli, and S. Ciraci, Appl. Phys. Lett. 92, 173118 (2008).
  * (40) H. Sevincli, M. Topsakal, and S. Ciraci, Phys. Rev. B 78, 245402 (2008).
  * (41) H. Haug and A. P. Jauho, Quantum Kinetics in Transport and Optics of Semiconductors (Springer, Heidelberg, 1996).
  * (42) D. M. T. Kuo, AIP Advances 10, 045222 (2020).
  * (43) T. T. Phung, R. Peters, A. Honecker, G. Trambly de Laissardiere, and J. Vahedi, Phys. Rev. B 102, 035160 (2020).
  * (44) R. S. Chen, G. L. Ding, Y. Zhou, and S. T. Han, J. Materials Chemistry C, 9 11407 (2021).
  * (45) G. D. Mahan and L. M. Woods, Phys. Rev. Lett. 80, 4016 (1998).
  * (46) Y. C. Chang, J. N. Schulman, G. Bastard, and Y. Guldner, Phys. Rev. B 31, 2557 (R)(1985).
  * (47) B. A. Bernevig, T. L. Hughes, and S. C. Zhang, Science 314, 1757 (2006)
  * (48) L. Z. Kou, S. C. Wu, C. Felser, T. Frauenheim, C. F. Chen, and B. H. Yan, ACS NANO 8, 10448-10454 (2014).
  * (49) H. Q. Huang and W. H. Duan, Nature Materials 15, 129 (2016).
  * (50) Z. L. Zhu, X. L. Cai, S. H. Yi, L. Chen, Y. W. Dai, C. Y. Niu, Z. X. Guo, M. H. Xie,F. Liu, J.H. Cho, Y. Jia, and Z. Y. Zhang, Phys. Rev. Lett. 119, 106101 (2017).
  * (51) L. Xian, A. PerezPaz, E. Bianco, P. M. Ajayan and A. Rubio, 2D materials. 4, 041003 (2017).
  * (52) G. D. Mahan and J. O. Sofo, Proc. Natl. Acad. Sci. USA 93, 7436 (1996).
  * (53) P. H. Chang, M. S. Bahramy, N. Nagaosa, and B. K. Nikolic, Nano Lett. 14, 3779 (2014).
  * (54) K. Wakabayashi, K Sasaki, T. Nakanishi, and T. Enoki, Sci. Technol. Adv. Mater. 11 054504 (2010).
  * (55) D. A. Areshkin, Daniel Gunlycke, and C. T. White, Nano Lett.7, 204 (2007).
  * (56) T. B. Martins, Antonio J. R. da Silva, R. H. Miwa and A. Fazzio, Nano Lett.8, 2293 (2008).
  * (57) T. C. Li and Shao-Ping Lu, Phys. Rev. B 77, 085408 (2008).
  * (58) P. C. Shen et al, Nature 593, 212 (2021).
  * (59) H. Zheng, H. J. Liu, X. J. Tan, H. Y. Lv, L. Pan, J. Shi, and X. F. Tang, Appl. Phys. Lett. 100, 093104 (2012).
  * (60) H. Sevincli and G. Cuniberti, Phys. Rev. B 81, 113401 (2010).
  * (61) P. G. Murphy and J. E. Moore, Phys. Rev. B 76, 155313 (2007).
  * (62) P. Samuelsson, S. Kheradsoud, and B. Sothmann, Phys. Rev. Lett. 118, 256801 (2017).
  * (63) F. Mazzamuto, V. Hung Nguyen, Y. Apertet, C. Caer, C. Chassat, J. Saint-Martin, and P. Dollfus, Phys. Rev. B 83, 235426 (2011).
  * (64) P. Dollfus, V. H. Nguyen, and J. Saint-Martin, J. Phys: Condens. Matter 27, 133204 (2015).
  * (65) V. T. Tran, J. Saint-Martin, and P. Dollfus, Nanotechnolgy 26, 495202 (2015).
  * (66) N. Merino-Diez, A. Garcia-Lekue, E. Carbonell-Sanroma, J. C. Li, M. Corso, L. Colazzo, F. Sedona, D. Sanchez-Portal, J. I. Pascual, and D. G de Oteyza, ACS NANO, 11, 11661 (2017).
  * (67) P. A. Lee and T. V. Ramakrishnan, Rev. Mod. Phys. 57, 287 (1985).
  * (68) D. M.-T. Kuo and Y. C. Chang, Phys. Rev. B 81, 205321 (2010).
  * (69) D. M. T. Kuo, C. C. Chen, and Y. C. Chang, Phys. Rev. B 95, 075432 (2017).
  * (70) F. Suarez, A. Nozariasbmarz, D. Vashaee, and M. C. Ozturk, Energy Environ Sci, 9, 2099 (2016).