# A DEMATEL-Based Completion Method for Incomplete Pairwise Comparison Matrix
in AHP

Xinyi Zhou  Yong Hu  Yong Deng<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>Felix T.S. Chan  Alessio Ishizaka  School of Computer
and Information Science, Southwest University, Chongqing, 400715, China.  Big
Data Decision Institute, Jinan University, Tianhe, Guangzhou, 510632, China.
Institute of Integrated Automation, School of Electronic and Information
Engineering, Xi’an Jiaotong University, Xi’an, Shaanxi, 710049, China.  School
of Engineering, Vanderbilt University, Nashville, TN, 37235, USA.  Department
of Industrial and Systems Engineering, The Hong Kong Polytechnic University,
Hung Hom, Hong Kong  University of Portsmouth, Portsmouth Business School,
Richmond Building, Portland Street, Portsmouth PO1 3DE, United Kingdom

###### Abstract

Pairwise comparison matrix as a crucial component of AHP, presents the
preference relations among alternatives. However, in many cases, the pairwise
comparison matrix is difficult to complete, which obstructs the subsequent
operations of the classical AHP. In this paper, based on DEMATEL which has
ability to derive the total relation matrix from direct relation matrix, a new
completion method for incomplete pairwise comparison matrix is proposed. The
proposed method provides a new perspective to estimate the missing values with
explicit physical meaning. Besides, the proposed method has low computational
cost. This promising method has a wide application in multi-criteria decision-
making.

###### keywords:

AHP , Incomplete pairwise comparison matrix , Completion method , DEMATEL

††journal: Annals of Operations Research

##  1 Introduction

Analytic hierarchy process (AHP) is a multi-criteria decision-making method
that helps the decision-makers facing a complex problem with multiple
conflicting and subjective criteria[1]. AHP has been widely applied in various
areas[2, 3, 4, 5]. Based on a hierarchical structure, the criteria weights can
be derived from the pairwise comparison matrix. Pairwise comparison matrix as
a crucial component of AHP, is commonly utilized to estimate preference values
of finite alternatives with respect to a given criterion.

However, the pairwise comparison matrix is difficult to complete in the
following cases, which obstructs the subsequent operations of the classical
AHP.

  * 1.

The experts lack the knowledge of one or more alternatives.

  * 2.

Partial data in pairwise comparison matrix has been lost.

  * 3.

Huge number of pairwise comparison requested: $(n^{2}-n)/2$ for n
alternatives. For example, 8 alternatives and 6 criteria require 183
entries[6].

Hitherto, it has attracted widespread attention over the issue of incomplete
pairwise comparison matrix in AHP. There are quantities of methods having been
developed to address the incomplete comparison matrices. Some of them use the
incomplete matrix to assign priority weight to each alternative without
completion[7, 8, 9, 10, 11]. Some methods are devoted to estimate the missing
values in matrix[12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22]. While these
methods have a good effect on allocating the priority weights to alternatives,
the methods that focus on matrix completion are more “well-rounded” to restore
the preference relations between each pair of alternatives, which are
inaccessible before completion.

After reviewing the existing completion methods for incomplete pairwise
comparison matrix in AHP, it can be found that the vast majority of these
methods focus on the optimal consistency of pairwise comparison matrix[12, 13,
14, 15, 16]. For example, Fedrizzi and Giove[12] proposed a new method to
complete the pairwise comparison matrix by minimizing its global
inconsistency. Ergu _et al._[13] extended the geometric mean induced bias
matrix to estimate the missing judgments and improve the consistency ratios by
the least absolute error method and the least square method. Based on the
graph theory, Bozoíki _et al._[14] developed a completion method which focuses
on the optimal consistency and estimates the missing values through iteration.
Besides, Hu and Tsai[17] proposed a well-known back propagation multi-layer
perception to estimate the missing comparisons of incomplete pairwise
comparison matrices in AHP. Alonso _et al._[18] defined a way to calculate the
consistency of incomplete pairwise comparison matrix and maintain its
consistency level after completion. Benítez _et al._[19] used the Frobenius
norm to obtain the most consistent/similar complete matrix with the initial
incomplete matrix. Gomez _et al._[20] proposed a completion method based on
the neural network. These methods are useful but have some limitations:

  * 1.

Most methods complete the pairwise comparison matrix with optimal consistency,
which sometimes can not restore the incomplete matrix to its original state,
especially when the initial matrix is not high consistent.

  * 2.

The methods are computationally expensive.

Taking into account these issues, a DEMATEL-based completion method for
incomplete pairwise comparison matrix in AHP is proposed. Based on DEMATEL
method which has ability to derive the total relation matrix from direct
relation matrix, the completion of incomplete pairwise comparison matrix can
be divided into three steps.

  * -

Firstly, turn the incomplete comparison matrix into direct relation matrix.

  * -

Secondly, DEMATEL method is adopted to obtain the total relation matrix.

  * -

Thirdly, transform the total relation matrix into pairwise comparison matrix
with reciprocal preference relations.

The proposed method provides a new perspective to estimate the missing values
in pairwise comparison matrix with explicit physical meaning (i.e. by
calculating the total relations among alternatives). Besides, the proposed
method has low computational cost.

The rest of this paper is organized as follows. Section 2 introduces some
preliminaries of this work. Section 3 illustrates the procedure of the
proposed method and an application for ranking tennis players is presented.
Section 4 discusses the superiority and efficiency of the proposed method and
section 5 ends the paper with the conclusion.

##  2 Preliminaries

###  2.1 Analytic Hierarchy Process (AHP)

AHP is a multi-criteria decision-making method developed by Saaty[23] and aims
at quantifying relative weights for a given set of criteria on a ratio
scale[24]. As a decision-making approach, it permits a hierarchical structure
of criteria, which provides users with a better focus on specific criteria and
sub-criteria when allocating the weights[25]. Besides, AHP has ability to
judge the consistency of pairwise comparison matrix to show its potential
conflicts in decision-making process. AHP has been widely used in supply chain
management[26, 27, 28], human reliability analysis[29], cloud computing[30],
energy planning[31, 32, 33], site selection[34, 35, 36], software
assessment[37, 38] and performance evaluation[39, 40]. AHP can be extended
further by other theories and methods such as fuzzy number[41, 42], choquet
integral[43], k-means algorithm[44], VIKOR[45] and TOPSIS[46, 47].

Generally, the process of applying AHP is divided into three steps. Firstly,
establish a hierarchical structure by recursively decomposing the decision
problem. Secondly, construct the pairwise comparison matrix to indicate the
relative importance of alternatives. A numerical rating including nine rank
scales is suggested, as shown in Tab.1. Thirdly, verify the consistency of
pairwise comparison matrix and calculate the priority weights of alternatives.
For the completeness of explanation, several basic concepts and formulas are
introduced as follows.

Table 1: Numerical rating in AHP Scale | Meaning
---|---
1 | Equal importance
3 | Moderate importance
5 | Strong importance
7 | Demonstrated importance
9 | Extreme importance
2, 4, 6, 8 | Intermediate values

###### Definition 1

(Pairwise comparison matrix) Assume
$\begin{Bmatrix}E_{1},E_{2},\cdots,E_{n}\end{Bmatrix}$ are $n$ alternatives
available for decision-making, the pairwise comparison matrix is indicated as
$M=\left(m_{ij}\right)_{n\times n}\left(i,j=1,2,\cdots,n\right)$, and
satisfies:

$m_{ij}=\left\\{\begin{matrix}\frac{1}{m{ji}}&i\neq j\\\
1&i=j\end{matrix}\right.$ (1)

where $m_{ij}$ represents the relative importance of $E_{i}$ over $E_{j}$.

Consistency checking is introduced in AHP to verify the usability of
constructed pairwise comparison matrix.

###### Definition 2

(Consistency ratio) For the pairwise comparison matrix $M_{n\times n}$, let
$\lambda_{\max}$ denote the largest eigenvalue of $M$, consistency index
($CI$) is defined as

$CI=\frac{{\lambda_{\max}-n}}{{n-1}}$ (2)

Based on $CI$, consistency ratio (CR) is defined as

$CR=\frac{{CI}}{{RI}}$ (3)

where $RI$ is the random consistency index related to the dimension of matrix,
listed in Tab.2.

Table 2: Random consistency index $RI$ n | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10
---|---|---|---|---|---|---|---|---|---|---
$RI$ | 0 | 0 | 0.52 | 0.89 | 1.12 | 1.26 | 1.36 | 1.41 | 1.46 | 1.49

If $CR<0.1$, the constructed pairwise comparison matrix is considered
acceptable and the priority weights of alternatives can be obtained through
the eigenvector as Def.3. Otherwise, the comparison matrix needs to be
reconstructed.

###### Definition 3

(Eigenvector) For the pairwise comparison matrix $M_{n\times n}$ with
acceptable consistency, suppose $\vec{w}=({w_{1},w_{2},\cdots,w_{n}})^{T}$ is
the eigenvector of $M$, whose $w_{i}(i=1,2,\cdots,n)$ is indicated as the
priority weight of the $i$-th alternative and calculated by

$M\vec{w}=\lambda_{\max}\vec{w}$ (4)

###  2.2 DEMATEL method

Decision-making trial and evaluation laboratory (DEMATEL) method was
originally developed by Battelle Memorial Institute of Geneva Research Center,
aimed at the fragmented and antagonistic phenomena of world societies and
searched for integrated solution[48, 49]. DEMATEL is built on the basis of
graph theory, effectively assessing the total relations and constructing a map
between system components in respect to their type and severity[50]. Reviewing
the former studies, DEMATEL has been successfully applied in many diverse
areas such as service quality analysis[51, 52], stock selection[53] and supply
chain management[54, 55, 56]. DEMATEL can be further extended by other
theories and methods such as fuzzy numbers[57, 58] and ANP[59, 60, 61, 62].

The procedure of DEMATEL is divided into five steps.

_Step 1: Define quality feature and establish measurement scale_

Quality feature is a set of influential characteristics that impact the
sophisticated system, which can be determined by literature review, brain-
storming, or expert evaluation. After defining the influential characteristics
in researching system, establish the measurement scale for the causal
relationship and pairwise comparison among influential characteristics. Four
levels 0,1,2,3 are suggested, respectively meaning “no impact”, “low impact”,
“high impact” and “extreme high impact”. In this step, factors and their
direct relations are displayed by a weighted and directed graph.

_Step 2: Extract the direct relation matrix of influential factors_

In this step, transformation from the weighted and directed graph into direct
relation matrix has been done. For $n$ influential characteristics
$F_{1},F_{2},\cdots,F_{n}$, direct relation matrix is denoted as
$D=\left(d_{ij}\right)_{n\times n}\left(i,j=1,2,\cdots,n\right)$, where
$d_{ij}$ is the direct relation of $F_{i}$ over $F_{j}$ based on the
measurement scale and satisfies $d_{ij}=0$ if $i=j$.

_Step 3: Normalize the direct relation matrix_

The normalized direct relations of factors are a mapping from $d_{ij}$ to
$\left[0,1\right]$, defined as below.

###### Definition 4

(Normalized direct relation matrix) For the framework of $n$ influential
characteristics $\begin{Bmatrix}F_{1},F_{2},\cdots,F_{n}\end{Bmatrix}$,
normalized matrix $N$ of direct relation matrix
$D=\left(d_{ij}\right)_{n\times n}\left(i,j=1,2,\cdots,n\right)$ is obtained
by

$N=\frac{D}{\max(\mathop{}\sum\limits_{j=1}^{n}{{d_{ij}},}\mathop{}\sum\limits_{i=1}^{n}{{d_{ij}}})}$
(5)

Note that there must be at least one row of matrix $D$ satisfying
$\sum_{j=1}^{n}d_{ij}<{\max(\mathop{}\sum\limits_{j=1}^{n}{{d_{ij}},}\mathop{}\sum\limits_{i=1}^{n}{{d_{ij}}})}$[63].
Hence, based on the studies of [64] and [65], sub-stochastic matrix can be
obtained through utilizing the normalized direct relation matrix $N$ and
absorbing state of Markov chain matrices, shown in (6).

$\mathop{\lim}\limits_{k\to\infty}N^{k}=O\quad
and\quad\mathop{\lim}\limits_{k\to\infty}(I+N+N^{2}+\cdots+N^{k})=N(I-N)^{-1}$
(6)

where $O$ is the null matrix and $I$ is the identity matrix.

_Step 4: Calculate the total relation matrix_

Due to the characteristic of normalized direct relation matrix $N$, total
relation matrix which contains direct and indirect relations among factors can
be derived from (7).

###### Definition 5

(Total relation matrix) For the framework of $n$ influential characteristics
$\begin{Bmatrix}F_{1},F_{2},\cdots,F_{n}\end{Bmatrix}$, assume
$N=\left(n_{ij}\right)_{n\times n}\left(i,j=1,2,\cdots,n\right)$ is the
normalized direct relation matrix, total relation matrix $T$ is defined as

$T=\mathop{\lim}\limits_{k\to\infty}(N+N^{2}+\cdots+N^{k})=N(I-N)^{-1}$ (7)

where $I$ is the identity matrix.

_Step 5: Construct the DEMATEL map_

Based on the sum of each row $R_{i}(i=1,2,\cdots,n)$ and column $C_{i}$
$(i=1,2,\cdots,n)$ of the total relation matrix $T_{n\times n}$,
$(R_{i}+C_{i})$ and $(R_{i}-C_{i})$ can be obtained. $(R_{i}+C_{i})$ is
defined as the prominence, indicating the importance of the i-th influential
factor. $(R_{i}-C_{i})$ classifies the i-th influential factor into the cause
or effect category in researching system. If $(R_{i}-C_{i})>0$, the factor is
regarded as a cause factor. If $(R_{i}-C_{i})<0$, the factor is an effect one.
The DEMATEL map is shown as Fig.1.

Figure 1: DEMATEL map Figure 2: Procedure of DEMATEL-based completion method

##  3 DEMATEL-based completion method

In this section, a DEMATEL-based completion method for incomplete pairwise
comparison matrix in AHP is proposed. Moreover, an application for ranking top
tennis players is presented.

###  3.1 Procedure of DEMATEL-based completion method

In this sub-section, the procedure of DEMATEL-based completion method is
presented and illustrated by an example. Fig.2 shows the process of the
proposed method.

Assume an $4\times 4$ incomplete pairwise comparison matrix $M$ and the
unavailable/missing values are displayed by ‘*’.

$M=\left[{\begin{array}[]{*{20}c}1&*&4&8\\\ &1&2&4\\\ 1/4&1/2&1&2\\\
1/8&1/4&1/2&1\end{array}}\right]$

_Step 1: From incomplete pairwise comparison matrix to direct relation
matrix._

Pairwise comparison matrix reflects the preference relations between each pair
of factors, where $m_{ij}(i,j=1,2,3,4)$ in $M$ indicates the relative
importance of factor $i$ to factor $j$. Hence, the known values in $M$ can be
put into the direct relation matrix $D$ without any transformation and the
unavailable ones(‘*’) are substituted by 0. The direct relation matrix $D$ is
as:

$D=\left[{\begin{array}[]{*{30}c}1&0&4&8\\\ 0&1&2&4\\\ 1/4&1/2&1&2\\\
1/8&1/4&1/2&1\end{array}}\right]$

_Step 2: From direct relation matrix to total relation matrix. _

First of all, normalize the direct relation matrix $D$ according to (5). Then,
based on the normalized matrix $N$ and (7), total relation matrix $T$ can be
obtained. The normalized direct relation matrix $N$ and total relation matrix
$T$ are calculated as:

$N=\left[{\begin{array}[]{*{30}c}0.0667&0.0000&0.2667&0.5333\\\
0.0000&0.0667&0.1333&0.2667\\\ 0.0167&0.0333&0.0667&0.1333\\\
0.0083&0.0167&0.0333&0.0667\end{array}}\right],T=\left[{\begin{array}[]{*{30}c}0.0835&0.0241&0.3371&0.6742\\\
0.0060&0.0835&0.1685&0.3371\\\ 0.0211&0.0421&0.0899&0.1798\\\
0.0105&0.0211&0.0449&0.0899\end{array}}\right]$

_Step 3: From total relation matrix to pairwise comparison matrix._

In fact, the total relation matrix has completed the unavailable values in
initial incomplete matrix. Nevertheless, in pairwise comparison matrix, the
multiplication of each pair of symmetric values along the diagonal line is
required to equal 1, yet the values in total relation matrix are between 0 and
1. Based on (8), it is feasible to accomplish this transformation from the
total relation matrix $T=(t_{ij})_{n\times n}(i,j=1,2,\cdots,n)$ to pairwise
comparison matrix $M_{c}=(c_{ij})_{n\times n}(i,j=1,2,\cdots,n)$.

$\left\\{{\begin{array}[]{*{20}l}\frac{t_{ij}}{c_{ij}}=\frac{t_{ji}}{c_{ji}}\\\
c_{ij}=\frac{1}{c_{ji}}\end{array}}\right.$ (8)

Hence, the pairwise comparison matrix transformed from $T_{4\times 4}$ is:

$M_{c}=\left[{\begin{array}[]{*{30}c}1.0000&2.0000&4.0000&8.0000\\\
0.5000&1.0000&2.0000&4.0000\\\ 0.2500&0.5000&1.0000&2.0000\\\
0.1250&0.2500&0.5000&1.0000\end{array}}\right]$

Finally, based on the proposed method, the completion for incomplete pairwise
comparison matrix $M$ has been accomplished as:

$M=\quad\left[{\begin{array}[]{*{20}c}1&*&4&8\\\ &1&2&4\\\ 1/4&1/2&1&2\\\
1/8&1/4&1/2&1\end{array}}\right]\rightarrow\left[{\begin{array}[]{*{20}c}1&\textbf{2}&4&8\\\
\textbf{1/2}&1&2&4\\\ 1/4&1/2&1&2\\\ 1/8&1/4&1/2&1\end{array}}\right]$

###  3.2 An application for ranking tennis players

Professional tennis as a global popular sport owns large numbers of viewers
around the world. The professional tennis associations such as ATP have been
collecting data about the tournaments and the players. There is a free
available database that collects the results of the top tennis players since
1973 (see at http://www.atpw-
orldtour.com/Players/Head-To-Head.aspx).

In original statistics, 25 tennis players who have been the champion on ATP
since 1973 are selected. The statistical competition results matrix is shown
in Tab.3. However, the statistical matrix is not complete since some of
players did not have the chance to compete with each other or the partial data
has been lost. Thus, DEMATEL-based completion method is adopted to rank the
tennis players based on these statistics. The process is divided into three
steps, which is shown in Fig.3.

Figure 3: Flowchart for ranking tennis players Figure 4: Rule of
transformation

Table 3: Statistical competition results matrix: the wins/total number of matches | Aga | Bec | Bor | Con | Cou | Djo | Edb | Fed | Fer | Hew | Kaf | Kue | Len | McE | Moy | Mus | Nad | Nas | New | Raf | Rio | Rod | Saf | Sam | Wil
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---
Agassi |  | 10/14 |  | 2/2 | 5/12 |  | 6/9 | 3/11 | 2/5 | 4/8 | 8/12 | 7/11 | 2/8 | 2/4 | 3/4 | 5/9 | 0/2 |  |  | 10/15 | 1/3 | 5/6 | 3/6 | 14/34 | 5/7
Becker | 4/14 |  |  | 6/6 | 6/7 |  | 25/35 |  |  | 1/1 | 4/6 |  | 10/21 | 8/10 | 2/4 | 2/3 |  | 1/1 |  | 2/3 | 3/5 |  | 0/1 | 7/19 | 7/10
Borg |  |  |  | 15/23 |  |  |  |  |  |  |  |  | 6/8 | 7/14 |  |  |  | 10/15 | 1/4 |  |  |  |  |  | 1/1
Connors | 0/2 | 0/6 | 8/23 |  | 0/3 |  | 6/12 |  |  |  |  |  | 13/34 | 14/34 |  |  |  | 12/27 | 2/4 |  |  |  |  | 0/2 | 0/5
Courier | 7/12 | 1/7 |  | 3/3 |  |  | 6/10 |  |  |  | 1/6 | 1/1 | 0/4 | 2/3 | 2/3 | 7/12 |  | 0/1 |  | 0/3 | 0/3 |  | 1/2 | 4/20 |
Djokovic |  |  |  |  |  |  |  | 15/31 | 2/3 | 6/7 |  |  |  |  | 2/4 |  | 17/39 |  |  |  |  | 4/9 | 0/2 |  |
Edberg | 3/9 | 10/35 |  | 6/12 | 4/10 |  |  |  |  |  | 1/3 |  | 14/27 | 6/13 | 1/1 | 10/10 |  |  |  | 3/3 | 1/1 |  |  | 6/14 | 9/20
Federer | 8/11 |  |  |  |  | 16/31 |  |  | 10/13 | 18/26 | 2/6 | 1/3 |  |  | 7/7 |  | 10/32 |  |  | 0/3 | 2/2 | 21/24 | 10/12 | 1/1 |
Ferrero | 3/5 |  |  |  |  | 1/3 |  | 3/13 |  | 4/10 | 1/3 | 3/5 |  |  | 8/14 |  | 2/9 |  |  | 2/3 | 3/4 | 0/5 | 6/12 |  |
Hewitt | 4/8 | 0/1 |  |  |  | 1/7 |  | 8/26 | 6/10 |  | 7/8 | 3/4 |  |  | 7/12 |  | 4/10 |  |  | 3/4 | 3/5 | 7/14 | 7/14 | 5/9 |
Kafelnikov | 4/12 | 2/6 |  |  | 5/6 |  | 2/3 | 4/6 | 2/3 | 1/8 |  | 5/12 |  |  | 3/6 | 1/5 |  |  |  | 3/5 | 6/8 |  | 2/4 | 2/13 | 1/2
Kuerten | 4/11 |  |  |  | 0/1 |  |  | 2/3 | 2/5 | 1/4 | 7/12 |  |  |  | 4/7 | 3/3 |  |  |  | 4/8 | 2/4 | 1/2 | 4/7 | 1/3 |
Lendl | 6/8 | 11/21 | 2/8 | 21/34 | 4/4 |  | 13/27 |  |  |  |  |  |  | 21/36 |  | 4/5 |  | 1/1 |  | 0/1 |  |  |  | 3/8 | 15/22
McEnroe | 2/4 | 2/10 | 7/14 | 20/34 | 1/3 |  | 7/13 |  |  |  |  |  | 15/36 |  |  |  |  | 6/9 | 1/2 |  |  |  |  | 0/3 | 7/13
Moya | 1/4 | 2/4 |  |  | 1/3 | 2/4 | 0/1 | 0/7 | 6/14 | 5/12 | 3/6 | 3/7 |  |  |  | 4/8 | 2/8 |  |  | 3/4 | 2/7 | 1/5 | 4/7 | 1/4 |
Muster | 4/9 | 1/3 |  |  | 5/12 |  | 0/10 |  |  |  | 4/5 | 0/3 | 1/5 |  | 4/8 |  |  |  |  | 0/3 | 3/4 |  | 0/1 | 2/11 | 0/2
Nad | 2/2 |  |  |  |  | 22/39 |  | 22/32 | 7/9 | 6/10 |  |  |  |  | 6/8 |  |  |  |  |  |  | 7/10 | 2/2 |  |
Nastase |  | 0/1 | 5/15 | 15/27 | 1/1 |  |  |  |  |  |  |  | 0/1 | 3/9 |  |  |  |  | 4/5 |  |  |  |  |  | 0/1
Newcombe |  |  | 3/4 | 2/4 |  |  |  |  |  |  |  |  |  | 1/2 |  |  |  | 1/5 |  |  |  |  |  |  |
Rafter | 5/15 | 1/3 |  |  | 3/3 |  | 0/3 | 3/3 | 1/3 | 1/4 | 2/5 | 4/8 | 1/1 |  | 1/4 | 3/3 |  |  |  |  | 2/3 |  | 1/1 | 4/16 | 1/3
Rios | 2/3 | 2/5 |  |  | 3/3 |  | 0/1 | 0/2 | 1/4 | 2/5 | 2/8 | 2/4 |  |  | 5/7 | 1/4 |  |  |  | 1/3 |  | 0/2 | 1/4 | 0/2 |
Roddick | 1/6 |  |  |  |  | 5/9 |  | 3/24 | 5/5 | 7/14 |  | 1/2 |  |  | 4/5 |  | 3/10 |  |  |  | 2/2 |  | 4/7 | 2/3 |
Safin | 3/6 | 1/1 |  |  | 1/2 | 2/2 |  | 2/12 | 6/12 | 7/14 | 2/4 | 3/7 |  |  | 3/7 | 1/1 | 0/2 |  |  | 0/1 | 3/4 | 3/7 |  | 4/7 |
Sampras | 20/34 | 12/19 |  | 2/2 | 16/20 |  | 8/14 | 0/1 |  | 4/9 | 11/13 | 2/3 | 5/8 | 3/3 | 3/4 | 9/11 |  |  |  | 12/16 | 2/2 | 1/3 | 3/7 |  | 2/3
Wilander | 2/7 | 3/10 | 0/1 | 5/5 |  |  | 11/20 |  |  |  | 1/2 |  | 7/22 | 6/13 |  | 2/2 |  | 1/1 |  | 2/3 |  |  |  | 1/3 |
_Note:_ Due to the space limitation, the names of tennis players have been
abbreviated to the first three letters in the horizontal row and their full
names can be seen in the vertical column.

Table 4: Incompete pairwise comparison matrix | Aga | Bec | Bor | Con | Cou | Djo | Edb | Fed | Fer | Hew | Kaf | Kue | Len | McE | Moy | Mus | Nad | Nas | New | Raf | Rio | Rod | Saf | Sam | Wil
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---
Agassi |  | 1.39 |  | 1.07 | 0.90 |  | 1.17 | 0.76 | 0.95 | 1.00 | 1.24 | 1.17 | 0.80 | 1.00 | 1.12 | 1.05 | 0.93 |  |  | 1.31 | 0.95 | 1.28 | 1.00 | 0.73 | 1.18
Becker | 0.72 |  |  | 1.38 | 1.38 |  | 2.28 |  |  | 1.04 | 1.18 |  | 0.95 | 1.43 | 1.00 | 1.05 |  | 1.03 |  | 1.05 | 1.05 |  | 0.97 | 0.77 | 1.24
Borg |  |  |  | 1.45 |  |  |  |  |  |  |  |  | 1.25 | 1.00 |  |  |  | 1.31 | 0.89 |  |  |  |  |  | 1.03
Connors | 0.93 | 0.73 | 0.69 |  | 0.88 |  | 1.00 |  |  |  |  |  | 0.66 | 0.73 |  |  |  | 0.86 | 1.00 |  |  |  |  | 0.93 | 0.78
Courier | 1.11 | 0.72 |  | 1.13 |  |  | 1.11 |  |  |  | 0.78 | 1.03 | 0.83 | 1.05 | 1.05 | 1.11 |  | 0.97 |  | 0.88 | 0.88 |  | 1.00 | 0.49 |
Djokovic |  |  |  |  |  |  |  | 0.95 | 1.06 | 1.38 |  |  |  |  | 1.00 |  | 0.77 |  |  |  |  | 0.95 | 0.93 |  |
Edberg | 0.85 | 0.44 |  | 1.00 | 0.90 |  |  |  |  |  | 0.95 |  | 1.05 | 0.95 | 1.03 | 1.89 |  |  |  | 1.13 | 1.03 |  |  | 0.90 | 0.90
Federer | 1.32 |  |  |  |  | 1.05 |  |  | 1.49 | 1.72 | 0.90 | 0.95 |  |  | 1.48 |  | 0.52 |  |  | 0.88 | 1.07 | 3.31 | 1.64 | 1.03 |
Ferrero | 1.05 |  |  |  |  | 0.95 |  | 0.67 |  | 0.90 | 0.95 | 1.05 |  |  | 1.11 |  | 0.75 |  |  | 1.05 | 1.12 | 0.78 | 1.00 | 0 |
Hewitt | 1.00 | 0.97 |  |  |  | 0.72 |  | 0.58 | 1.11 |  | 1.49 | 1.12 |  |  | 1.11 |  | 0.90 |  |  | 1.12 | 1.05 | 1.00 | 1.00 | 1.05 |
Kafelnikov | 0.81 | 0.90 |  |  | 1.28 |  | 1.05 | 1.11 | 1.05 | 0.67 |  | 0.90 |  |  | 1.00 | 0.84 |  |  |  | 1.05 | 1.25 |  | 1.00 | 0.57 | 1.00
Kuerten | 0.85 |  |  |  | 0.97 |  |  | 1.05 | 0.95 | 0.89 | 1.11 |  |  |  | 1.05 | 1.13 |  |  |  | 1.00 | 1.00 | 1.00 | 1.05 | 0.95 |
Lendl | 1.25 | 1.05 | 0.80 | 1.52 | 1.20 |  | 0.95 |  |  |  |  |  |  | 1.36 |  | 1.19 |  | 1.03 |  | 0.97 |  |  |  | 0.90 | 1.54
McEnroe | 1.00 | 0.70 | 1.00 | 1.36 | 0.95 |  | 1.05 |  |  |  |  |  | 0.73 |  |  |  |  | 1.17 | 1.00 |  |  |  |  | 0.88 | 1.05
Moya | 0.89 | 1.00 |  |  | 0.95 | 1.00 | 0.97 | 0.67 | 0.90 | 0.90 | 1.00 | 0.95 |  |  |  | 1.00 | 0.80 |  |  | 1.12 | 0.85 | 0.84 | 1.05 | 0.89 |
Muster | 0.95 | 0.95 |  |  | 0.90 |  | 0.53 |  |  |  | 1.19 | 0.88 | 0.84 |  | 1.00 |  |  |  |  | 0.88 | 1.12 |  | 0.97 | 0.65 | 0.93
Nadal | 1.07 |  |  |  |  | 1.29 |  | 1.91 | 1.34 | 1.11 |  |  |  |  | 1.25 |  |  |  |  |  |  | 1.24 | 1.07 |  |
Nastase |  | 0.97 | 0.77 | 1.17 | 1.03 |  |  |  |  |  |  |  | 0.96 | 0.85 |  |  |  |  | 1.19 |  |  |  |  |  | 0.97
Newcombe |  |  | 1.12 | 1.00 |  |  |  |  |  |  |  |  |  | 1.00 |  |  |  | 0.84 |  |  |  |  |  |  |
Rafter | 0.77 | 0.95 |  |  | 1.13 |  | 0.88 | 1.13 | 0.95 | 0.89 | 0.95 | 1.00 | 1.04 |  | 0.89 | 1.13 |  |  |  |  | 1.05 |  | 1.03 | 0.64 | 0.95
Rios | 1.05 | 0.95 |  |  | 1.13 |  | 0.97 | 0.93 | 0.89 | 0.95 | 0.80 | 1.00 |  |  | 1.18 | 0.89 |  |  |  | 0.95 |  | 0.93 | 0.89 | 0.93 |
Roddick | 0.78 |  |  |  |  | 1.05 |  | 0.30 | 1.28 | 1.00 |  | 1.00 |  |  | 1.19 |  | 0.80 |  |  |  | 1.07 |  | 1.05 | 1.05 |
Safin | 1.00 | 1.03 |  |  | 1.00 | 1.07 |  | 0.61 | 1.00 | 1.00 | 1.00 | 0.95 |  |  | 0.95 | 1.03 | 0.93 |  |  | 0.97 | 1.12 | 0.95 |  | 1.05 |
Sampras | 1.36 | 1.30 |  | 1.07 | 2.04 |  | 1.11 | 0.97 |  | 0.95 | 1.77 | 1.05 | 1.11 | 1.13 | 1.12 | 1.53 |  |  |  | 1.57 | 1.07 | 0.95 | 0.95 |  | 1.05
Wilander | 0.85 | 0.80 | 0.97 | 1.28 |  |  | 1.11 |  |  |  | 1.00 |  | 0.65 | 0.95 |  | 1.07 |  | 1.03 |  | 1.05 |  |  |  | 0.95 |
_Note:_ Due to the space limitation, the names of tennis players have been
abbreviated to the first three letters in the horizontal row and their full
names can be seen in the vertical column.

Table 5: Complete pairwise comparison matrix | Aga | Bec | Bor | Con | Cou | Djo | Edb | Fed | Fer | Hew | Kaf | Kue | Len | McE | Moy | Mus | Nad | Nas | New | Raf | Rio | Rod | Saf | Sam | Wil
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---
Agassi | 1.00 | 1.39 | 0.89 | 1.07 | 0.90 | 1.00 | 1.17 | 0.76 | 0.95 | 1.00 | 1.24 | 1.17 | 0.80 | 1.00 | 1.12 | 1.05 | 0.93 | 1.06 | 1.08 | 1.31 | 0.95 | 1.28 | 1.00 | 0.73 | 1.18
Becker | 0.72 | 1.00 | 0.99 | 1.38 | 1.38 | 0.96 | 2.28 | 0.83 | 1.04 | 1.04 | 1.18 | 1.04 | 0.95 | 1.43 | 1.00 | 1.05 | 0.81 | 1.03 | 1.26 | 1.05 | 1.05 | 1.04 | 0.97 | 0.77 | 1.24
Borg | 1.13 | 1.01 | 1.00 | 1.45 | 1.25 | 1.08 | 1.22 | 0.94 | 1.15 | 1.06 | 1.19 | 1.16 | 1.25 | 1.00 | 1.20 | 1.30 | 0.92 | 1.31 | 0.89 | 1.19 | 1.17 | 1.19 | 1.13 | 0.99 | 1.03
Connors | 0.93 | 0.73 | 0.69 | 1.00 | 0.88 | 0.84 | 1.00 | 0.72 | 0.89 | 0.82 | 0.94 | 0.90 | 0.66 | 0.73 | 0.92 | 1.00 | 0.72 | 0.86 | 1.00 | 0.92 | 0.89 | 0.93 | 0.86 | 0.93 | 0.78
Courier | 1.11 | 0.72 | 0.80 | 1.13 | 1.00 | 0.89 | 1.11 | 0.74 | 0.92 | 0.83 | 0.78 | 1.03 | 0.83 | 1.05 | 1.05 | 1.11 | 0.75 | 0.97 | 1.04 | 0.88 | 0.88 | 0.92 | 1.00 | 0.49 | 0.94
Djokovic | 1.00 | 1.04 | 0.92 | 1.20 | 1.12 | 1.00 | 1.10 | 0.95 | 1.05 | 1.38 | 1.11 | 1.04 | 0.95 | 1.11 | 1.00 | 1.12 | 0.77 | 1.09 | 1.14 | 1.09 | 1.08 | 0.95 | 0.93 | 0.93 | 1.10
Edberg | 0.85 | 0.44 | 0.82 | 1.00 | 0.90 | 0.91 | 1.00 | 0.79 | 0.96 | 0.86 | 0.95 | 0.97 | 1.05 | 0.95 | 1.03 | 1.89 | 0.75 | 0.91 | 0.97 | 1.13 | 1.03 | 0.97 | 0.94 | 0.90 | 0.90
Federer | 1.32 | 1.21 | 1.06 | 1.39 | 1.36 | 1.05 | 1.26 | 1.00 | 1.49 | 1.72 | 0.90 | 0.95 | 1.10 | 1.29 | 1.48 | 1.32 | 0.52 | 1.27 | 1.32 | 0.88 | 1.07 | 3.31 | 1.64 | 1.03 | 1.23
Ferrero | 1.05 | 0.96 | 0.87 | 1.12 | 1.09 | 0.95 | 1.04 | 0.67 | 1.00 | 0.90 | 0.95 | 1.05 | 0.90 | 1.03 | 1.11 | 1.07 | 0.75 | 1.03 | 1.07 | 1.05 | 1.12 | 0.78 | 1.00 | 0.82 | 1.03
Hewitt | 1.00 | 0.97 | 0.95 | 1.22 | 1.20 | 0.72 | 1.17 | 0.58 | 1.11 | 1.00 | 1.49 | 1.12 | 0.97 | 1.14 | 1.11 | 1.17 | 0.90 | 1.10 | 1.17 | 1.12 | 1.05 | 1.00 | 1.00 | 1.05 | 1.13
Kafelnikov | 0.81 | 0.90 | 0.84 | 1.07 | 1.28 | 0.90 | 1.05 | 1.11 | 1.05 | 0.67 | 1.00 | 0.90 | 0.83 | 0.97 | 1.00 | 0.84 | 0.74 | 0.98 | 1.02 | 1.05 | 1.25 | 1.00 | 1.00 | 0.57 | 1.00
Kuerten | 0.85 | 0.96 | 0.86 | 1.11 | 0.97 | 0.96 | 1.03 | 1.05 | 0.95 | 0.89 | 1.11 | 1.00 | 0.89 | 1.03 | 1.05 | 1.13 | 0.78 | 1.02 | 1.06 | 1.00 | 1.00 | 1.00 | 1.05 | 0.95 | 1.03
Lendl | 1.25 | 1.05 | 0.80 | 1.52 | 1.20 | 1.05 | 0.95 | 0.91 | 1.12 | 1.03 | 1.20 | 1.12 | 1.00 | 1.36 | 1.15 | 1.19 | 0.91 | 1.03 | 1.18 | 0.97 | 1.13 | 1.15 | 1.10 | 0.90 | 1.54
McEnroe | 1.00 | 0.70 | 1.00 | 1.36 | 0.95 | 0.90 | 1.05 | 0.78 | 0.97 | 0.88 | 1.03 | 0.97 | 0.73 | 1.00 | 1.00 | 1.10 | 0.78 | 1.17 | 1.00 | 1.01 | 0.96 | 1.00 | 0.93 | 0.88 | 1.05
Moya | 0.89 | 1.00 | 0.84 | 1.08 | 0.95 | 1.00 | 0.97 | 0.67 | 0.90 | 0.90 | 1.00 | 0.95 | 0.87 | 1.00 | 1.00 | 1.00 | 0.80 | 0.98 | 1.03 | 1.12 | 0.85 | 0.84 | 1.05 | 0.89 | 0.99
Muster | 0.95 | 0.95 | 0.77 | 1.00 | 0.90 | 0.89 | 0.53 | 0.76 | 0.93 | 0.86 | 1.19 | 0.88 | 0.84 | 0.91 | 1.00 | 1.00 | 0.74 | 0.91 | 0.94 | 0.88 | 1.12 | 0.93 | 0.97 | 0.65 | 0.93
Nadal | 1.07 | 1.23 | 1.08 | 1.39 | 1.34 | 1.29 | 1.33 | 1.91 | 1.34 | 1.11 | 1.36 | 1.29 | 1.10 | 1.28 | 1.25 | 1.35 | 1.00 | 1.28 | 1.33 | 1.34 | 1.32 | 1.24 | 1.07 | 1.12 | 1.32
Nastase | 0.94 | 0.97 | 0.77 | 1.17 | 1.03 | 0.92 | 1.10 | 0.79 | 0.97 | 0.91 | 1.02 | 0.98 | 0.96 | 0.85 | 1.02 | 1.10 | 0.78 | 1.00 | 1.19 | 1.01 | 0.99 | 0.99 | 0.96 | 0.80 | 0.97
Newcombe | 0.93 | 0.79 | 1.12 | 1.00 | 0.96 | 0.88 | 1.03 | 0.76 | 0.94 | 0.86 | 0.98 | 0.94 | 0.85 | 1.00 | 0.97 | 1.06 | 0.75 | 0.84 | 1.00 | 0.97 | 0.94 | 0.96 | 0.91 | 0.82 | 0.96
Rafter | 0.77 | 0.95 | 0.84 | 1.09 | 1.13 | 0.92 | 0.88 | 1.13 | 0.95 | 0.89 | 0.95 | 1.00 | 1.04 | 0.99 | 0.89 | 1.13 | 0.74 | 0.99 | 1.03 | 1.00 | 1.05 | 1.02 | 1.03 | 0.64 | 0.95
Rios | 1.05 | 0.95 | 0.85 | 1.12 | 1.13 | 0.92 | 0.97 | 0.93 | 0.89 | 0.95 | 0.80 | 1.00 | 0.89 | 1.04 | 1.18 | 0.89 | 0.76 | 1.01 | 1.07 | 0.95 | 1.00 | 0.93 | 0.89 | 0.93 | 1.00
Roddick | 0.78 | 0.96 | 0.84 | 1.08 | 1.08 | 1.05 | 1.03 | 0.30 | 1.28 | 1.00 | 1.00 | 1.00 | 0.87 | 1.00 | 1.19 | 1.07 | 0.80 | 1.01 | 1.04 | 0.98 | 1.07 | 1.00 | 1.05 | 1.05 | 1.00
Safin | 1.00 | 1.03 | 0.88 | 1.16 | 1.00 | 1.07 | 1.06 | 0.61 | 1.00 | 1.00 | 1.00 | 0.95 | 0.91 | 1.07 | 0.95 | 1.03 | 0.93 | 1.04 | 1.10 | 0.97 | 1.12 | 0.95 | 1.00 | 1.05 | 1.05
Sampras | 1.36 | 1.30 | 1.01 | 1.07 | 2.04 | 1.07 | 1.11 | 0.97 | 1.22 | 0.95 | 1.77 | 1.05 | 1.11 | 1.13 | 1.12 | 1.53 | 0.89 | 1.25 | 1.22 | 1.57 | 1.07 | 0.95 | 0.95 | 1.00 | 1.05
Wilander | 0.85 | 0.80 | 0.97 | 1.28 | 1.06 | 0.91 | 1.11 | 0.81 | 0.97 | 0.88 | 1.00 | 0.97 | 0.65 | 0.95 | 1.01 | 1.07 | 0.76 | 1.03 | 1.05 | 1.05 | 1.00 | 1.00 | 0.96 | 0.95 | 1.00
_Note:_ Due to the space limitation, the names of tennis players have been
abbreviated to the first three letters in the horizontal row and their full
names can be seen in the vertical column.

Firstly, the statistical competition results matrix (i.e. with additive
preference relations) is converted into the pairwise comparison matrix (i.e.
with multiplicative preference relations) following the method described in
[66]. Fig.4 illustrates this transformation process. The converted pairwise
comparison matrix is shown in Tab.4.

Secondly, DEMATEL-based completion method is adopted to complete the pairwise
comparison matrix. Following the procedure of the proposed method, the
complete pairwise comparison matrix is obtained, which is shown in Tab.5.

Thirdly, based on the complete pairwise comparison matrix, the priorities of
tennis players can be calculated via (4). Tab.6 shows the final ranking of the
25 tennis players.

Table 6: Ranking of 25 tennis players Rank | Player | Priority | Rank | Player | Priority
---|---|---|---|---|---
1 | Nadal | 0.0503 | 14 | Rios | 0.0380
2 | Federer | 0.0503 | 15 | McEnroe | 0.0380
3 | Sampras | 0.0467 | 16 | Nastase | 0.0380
4 | Borg | 0.0444 | 17 | Rafter | 0.0378
5 | Lendl | 0.0436 | 18 | Wilander | 0.0378
6 | Becker | 0.0430 | 19 | Kafelnikov | 0.0375
7 | Hewitt | 0.0414 | 20 | Edberg | 0.0374
8 | Djokovic | 0.0412 | 21 | Moya | 0.0371
9 | Agassi | 0.0409 | 22 | Newcombe | 0.0365
10 | Safin | 0.0392 | 23 | Courier | 0.0360
11 | Kuerten | 0.0390 | 24 | Muster | 0.0353
12 | Roddick | 0.0384 | 25 | Connors | 0.0339
13 | Ferrero | 0.0383 |  |  |

##  4 Discussion

The result in last section indicates that the proposed method has ability to
complete the pairwise comparison matrices in AHP. In this section, we would
further analyze the superiority and efficiency of the proposed method.

The proposed method has the following advantages:

  * 1.

The proposed method estimates the missing values in pairwise comparison matrix
from a new perspective (i.e. by calculating the total relations among
alternatives).

  * 2.

The proposed method has low computational cost.

After reviewing the prior studies, an assessment for the computational cost of
the existing completion methods and the proposed method is done. The methods
for comparison satisfy:

  * 1.

The methods belong to completion methods, which have ability to estimate the
missing values in pairwise comparison matrix.

  * 2.

The pairwise comparison matrices, which can be completed by these completion
methods, present reciprocal preference relations among alternatives. In other
words, the matrices request that the multiplication of each pair of symmetric
values along the diagonal line is equal to 1.

As for these completion methods, some of them estimate the missing values
through iteration, or based on the complex theories/methods, which causes the
methods expensive computational cost. For example, Bozoíki _et al._[14]
developed a completion method which estimates the missing values in pairwise
comparison matrix through iteration and optimizes the matrix consistency based
on graph theory. Gomez _et al._[20] proposed a completion method based on the
neural network. Neural network as one of the branches of artificial
intelligence leads to the complexity of the completion method. As for the
completion methods without iteration or complex theories/algorithms such as
the methods listed in Tab.7, time complexity is used to accurately reflect
their computational cost. Assume there is an $m\times m$ incomplete pairwise
comparison matrix, which has n missing values, the time complexity of these
methods is shown in Tab.8.

Table 7: Methods for comparison |  Description |  Source
---|---|---
Method 1 |  | A completion method extends the geometric mean induced bias matrix to est-
---
imate the missing judgments and improves the consistency ratios by the least
absolute error method and the least square method.
Ergu _et al._[13]
Method 2 |  | A consistency-based completion procedure, which takes a measure to maintain
---
the consistency level of matrices before and after completion.
Alonso _et al._[18]
Method 3 |  | A completion method to obtain the most consistent/similar matrix with the in-
---
itial incomplete matrix based on Frobenius norm.
Benítez _et al._[19]
Table 8: Time complexity of methods | Method 1 | Method 2 | Method 3 | The proposed method
---|---|---|---|---
Time complexity | $m^{3}$ | $m^{3}$ | $m^{2}+mn$ | $m^{2}$

It can be seen from the Tab.8 that the time complexity of the proposed method
is lower than the others. In addition, the most existing completion methods
have the following two operations:

  * 1.

Divide the matrix into upper and lower triangular matrix. The pairwise
comparison matrix in AHP has a property that the multiplication of each pair
of symmetric values along the diagonal line is equal to 1. These methods
firstly focus on the completion of upper/lower triangle of incomplete matrix
and use the property to accomplish the other triangular matrix.

  * 2.

Traverse the pairwise comparison matrix to obtain the positions of unavailable
values, then estimate these missing values in sequence.

However, the proposed method completes the matrix from a holistic viewpoint
without dividing the matrix into the upper and lower triangular matrix.
Furthermore, based on DEMATEL, all the missing values in matrix can be
simultaneously estimated. Hence, following the procedure of the proposed
method, the completion of incomplete pairwise comparison matrix becomes
straightforward and efficient.

##  5 Conclusion

Pairwise comparison matrix plays a pivotal role in AHP. However, in many
cases, only partial information in pairwise comparison matrix is available,
which obstructs the subsequent operations of the classical AHP. In this paper,
we propose a new completion method for incomplete pairwise comparison matrix
in AHP. Based on DEMATEL, the proposed method has ability to complete the
matrix by obtaining the total (direct and indirect) relations from the direct
preference relations between each pair of alternatives. Furthermore, an
application of the proposed method for ranking tennis players is presented.
From the ranking result and the comparison with several representative
completion methods, the proposed method is proved to be effective and
efficient. The proposed method provides a new perspective to complete the
matrix with explicit physical meaning. Besides, the proposed method has low
computation cost. This promising method has a wide application in multi-
criteria decision-making. In our further study, we would extend the proposed
method for the incomplete fuzzy pairwise comparison matrices.

## Acknowledgment

The work is partially supported by National Natural Science Foundation of
China (Grant Nos. 61573290,61503237), China State Key Laboratory of Virtual
Reality Technology and Systems, Beihang University (Grant No.BUAA-VR-14KF-02).

## References

## References

  * [1] A. Ishizaka, A. Labib, Analytic hierarchy process and expert choice: Benefits and limitations, OR Insight 22 (4) (2009) 201–220.
  * [2] B. Zhu, Z. Xu, R. Zhang, M. Hong, Hesitant analytic hierarchy process, European Journal of Operational Research 250 (2) (2016) 602–614.
  * [3] C. T. Lin, Y. H. Goay, Measuring National Tourism Organization Abroad Office Competitiveness, International Journal of Tourism Research 17 (2) (2015) 118–129.
  * [4] V. Yepes, T. García-Segura, J. Moreno-Jiménez, A cognitive approach for the multi-objective optimization of RC structural problems, Archives of Civil and Mechanical Engineering 15 (4) (2015) 1024–1036.
  * [5] Y.-M. Wang, K.-S. Chin, A linear goal programming approach to determining the relative importance weights of customer requirements in quality function deployment, Information Sciences 181 (24) (2011) 5523–5533.
  * [6] A. Ishizaka, A. Labib, Review of the main developments in the analytic hierarchy process, Expert systems with applications 38 (11) (2011) 14336–14345.
  * [7] P. T. Harker, Incomplete pairwise comparisons in the analytic hierarchy process, Mathematical Modelling 9 (11) (1987) 837–848.
  * [8] F. J. Carmone, A. Kara, S. H. Zanakis, A Monte Carlo investigation of incomplete pairwise comparison matrices in AHP, European Journal of Operational Research 102 (3) (1997) 538–553.
  * [9] Z. Hua, B. Gong, X. Xu, A DS–AHP approach for multi-attribute decision making problem with incomplete information, Expert systems with applications 34 (3) (2008) 2221–2227.
  * [10] Y. Xu, R. Patnayakuni, H. Wang, Logarithmic least squares method to priority for group decision making with incomplete fuzzy preference relations, Applied Mathematical Modelling 37 (4) (2013) 2139–2152.
  * [11] E. Dopazo, M. Ruiz-Tagle, A parametric GP model dealing with incomplete information for group decision-making, Applied Mathematics and Computation 218 (2) (2011) 514–519.
  * [12] M. Fedrizzi, S. Giove, Incomplete pairwise comparison and consistency optimization, European Journal of Operational Research 183 (1) (2007) 303–313.
  * [13] D. Ergu, G. Kou, Y. Peng, M. Zhang, Estimating the missing values for the incomplete decision matrix and consistency optimization in emergency management, Applied Mathematical Modelling 40 (1) (2016) 254–267.
  * [14] S. Bozóki, J. Fülöp, L. Rónyai, On optimal completion of incomplete pairwise comparison matrices, Mathematical and Computer Modelling 52 (1) (2010) 318–333.
  * [15] K. Chen, G. Kou, J. M. Tarn, Y. Song, Bridging the gap between missing and inconsistent values in eliciting preference from pairwise comparison matrices, Annals of Operations Research 235 (1) (2015) 155–175.
  * [16] G. Zhang, Y. Dong, Y. Xu, Linear optimization modeling of consistency issues in group decision making based on fuzzy preference relations, Expert Systems with Applications 39 (3) (2012) 2415–2420.
  * [17] Y.-C. Hu, J.-F. Tsai, Backpropagation multi-layer perceptron for incomplete pairwise comparison matrices in analytic hierarchy process, Applied Mathematics and Computation 180 (1) (2006) 53–62.
  * [18] S. Alonso, F. Chiclana, F. Herrera, E. Herrera-Viedma, J. Alcalá-Fdez, C. Porcel, A consistency-based procedure to estimate missing pairwise preference values, International Journal of Intelligent Systems 23 (2) (2008) 155–175.
  * [19] J. Benítez, X. Delgado-Galván, J. Izquierdo, R. Pérez-García, Consistent completion of incomplete judgments in decision making using AHP, Journal of Computational and Applied Mathematics 290 (2015) 412–422.
  * [20] J. A. Gomez-Ruiz, M. Karanik, J. I. Peláez, Estimation of missing judgments in AHP pairwise matrices using a neural network-based model, Applied Mathematics and Computation 216 (10) (2010) 2959–2975.
  * [21] G. Büyüközkan, G. Çifçi, A new incomplete preference relations based approach to quality function deployment, Information Sciences 206 (2012) 30–41.
  * [22] F. J. Cabrerizo, I. J. Pérez, E. Herrera-Viedma, Managing the consensus in group decision making in an unbalanced fuzzy linguistic context with incomplete information, Knowledge-Based Systems 23 (2) (2010) 169–181.
  * [23] T. L. Saaty, The analytic hierarchy process: planning, priority setting, resources allocation, New York: McGraw, 1980.
  * [24] Z. Xu, Deviation square priority method for distinct preference structures based on generalized multiplicative consistency, Fuzzy Systems, IEEE Transactions on 23 (4) (2015) 1164–1180.
  * [25] A. Ishizaka, P. Nemery, Multi-criteria decision analysis: methods and software, John Wiley & Sons, 2013.
  * [26] F. T. Chan, N. Kumar, Global supplier development considering risk factors using fuzzy extended AHP-based approach, Omega 35 (4) (2007) 417–431.
  * [27] X. Deng, Y. Hu, Y. Deng, S. Mahadevan, Supplier selection using AHP methodology extended by D numbers, Expert Systems with Applications 41 (1) (2014) 156–167.
  * [28] F. T. Chan, H. Chan, H. C. Lau, R. W. Ip, An AHP approach in benchmarking logistics performance of the postal industry, Benchmarking: An International Journal 13 (6) (2006) 636–661.
  * [29] X. Su, S. Mahadevan, P. Xu, Y. Deng, Dependence assessment in human reliability analysis using evidence theory and AHP, Risk Analysis 35 (7) (2015) 1296–1316.
  * [30] S. Liu, F. T. Chan, W. Ran, Decision making for the selection of cloud vendor: An improved approach under group decision-making with integrated weights and objective/subjective attributes, Expert Systems with Applications 55 (2016) 37–47.
  * [31] A. Ishizaka, S. Siraj, P. Nemery, Which energy mix for the UK (United Kingdom)? An evolutive descriptive mapping with the integrated GAIA (graphical analysis for interactive aid)–AHP (analytic hierarchy process) visualization tool, Energy 95 (2016) 602–611.
  * [32] S. C. Onar, B. Oztaysi, İ. Otay, C. Kahraman, Multi-expert wind energy technology selection using interval-valued intuitionistic fuzzy sets, Energy 90 (1) (2015) 274–285.
  * [33] C. Kahraman, İ. Kaya, S. Cebi, A comparative analysis for multiattribute selection among renewable energy alternatives using fuzzy axiomatic design and fuzzy analytic hierarchy process, Energy 34 (10) (2009) 1603–1616.
  * [34] M. Erdoğan, İ. Kaya, A combined fuzzy approach to determine the best region for a nuclear power plant in Turkey, Applied Soft Computing 39 (2016) 84–93.
  * [35] A. H. Lee, H.-Y. Kang, C.-Y. Lin, K.-C. Shen, An Integrated Decision-Making Model for the Location of a PV Solar Plant, Sustainability 7 (10) (2015) 13522–13541.
  * [36] A. Pourahmad, A. Hosseini, A. Banaitis, H. Nasiri, N. Banaitienė, G.-H. Tzeng, Combination of fuzzy-AHP and DEMATEL-ANP with gis in a new hybrid mcdm model used for the selection of the best space for leisure in a blighted urban site, Technological and Economic Development of Economy 21 (5) (2015) 773–796.
  * [37] Y. Hu, X. Zhang, E. Ngai, R. Cai, M. Liu, Software project risk analysis using bayesian networks with causality constraints, Decision Support Systems 56 (2013) 439–449.
  * [38] Y. Hu, B. Feng, X. Mo, X. Zhang, E. Ngai, M. Fan, M. Liu, Cost-sensitive and ensemble-based prediction model for outsourced software project risk prediction, Decision Support Systems 72 (2015) 11–23.
  * [39] E. K. Zavadskas, Z. Turskis, J. Antucheviciene, Selecting a Contractor by Using a Novel Method for Multiple Attribute Analysis: Weighted Aggregated Sum Product Assessment with Grey Values (WASPAS-G), Studies in Informatics and Control 24 (2) (2015) 141–150.
  * [40] C.-R. Wu, C.-T. Lin, P.-H. Tsai, Evaluating business performance of wealth management banks, European Journal of Operational Research 207 (2) (2010) 971–979.
  * [41] F. T. Chan, N. Kumar, M. Tiwari, H. Lau, K. Choy, Global supplier selection: a fuzzy-AHP approach, International Journal of Production Research 46 (14) (2008) 3825–3857.
  * [42] Y.-M. Wang, Y. Luo, Z. Hua, On the extent analysis method for fuzzy AHP and its applications, European Journal of Operational Research 186 (2) (2008) 735–747.
  * [43] S. Corrente, S. Greco, A. Ishizaka, Combining analytical hierarchy process and choquet integral within non-additive robust ordinal regression, Omega 61 (2016) 2–18.
  * [44] F. Lolli, A. Ishizaka, R. Gamberini, New AHP-based approaches for multi-criteria inventory classification, International Journal of Production Economics 156 (2014) 62–74.
  * [45] G. Büyüközkan, A. Görener, Evaluation of product development partners using an integrated AHP-VIKOR model, Kybernetes 44 (2) (2015) 220–237.
  * [46] M. Erdogan, I. Kaya, Evaluating Alternative-Fuel Busses for Public Transportation in Istanbul Using Interval Type-2 Fuzzy AHP and TOPSIS, Journal of multiple-valued logic and soft computing 26 (6) (2016) 625–642.
  * [47] C. Kahraman, A. Suder, E. T. Bekar, Fuzzy multiattribute consumer choice among health insurance options, Technological and Economic Development of Economy 22 (1) (2016) 1–20.
  * [48] A. Gabus, E. Fontela, World problems, an invitation to further thought within the framework of DEMATEL, Battelle Geneva Research Center, Geneva, Switzerland, 1972.
  * [49] E. Fontela, A. Gabus, The DEMATEL observer, Battelle Geneva Research Centre, Geneva, 1976.
  * [50] G.-H. Tzeng, C.-H. Chiang, C.-W. Li, Evaluating intertwined effects in e-learning programs: A novel hybrid MCDM model based on factor analysis and DEMATEL, Expert systems with Applications 32 (4) (2007) 1028–1044.
  * [51] J.-I. Shieh, H.-H. Wu, K.-K. Huang, A DEMATEL method in identifying key success factors of hospital service quality, Knowledge-Based Systems 23 (3) (2010) 277–282.
  * [52] M.-L. Tseng, A causal and effect decision making model of service quality expectation using grey-fuzzy DEMATEL approach, Expert Systems with Applications 36 (4) (2009) 7738–7748.
  * [53] K.-Y. Shen, G.-H. Tzeng, Combined soft computing model for value stock selection based on fundamental analysis, Applied Soft Computing 37 (2015) 142–155.
  * [54] J. J. H. Liou, J. Tamosaitiene, E. K. Zavadskas, G.-H. Tzeng, New hybrid COPRAS-G MADM Model for improving and selecting suppliers in green supply chain management, International Journal of Production Research 54 (1, SI) (2016) 114–134.
  * [55] H.-H. Wu, S.-Y. Chang, A case study of using DEMATEL method to identify critical factors in green supply chain management, Applied Mathematics and Computation 256 (2015) 394–403.
  * [56] K.-J. Wu, C.-J. Liao, M.-L. Tseng, A. S. F. Chiu, Exploring decisive factors in green supply chain practices under uncertainty, International Journal of Production Economics 159 (SI) (2015) 147–157.
  * [57] A. Mentes, H. Akyildiz, M. Yetkin, N. Turkoglu, A FSA based fuzzy DEMATEL approach for risk assessment of cargo ships at coasts and open seas of Turkey, Safety Science 79 (2015) 1–10.
  * [58] S.-B. Tsai, M.-F. Chien, Y. Xue, L. Li, X. Jiang, Q. Chen, J. Zhou, L. Wang, Using the Fuzzy DEMATEL to Determine Environmental Performance: A Case of Printed Circuit Board Industry in Taiwan, PloS one 10 (6) (2015) e0129153.
  * [59] J. J. H. Liou, Building an effective system for carbon reduction management, Journal of Cleaner Production 103 (2015) 353–361.
  * [60] G.-H. Tzeng, C.-Y. Huang, Combined DEMATEL technique with hybrid MCDM methods for creating the aspired intelligent global manufacturing & logistics systems, Annals of Operations Research 197 (1) (2012) 159–190.
  * [61] W.-H. Tsai, W.-C. Chou, Selecting management systems for sustainable development in SMEs: A novel hybrid model based on DEMATEL, ANP, and ZOGP, Expert Systems with Applications 36 (2) (2009) 1444–1458.
  * [62] W.-W. Wu, Choosing knowledge management strategies by using a combined ANP and DEMATEL approach, Expert Systems with Applications 35 (3) (2008) 828–835.
  * [63] C.-J. Lin, W.-W. Wu, A causal analytical method for group decision-making under fuzzy environment, Expert Systems with Applications 34 (1) (2008) 205–213.
  * [64] A. Papoulis, S. U. Pillai, Probability, random variables, and stochastic processes, McGraw-Hill, 1985.
  * [65] R. Goodman, Introduction to stochastic models, Courier Corporation, 1988.
  * [66] S. Bozóki, L. Csató, J. Temesi, An application of incomplete pairwise comparison matrices for ranking top tennis players, European Journal of Operational Research 248 (1) (2016) 211–218.


# Berry Curvature Spectroscopy from Bloch Oscillations

Christophe De Beule  Department of Physics and Astronomy, University of
Pennsylvania, Philadelphia PA 19104  Department of Physics and Materials
Science, University of Luxembourg, L-1511 Luxembourg, Luxembourg  E. J. Mele
Department of Physics and Astronomy, University of Pennsylvania, Philadelphia
PA 19104

###### Abstract

We demonstrate that the Berry curvature of an isolated Bloch miniband in two-
dimensional superlattices can be probed by the _dressed_ linear optical
response when a uniform static field is applied to the system. In particular,
when the static field is sufficiently strong such that full Bloch oscillations
occur before the crystal momentum relaxes to equilibrium, the optical response
of the dressed system becomes resonant at the Bloch frequencies. The latter
are in the THz regime when the superlattice periodicity is of the order of
$10\;\text{nm}$. Using a band-projected semiclassical theory, we define a
dressed optical conductivity and find that the height of the resonances in the
dressed _Hall_ conductivity are proportional to the Fourier components of the
Berry curvature. We illustrate our results with a low-energy model on an
effective honeycomb lattice.

Nonlinear optical responses are becoming an increasingly important tool to
investigate the spectral and geometric properties of electron Bloch bands in
low-dimensional materials [1, 2, 3]. In particular, the nonlinear Hall effect
[4] which probes the multipoles of the Berry curvature of the band at
successive orders in the driving field [5]. The distribution of the Berry
curvature can be studied by nonlinear response even in crystals with time-
reversal symmetry, which show no linear Hall effect, since time-reversal
symmetry only precludes odd powers of the field in the Hall response.
Recently, the advent of moiré [6, 7, 8] and other two-dimensional (2D)
superlattice materials [9, 10, 11] has opened up the prospect of studying
responses at nonperturbative order in the driving field [12, 13, 14]. These
systems can host spectrally isolated and flattened minibands, and nonlinear
responses have already been used to study their properties [15, 16, 17, 18,
19, 20, 21, 22]. Moreover, because the real-space periodicity of such
superlattices can be much larger than the underlying atomic scale, with
lattice constants ranging between $1$–$100$ nm, the momentum space
superlattice Brillouin zone (BZ) is relatively small. Under an applied
electric field, it therefore becomes possible that an electron traverses the
_entire_ zone, i.e., performs a full Bloch oscillation [23, 24], before
relaxing back to equilibrium by scattering processes. To quantify this,
consider an applied uniform electric field of the form

$\bm{E}(t)=\bm{E}_{0}+\bm{E}_{1}(t),$ (1)

which has a static component
$\bm{E}_{0}=E_{0}\left(\cos\theta_{0},\,\sin\theta_{0}\right)$ and an
oscillating component $\bm{E}_{1}(t)$. The latter acts as a weak probe for the
system that is _dressed_ by the static field. Here, the nonperturbative regime
is defined by the condition $\omega_{B}\tau\gg 1$ [12, 13, 14] where
$\omega_{B}=eE_{0}L/\hbar$ is the characteristic frequency of Bloch
oscillations, which we call the Bloch frequency, and $\tau$ is the momentum-
relaxation time with $L$ the lattice constant. If we estimate
$\tau=1\;\text{ps}$ we find that
$\omega_{B}\tau\approx\tfrac{1.5E_{0}}{\text{kV/cm}}\tfrac{L}{10\;\text{nm}}$
such that $\omega_{B}\tau$ can become large in superlattice materials for
reasonable field strengths [12, 13].

Figure 1: (a) 2D superlattice (e.g., a moiré) subjected to a static uniform
in-plane electric field $\bm{E}_{0}$ and probed by monochromatic light of
frequency $\omega$. (b) Berry curvature $\Omega_{\bm{k}}$ in the first-shell
approximation for a system with $D_{3}$ or $C_{3v}$ symmetry. (c) Imaginary
part of the dressed optical Hall conductivity $\sigma_{H}(\omega,\bm{E}_{0})$
for the Berry curvature shown in (b) as a function of $\omega/\omega_{B}$ for
$\theta_{0}=15^{\circ}$ and different values of $\omega_{B}\tau$. (d)
$\text{Im}\,\sigma_{H}$ in units $e^{2}\Omega_{1}f_{1}^{0}/2\hbar V_{c}$ for
$\omega_{B}\tau=15$ as a function of the frequency and the field direction
$\theta_{0}$. The resonant frequencies for the first shell
$\omega_{n}^{(1)}=|e\bm{E}_{0}\cdot\bm{L}_{n}/\hbar|$ are indicated.

In this work, we study the dressed time-dependent response of time-reversal-
invariant 2D superlattices with lattice constants $L\sim 10\;\text{nm}$, that
are subjected to a uniform electric field of the form given in Eq. (1). This
setup is illustrated in Fig. 1(a). When the static field is in the
nonperturbative regime, we find an optical response, linear in the
_oscillating_ component, that is resonant at the Bloch frequencies. For the
studied systems, the latter are on the order of $10\;\text{THz}$. Moreover, we
show that the peak heights of the resonances in the dressed optical _Hall_
conductivity are proportional to the Fourier components of the Berry
curvature. Hence our approach is in some sense dual to probing the momentum-
space distribution of the Berry curvature via its multipoles at successive
harmonics [25] and complementary to other methods that study orbital moments
with circular dichroism [26]. In contrast, in our proposal, all information on
the Berry curvature is contained in the dressed _linear_ optical response.

_Semiclassical theory._ Our starting point is the band-projected semiclassical
theory of electron dynamics for a 2D crystal in a uniform electric field
$\bm{E}(t)$. The equations of motion for the central position and crystal
momentum of a wave packet constructed from the Bloch states of an energy band
$\varepsilon_{n\bm{k}}$ are given by [27, 28]

$\displaystyle\hbar\dot{\bm{r}}_{n\bm{k}}$
$\displaystyle=\nabla_{\bm{k}}\varepsilon_{n\bm{k}}-\hbar\dot{\bm{k}}\times\Omega_{n\bm{k}}\hat{z},$
(2) $\displaystyle\hbar\dot{\bm{k}}$ $\displaystyle=-e\bm{E}(t),$ (3)

where $-e$ is the electron charge and
$\Omega_{n\bm{k}}=-2\,\text{Im}\,\langle\partial_{k_{x}}u_{n\bm{k}}|\partial_{k_{y}}u_{n\bm{k}}\rangle_{\text{cell}}$
is the Berry curvature 111To construct a wave packet
$\left|W(t)\right>=\int_{\bm{k}}c_{\bm{k}}(t)\left|\Psi_{\bm{k}}\right>$, the
Bloch states
$\left|\Psi_{\bm{k}}\right>=e^{i\bm{k}\cdot\hat{\bm{r}}}\left|u_{\bm{k}}\right>$
should be smooth on the BZ torus, i.e.,
$\left|\Psi_{\bm{k}+\bm{G}}\right>=\left|\Psi_{\bm{k}}\right>$ with $\bm{G}$ a
reciprocal lattice vector. This is periodic gauge [33] and yields
$\left|u_{\bm{k}+\bm{G}}\right>=e^{-i\bm{G}\cdot\hat{\bm{r}}}\left|u_{\bm{k}}\right>$,
in contrast to
$\left|\tilde{u}_{\bm{k}+\bm{G}}\right>=\left|\tilde{u}_{\bm{k}}\right>$ for
which the Bloch Hamiltonian is periodic (Bloch form). The Berry curvature is
generally different in both gauges.. The band-projected theory holds as long
as interband transitions can be neglected. These can arise both from optical
transitions and electric breakdown (Zener tunneling) [30]. The former are
absent for frequencies well below the energy gap to the other energy bands
$\varepsilon_{\text{gap}}$, while the absence of the latter can be estimated
by the condition that
$\varepsilon_{\text{gap}}^{2}/\varepsilon_{\text{width}}\gg eE_{0}L$ where
$\varepsilon_{\text{width}}$ is the bandwidth [13, 14].

In the following, we drop the band index $n$ since we consider a single band.
The current is then given by

$\bm{j}(t)=-e\int_{\bm{k}}\,\dot{\bm{r}}_{\bm{k}}(t)f_{\bm{k}}(t),$ (4)

with $\int_{\bm{k}}=\int_{\text{BZ}}d^{2}\bm{k}/(2\pi)^{2}$ and where
$f_{\bm{k}}(t)$ is the nonequilibrium occupation of the electrons in the band.
The latter is obtained from the Boltzmann transport equation in the
relaxation-time approximation:

$\tau\partial_{t}f_{\bm{k}}-\frac{e\tau}{\hbar}\bm{E}(t)\cdot\nabla_{\bm{k}}f_{\bm{k}}=f_{\bm{k}}^{0}-f_{\bm{k}},$
(5)

where $\tau$ is the momentum-relaxation time and
$f_{\bm{k}}^{0}=n_{F}(\varepsilon_{\bm{k}}-\mu)$ with $n_{F}$ the Fermi
function and $\mu$ the chemical potential. Because the system has
translational symmetry, the occupation function is periodic in momentum space:
$f_{\bm{k}}=\sum_{\bm{R}}f_{\bm{R}}e^{i\bm{k}\cdot\bm{R}}$ where the sum runs
over lattice vectors $\bm{R}$ with
$f_{\bm{R}}=V_{c}\int_{\bm{k}}f_{\bm{k}}e^{-i\bm{k}\cdot\bm{R}}$. Plugging
this expansion in Eq. (S5) we obtain an ordinary differential equation with
the steady-state solution [31]

$f_{\bm{R}}(t)=f^{0}_{\bm{R}}\int_{0}^{\infty}ds\,e^{-s}\exp\left[\frac{ie}{\hbar}\int_{t-s\tau}^{t}dt^{\prime}\,\bm{E}(t^{\prime})\cdot\bm{R}\right],$
(6)

as shown in the Supplemental Material (SM) [32]. The occupation $f_{\bm{k}}$
is thus given by a weighted sum of displaced Fermi functions where the drift
due to the electric field is determined by the accumulated momentum between
collisions at time $t-s\tau$ and time $t$. The exponential weight $e^{-s}$
reflects the fact that scattering is modeled as a Poisson process in the
relaxation-time approximation.

The current in Eq. (4) can be decomposed into two terms as
$\bm{j}(t)=\bm{j}_{\text{Bloch}}(t)+\bm{j}_{\text{geom}}(t)$ where

$\displaystyle\bm{j}_{\text{Bloch}}(t)$ $\displaystyle=\frac{ie}{\hbar
V_{c}}\sum_{\bm{R}}\bm{R}\varepsilon_{-\bm{R}}f_{\bm{R}}(t),$ (7)
$\displaystyle\bm{j}_{\text{geom}}(t)$
$\displaystyle=\hat{z}\times\frac{e^{2}}{\hbar
V_{c}}\sum_{\bm{R}}\Omega_{-\bm{R}}\bm{E}(t)f_{\bm{R}}(t),$ (8)

where $V_{c}$ is the unit cell area and we made use of the expansions of the
band dispersion and the Berry curvature, as well as
$V_{c}\int_{\bm{k}}e^{i\bm{k}\cdot\bm{R}}=\delta_{\bm{R},\bm{0}}$. The Bloch
current $\bm{j}_{\text{Bloch}}$ originates from the band dispersion while the
geometric current $\bm{j}_{\text{geom}}$ originates from the anomalous
velocity due to the Berry curvature in Eq. (2).

_Dressed optical conductivity._ We now consider probing the system by
monochromatic light of frequency $\omega$ at normal incidence. In the
electric-dipole approximation, the electric field of the light can be written
as

$\bm{E}_{1}(t)=\bm{\mathcal{E}}_{1}e^{i\omega
t}+\bm{\mathcal{E}}_{1}^{*}e^{-i\omega t},$ (9)

where $\bm{\mathcal{E}}_{1}\in\mathds{C}^{2}$ gives the polarization. To
investigate the response of the system at frequency $\omega$, we expand the
lattice Fourier components of the distribution function in terms of their
frequency components. We have
$f_{\bm{R}}(t)=\sum_{m=-\infty}^{\infty}f_{\bm{R},m}e^{im\omega t}$ where
$f_{\bm{R},m}=(\omega/2\pi)\int_{0}^{2\pi/\omega}dt\,f_{\bm{R}}(t)e^{-im\omega
t}$ with $f_{\bm{R},-m}=f_{-\bm{R},m}^{*}$. The frequency components of the
currents become

$\displaystyle\bm{j}_{\text{Bloch}}^{(m)}$ $\displaystyle=\frac{ie}{\hbar
V_{c}}\sum_{\bm{R}}\bm{R}\,\varepsilon_{-\bm{R}}f_{\bm{R},m},$ (10)
$\displaystyle\begin{split}\bm{j}_{\text{geom}}^{(m)}&=\hat{z}\times\frac{e^{2}}{\hbar
V_{c}}\sum_{\bm{R}}\Omega_{-\bm{R}}\\\
&\cdot\left(\bm{E}_{0}f_{\bm{R},m}+\bm{\mathcal{E}}_{1}f_{\bm{R},m-1}+\bm{\mathcal{E}}_{1}^{*}f_{\bm{R},m+1}\right).\end{split}$
(11)

Since we are interested in the linear response _dressed_ by the static part of
the field, we expand Eq. (6) in orders of
$|e\bm{\mathcal{E}}_{1}\cdot\bm{R}/\hbar\omega|$ while retaining all orders in
$\bm{E}_{0}$. Up to first order, the only nonzero terms are given by

$\displaystyle f_{\bm{R},0}$
$\displaystyle=\frac{f^{0}_{\bm{R}}}{1-i\omega_{\bm{R}}\tau},$ (12)
$\displaystyle f_{\bm{R},1}$
$\displaystyle=\frac{f^{0}_{\bm{R}}}{1-i\omega_{\bm{R}}\tau}\frac{e\bm{\mathcal{E}}_{1}\cdot\bm{R}/\hbar}{\omega-\omega_{\bm{R}}-\frac{i}{\tau}}=f_{-\bm{R},-1}^{*},$
(13)

with $\omega_{\bm{R}}=e\bm{E}_{0}\cdot\bm{R}/\hbar$.

Figure 2: Roses for the real (a) and imaginary (b) part of the dressed optical
Hall conductivity $\sigma_{H}(\omega,\bm{E}_{0})$ for the Berry curvature
shown in Fig. 1(b) with $\omega_{B}\tau=15$. The angle corresponds to the
direction of the static electric field $\theta_{0}$ and the color scale gives
the frequency $\omega$ of the oscillating field.

The response at frequency $\omega$ can then be written as
$j_{a}^{(1)}=\sigma_{ab}\mathcal{E}_{1b}$ where $a,b=x,y$ and summation over
repeated indices is implied. This leads us to the main result of this work:
the _dressed_ optical conductivity

$\displaystyle\sigma_{ab}(\omega,\bm{E}_{0})=\frac{ie^{2}}{\hbar^{2}V_{c}}\sum_{\bm{R}}\frac{R_{a}R_{b}\varepsilon_{-\bm{R}}f_{\bm{R}}^{0}}{\left(1-i\omega_{\bm{R}}\tau\right)\left(\omega-\omega_{\bm{R}}-\frac{i}{\tau}\right)}$
(14) $\displaystyle-\frac{e^{2}}{\hbar
V_{c}}\sum_{\bm{R}}\frac{\Omega_{-\bm{R}}f_{\bm{R}}^{0}}{1-i\omega_{\bm{R}}\tau}\left[\epsilon_{ab}+\frac{e\epsilon_{ac}E_{0c}R_{b}/\hbar}{\omega-\omega_{\bm{R}}-\frac{i}{\tau}}\right],$

where $\epsilon_{ab}$ is the permutation symbol and
$\sigma_{ab}(\omega,\bm{E}_{0})^{*}=\sigma_{ab}(-\omega,\bm{E}_{0})$ such that
the real (imaginary) part is even (odd) in $\omega$. As a check, we undress
the conductivity by setting $E_{0}=0$. In this case, the two terms in Eq. (14)
reduce to the Drude and anomalous Hall conductivity, respectively.
Importantly, the dressed linear Hall response does _not_ vanish when time-
reversal symmetry is conserved, because it is effectively a compound nonlinear
response in the fields $\bm{E}_{0}$ and $\bm{E}_{1}(t)$.

Let us now focus on the case where $\bm{E}_{0}$ is finite and consider the
dressed longitudinal $\sigma_{L}=\delta_{ab}\sigma_{ab}/2$ and Hall
$\sigma_{H}=\epsilon_{ab}\sigma_{ab}/2$ conductivities, which transform as a
scalar and pseudoscalar under coordinate transformations, respectively. We
obtain

$\displaystyle\sigma_{L}$
$\displaystyle=\frac{ie^{2}}{2\hbar^{2}V_{c}}\sum_{\bm{R}}\frac{R^{2}\varepsilon_{-\bm{R}}f^{0}_{\bm{R}}}{\left(1-i\omega_{\bm{R}}\tau\right)\left(\omega-\omega_{\bm{R}}-\frac{i}{\tau}\right)},$
(15) $\displaystyle\sigma_{H}$ $\displaystyle=-\frac{e^{2}}{\hbar
V_{c}}\sum_{\bm{R}}\frac{\Omega_{-\bm{R}}f^{0}_{\bm{R}}}{1-i\omega_{\bm{R}}\tau}\left(1+\frac{1}{2}\frac{\omega_{\bm{R}}}{\omega-\omega_{\bm{R}}-\frac{i}{\tau}}\right),$
(16)

which for $\omega_{B}\tau\gg 1$ simplify to

$\displaystyle\sigma_{L}(\omega,\bm{E}_{0})$
$\displaystyle=-\frac{e^{2}}{h}\frac{\pi}{\tau
V_{c}}\sum_{\bm{R}}\frac{R^{2}\varepsilon_{-\bm{R}}f^{0}_{\bm{R}}}{\hbar\omega_{\bm{R}}\left(\omega-\omega_{\bm{R}}-\frac{i}{\tau}\right)},$
(17) $\displaystyle\sigma_{H}(\omega,\bm{E}_{0})$
$\displaystyle=-\frac{e^{2}}{h}\frac{\pi}{\tau
V_{c}}\sum_{\bm{R}}\frac{i\Omega_{-\bm{R}}f^{0}_{\bm{R}}}{\omega-\omega_{\bm{R}}-\frac{i}{\tau}}.$
(18)

For crystals with time-reversal symmetry, the band dispersion (Berry
curvature) is an even (odd) function of momentum, such that
$\varepsilon_{\bm{R}}$ and $f_{\bm{R}}^{0}$ are real, while $\Omega_{\bm{R}}$
is imaginary. In this case, $\text{Im}\,\sigma_{L}$ and
$\text{Im}\,\sigma_{H}$ are given by a series of Lorentzians centered at the
Bloch frequencies $\omega_{\bm{R}}$ in the nonperturbative limit. The height
of these resonances is proportional to $\varepsilon_{\bm{R}}$ and
$\Omega_{\bm{R}}$, respectively, and independent of the relaxation time
$\tau$. Conversely, the real part of the dressed conductivity vanishes at
resonance. Physically, we can interpret this as the transfer of energy from
the incident light into coherent Bloch motion of the electrons. And since we
assume here that $\omega_{B}\tau\gg 1$, the system is essentially
collisionless on the time scale set by the Bloch oscillations.

These results can thus potentially be used to map out the distribution of the
Berry curvature in systems with time-reversal symmetry by measuring the
resonances in the dressed optical Hall conductivity in the nonperturbative
regime where $\omega_{B}\tau\gg 1$.

_First-shell approximation._ It is instructive to first evaluate the dressed
optical conductivity by only taking into account the leading-order terms in
the sum over the lattice vectors. For concreteness, we consider a system with
point group $D_{3}$ or $C_{3v}$ which lacks inversion or $\mathcal{C}_{2z}$
rotation symmetry. In this case, the Berry curvature is generally nonzero even
though the Chern number of the band vanishes. In the first-shell
approximation, we only take into account the shortest nonzero lattice vectors
such that
$\varepsilon_{\bm{k}}=\varepsilon_{1}\sum_{n=1}^{3}\cos(\bm{k}\cdot\bm{L}_{n})$
up to an overall constant and
$\Omega_{\bm{k}}=\Omega_{1}\sum_{n=1}^{3}\sin(\bm{k}\cdot\bm{L}_{n})$ where
$\varepsilon_{1}$ and $\Omega_{1}$ are real parameters that depend on the
details of the system, and $\bm{L}_{1}=L(1/2,\sqrt{3}/2)$,
$\bm{L}_{2}=(-L,0)$, and $\bm{L}_{3}=-(\bm{L}_{1}+\bm{L}_{2})$ are related by
$\mathcal{C}_{3z}$ rotation symmetry [13, 14].

The imaginary part of the dressed optical Hall conductivity is shown in Fig.
1(c) as a function of $\omega$ for $\theta_{0}=15^{\circ}$ and different
values of $\omega_{B}\tau$. There are three resonances in this case because
the first coordination shell supports three Bloch frequencies
$\omega_{n}^{(1)}=|e\bm{E}_{0}\cdot\bm{L}_{n}|$ which are nondegenerate for a
general field direction. The height of these resonances is approximately equal
due to $\mathcal{C}_{3z}$ and time-reversal symmetry and saturates to
$e^{2}\Omega_{1}f_{1}^{0}/2\hbar V_{c}$ in the limit $\omega_{B}\tau\gg 1$,
where $f_{1}^{0}=f_{\bm{R}=\bm{L}_{n}}^{0}$. Notice that the resonances are
only well-defined for $\omega_{B}\tau\gtrsim 10$. The dependence on the
direction of the static field is shown in Fig. 1(c). Here we show
$\text{Im}\,\sigma_{H}$ for $\omega_{B}\tau=15$ as a function of $\omega$ and
$\theta_{0}$. As we rotate the static field, the resonances move along the
curves $\omega=\pm\omega_{B}\cos(\theta_{0}-\theta_{n})$ with
$\theta_{n}=\\{\pi/3,\pi,-\pi/3\\}$. For the special case $\theta_{0}=m\pi/3$
($m\in\mathbb{Z}$) two Bloch frequencies coincide and the peaks are doubled.
On the contrary, for $\theta_{0}=(2m+1)\pi/6$ the response vanishes due to
$\mathcal{M}_{x}$ ($x\mapsto-x$) mirror symmetry. These features can also be
seen in the rose plots of Fig. 2. Here we clearly see that the strongest
resonance occurs when two lattice vectors have the same projection along the
static field. Away from these directions, the resonance splits into two peaks
that shift to higher and lower frequencies.

_Low-energy model._ Going beyond the first-shell approximation, we now
consider a low-energy model defined on an effective honeycomb lattice with one
orbital per site, and with nearest-neighbor hopping amplitude $t>0$ and a
sublattice-staggering potential $m$. The Bloch Hamiltonian is given by

$\displaystyle\mathcal{H}(\bm{k})=\bm{d}(\bm{k})\cdot\bm{\sigma},$ (19)
$\displaystyle\bm{d}(\bm{k})=\left(-t\,\text{Re}\,g_{\bm{k}},\,-t\,\text{Im}\,g_{\bm{k}},\,m\right),$
(20)

where $\bm{\sigma}=\left(\sigma_{x},\sigma_{y},\sigma_{z}\right)$ are the
Pauli matrices and
$g_{\bm{k}}=e^{-i\bm{k}\cdot\bm{\tau}}\left[1+e^{i\bm{k}\cdot\bm{L}_{1}}+e^{i\bm{k}\cdot\left(\bm{L}_{1}+\bm{L}_{2}\right)}\right]$
with $\bm{\tau}=L\hat{y}/\sqrt{3}$ the relative separation of the two
sublattices. Note that we work in periodic gauge for which the semiclassical
equation given in Eq. (2) is valid [28, 33]. This model has time-reversal
symmetry with point group $C_{3v}$ generated by $\mathcal{C}_{3z}$ and
$\mathcal{M}_{x}$, and can be seen as a minimal low-energy model for moirés
such as hBN-aligned twisted bilayer graphene [34, 35] or twisted double
bilayer graphene [36, 37], as well as other superlattice systems belonging to
the same symmetry class such as periodically-buckled graphene with a $C_{3v}$
height profile [38, 39, 40, 14].

Figure 3: Dressed optical conductivities $\sigma_{L}(\omega,\bm{E}_{0})$ and
$\sigma_{H}(\omega,\bm{E}_{0})$ for the valence band of the effective two-band
model with $m/t=0.5$ where $\omega_{B}\tau=15$, $\theta_{0}=0^{\circ}$, and
$k_{B}T/t=0.004$. The color scale gives the filling $\nu$ [see inset of (b)].
(a, b) Real and imaginary part of the longitudinal conductivity. (c, d) Real
and imaginary part of the Hall conductivity. Dashed vertical lines give the
position of the resonances $\omega_{\bm{R}}$ and the inset in (a) and (d)
shows the relative magnitude and phase of $\varepsilon_{\bm{R}-}$ and
$\Omega_{\bm{R}-}$, respectively.

The model gives two energy bands $\varepsilon_{\bm{k}\pm}=\pm|\bm{d}(\bm{k})|$
that are separated by a gap $|2m|$ at the zone corners. Because
$\mathcal{C}_{2z}$ symmetry is broken by the sublattice potential, the Berry
curvature is nonzero and given by

$\Omega_{\bm{k}\pm}=\pm\frac{mt^{2}V_{c}}{6|\bm{d}(\bm{k})|^{3}}\sum_{n=1}^{3}\sin\left(\bm{k}\cdot\bm{L}_{n}\right),$
(21)

with $V_{c}=\sqrt{3}L^{2}/2$. In the limit $|m/t|\gg 1$, we have
$|\bm{d}(\bm{k})|\simeq|m|$ and the first shell dominates with
$\Omega_{1}=\pm\operatorname{sgn}(m)V_{c}\,t^{2}/6m^{2}$. However, in general
many shells contribute, as illustrated in Fig. 3 where we show
$\sigma_{L}(\omega,\bm{E}_{0})$ in panels (a) and (b), and
$\sigma_{H}(\omega,\bm{E}_{0})$ in panels (c) and (d) for $m/t=0.5$ and
different fillings $\nu$ of the valence band. Here the static field lies along
the $x$ direction and $k_{B}T/t\ll 1$. Note $\sigma_{L}$ decays faster with
frequency than $\sigma_{H}$ because the first shell of the dispersion is
dominant [see inset of Fig. 3(a)] as well as the additional factor of
$1/\omega$ in Eq. (17). Note that the filling $\nu$ enters only through the
Fourier components of the Fermi function $f_{\bm{R}}^{0}$ which modulate the
height of the peaks in the imaginary part of the conductivities and can change
sign as a function of $\nu$, see Fig. 3(d).

In conclusion, we developed a band-projected semiclassical theory for the
optical response of a superlattice material that is dressed by a uniform
static field. When the static field is sufficiently strong, which is achieved
for field strengths of the order of $10\;\text{kV/cm}$ for a lattice constant
of order $10\;\text{nm}$, the dressed system becomes resonant at the
frequencies of Bloch oscillations which are in the $10\;\text{THz}$ regime. We
quantified this effect by defining a dressed optical conductivity whose
imaginary part displays peaks at the Bloch frequencies, while the real part
vanishes at resonance. In particular, the height of the resonances in the
optical Hall conductivity probes the Berry curvature distribution of the
electronic band. The dressed optical conductivity can for example be obtained
from THz Faraday rotation and ellipticity spectroscopy measurements [41, 42].
Our work thus provides a novel route to probe the Berry curvature in
superlattice materials with time-reversal symmetry.

###### Acknowledgements.

We thank V. T. Phong for discussions. This research was funded in whole, or in
part, by the Luxembourg National Research Fund (FNR) (project No. 16515716).
CDB and EJM are supported by the Department of Energy under grant DE-
FG02-84ER45118.

## References

  * Morimoto and Nagaosa [2016] T. Morimoto and N. Nagaosa, Topological nature of nonlinear optical effects in solids, Science Advances 2, e1501524 (2016).
  * Wu  _et al._ [2017] L. Wu, S. Patankar, T. Morimoto, N. L. Nair, E. Thewalt, A. Little, J. G. Analytis, J. E. Moore, and J. Orenstein, Giant anisotropic nonlinear optical response in transition metal monopnictide weyl semimetals, Nature Physics 13, 350 (2017).
  * Ahn  _et al._ [2022] J. Ahn, G.-Y. Guo, N. Nagaosa, and A. Vishwanath, Riemannian geometry of resonant optical responses, Nature Physics 18, 290 (2022).
  * Sodemann and Fu [2015] I. Sodemann and L. Fu, Quantum Nonlinear Hall Effect Induced by Berry Curvature Dipole in Time-Reversal Invariant Materials, Phys. Rev. Lett. 115, 216806 (2015).
  * Zhang  _et al._ [2023] C.-P. Zhang, X.-J. Gao, Y.-M. Xie, H. C. Po, and K. T. Law, Higher-order nonlinear anomalous Hall effects induced by Berry curvature multipoles, Phys. Rev. B 107, 115142 (2023).
  * Andrei and MacDonald [2020] E. Y. Andrei and A. H. MacDonald, Graphene bilayers with a twist, Nature Materials 19, 1265 (2020).
  * Andrei  _et al._ [2021] E. Y. Andrei, D. K. Efetov, P. Jarillo-Herrero, A. H. MacDonald, K. F. Mak, T. Senthil, E. Tutuc, A. Yazdani, and A. F. Young, The marvels of moiré materials, Nat. Rev. Mater. 6, 201 (2021).
  * Mak and Shan [2022] K. F. Mak and J. Shan, Semiconductor moiré materials, Nature Nanotechnology 17, 686 (2022).
  * R. [2005] R. R.,  _Superlattice to Nanoelectronics_ (Elsevier Science, 2005).
  * Forsythe  _et al._ [2018] C. Forsythe, X. Zhou, K. Watanabe, T. Taniguchi, A. Pasupathy, P. Moon, M. Koshino, P. Kim, and C. R. Dean, Band structure engineering of 2D materials using patterned dielectric superlattices, Nat. Nanotechnol. 13, 566 (2018).
  * Mao  _et al._ [2020] J. Mao, S. P. Milovanović, M. Anđelković, X. Lai, Y. Cao, K. Watanabe, T. Taniguchi, L. Covaci, F. M. Peeters, A. K. Geim, Y. Jiang, and E. Y. Andrei, Evidence of flat bands and correlated states in buckled graphene superlattices, Nature 584, 215 (2020).
  * Fahimniya  _et al._ [2021] A. Fahimniya, Z. Dong, E. I. Kiselev, and L. Levitov, Synchronizing Bloch-Oscillating Free Carriers in Moiré Flat Bands, Phys. Rev. Lett. 126, 256803 (2021).
  * Phong and Mele [2022a] V. T. Phong and E. J. Mele, Quantum Geometric Oscillations in Two-Dimensional Flat-Band Solids (2022a), arXiv:2212.12852 [cond-mat.mes-hall] .
  * Beule  _et al._ [2023] C. D. Beule, V. T. Phong, and E. J. Mele, Rose patterns in the nonperturbative current response of two-dimensional superlattices (2023), arXiv:2305.03013 [cond-mat.mes-hall] .
  * Pantaleón  _et al._ [2021] P. A. Pantaleón, T. Low, and F. Guinea, Tunable large Berry dipole in strained twisted bilayer graphene, Phys. Rev. B 103, 205403 (2021).
  * He and Weng [2021] Z. He and H. Weng, Giant nonlinear Hall effect in twisted bilayer WTe2, npj Quantum Materials 6, 101 (2021).
  * Sinha  _et al._ [2022] S. Sinha, P. C. Adak, A. Chakraborty, K. Das, K. Debnath, L. D. V. Sangani, K. Watanabe, T. Taniguchi, U. V. Waghmare, A. Agarwal, and M. M. Deshmukh, Berry curvature dipole senses topological transition in a moiré superlattice, Nature Physics 18, 765 (2022).
  * Chakraborty  _et al._ [2022] A. Chakraborty, K. Das, S. Sinha, P. C. Adak, M. M. Deshmukh, and A. Agarwal, Nonlinear anomalous Hall effects probe topological phase-transitions in twisted double bilayer graphene, 2D Materials 9, 045020 (2022).
  * Zhang  _et al._ [2022] C.-P. Zhang, J. Xiao, B. T. Zhou, J.-X. Hu, Y.-M. Xie, B. Yan, and K. T. Law, Giant nonlinear Hall effect in strained twisted bilayer graphene, Phys. Rev. B 106, L041111 (2022).
  * Pantaleón  _et al._ [2022] P. A. Pantaleón, V. o. T. Phong, G. G. Naumis, and F. Guinea, Interaction-enhanced topological Hall effects in strained twisted bilayer graphene, Phys. Rev. B 106, L161101 (2022).
  * Duan  _et al._ [2022] J. Duan, Y. Jian, Y. Gao, H. Peng, J. Zhong, Q. Feng, J. Mao, and Y. Yao, Giant Second-Order Nonlinear Hall Effect in Twisted Bilayer Graphene, Phys. Rev. Lett. 129, 186801 (2022).
  * Zhong  _et al._ [2023] J. Zhong, J. Duan, S. Zhang, H. Peng, Q. Feng, Y. Hu, Q. Wang, J. Mao, J. Liu, and Y. Yao, Effective manipulation and realization of a colossal nonlinear Hall effect in an electric-field tunable moiré system (2023), arXiv:2301.12117 [cond-mat.mes-hall] .
  * Bloch [1929] F. Bloch, Über die Quantenmechanik der Elektronen in Kristallgittern, Zeitschrift für Physik 52, 555 (1929).
  * Leo  _et al._ [1992] K. Leo, P. H. Bolivar, F. Brüggemann, R. Schwedler, and K. Köhler, Observation of Bloch oscillations in a semiconductor superlattice, Solid State Communications 84, 943 (1992).
  * Luu and Wörner [2018] T. T. Luu and H. J. Wörner, Measurement of the Berry curvature of solids using high-harmonic spectroscopy, Nature Communications 9, 916 (2018).
  * Schüler  _et al._ [2020] M. Schüler, U. D. Giovannini, H. Hübener, A. Rubio, M. A. Sentef, and P. Werner, Local Berry curvature signatures in dichroic angle-resolved photoelectron spectroscopy from two-dimensional materials, Science Advances 6, eaay2730 (2020).
  * Chang and Niu [1995] M.-C. Chang and Q. Niu, Berry Phase, Hyperorbits, and the Hofstadter Spectrum, Phys. Rev. Lett. 75, 1348 (1995).
  * Sundaram and Niu [1999] G. Sundaram and Q. Niu, Wave-packet dynamics in slowly perturbed crystals: Gradient corrections and Berry-phase effects, Phys. Rev. B 59, 14915 (1999).
  * Note [1] To construct a wave packet $\left|W(t)\right>=\intop\ilimits@_{\bm{k}}c_{\bm{k}}(t)\left|\Psi_{\bm{k}}\right>$, the Bloch states $\left|\Psi_{\bm{k}}\right>=e^{i\bm{k}\cdot\hat{\bm{r}}}\left|u_{\bm{k}}\right>$ should be smooth on the BZ torus, i.e., $\left|\Psi_{\bm{k}+\bm{G}}\right>=\left|\Psi_{\bm{k}}\right>$ with $\bm{G}$ a reciprocal lattice vector. This is periodic gauge [33] and yields $\left|u_{\bm{k}+\bm{G}}\right>=e^{-i\bm{G}\cdot\hat{\bm{r}}}\left|u_{\bm{k}}\right>$, in contrast to $\left|\tilde{u}_{\bm{k}+\bm{G}}\right>=\left|\tilde{u}_{\bm{k}}\right>$ for which the Bloch Hamiltonian is periodic (Bloch form). The Berry curvature is generally different in both gauges.
  * Ashcroft and Mermin [1976] N. W. Ashcroft and N. D. Mermin,  _Solid State Physics_ (Saunders College Publishing, 1976).
  * Mikhailov [2017] S. A. Mikhailov, Nonperturbative quasiclassical theory of the nonlinear electrodynamic response of graphene, Phys. Rev. B 95, 085432 (2017).
  * [32] See Supplemental Material at [link] for a detailed calculation of the occupation function and the dressed optical conductivity.
  * Vanderbilt [2018] D. Vanderbilt,  _Berry Phases in Electronic Structure Theory_ (Cambridge University Press, 2018).
  * Zhang  _et al._ [2019] Y.-H. Zhang, D. Mao, and T. Senthil, Twisted bilayer graphene aligned with hexagonal boron nitride: Anomalous Hall effect and a lattice model, Phys. Rev. Res. 1, 033126 (2019).
  * Lewandowski and Levitov [2019] C. Lewandowski and L. Levitov, Intrinsically undamped plasmon modes in narrow electron bands, Proceedings of the National Academy of Sciences 116, 20869 (2019).
  * Koshino [2019] M. Koshino, Band structure and topological properties of twisted double bilayer graphene, Phys. Rev. B 99, 235406 (2019).
  * Chebrolu  _et al._ [2019] N. R. Chebrolu, B. L. Chittari, and J. Jung, Flat bands in twisted double bilayer graphene, Phys. Rev. B 99, 235417 (2019).
  * Milovanović  _et al._ [2020] S. P. Milovanović, M. Anđelković, L. Covaci, and F. M. Peeters, Band flattening in buckled monolayer graphene, Phys. Rev. B 102, 245427 (2020).
  * Phong and Mele [2022b] V. T. Phong and E. J. Mele, Boundary Modes from Periodic Magnetic and Pseudomagnetic Fields in Graphene, Phys. Rev. Lett. 128, 176406 (2022b).
  * Gao  _et al._ [2022] Q. Gao, J. Dong, P. Ledwith, D. Parker, and E. Khalaf, Untwisting moiré physics: Almost ideal bands and fractional Chern insulators in periodically strained monolayer graphene (2022), arXiv:2211.00658 [cond-mat.mes-hall] .
  * Spielman  _et al._ [1994] S. Spielman, B. Parks, J. Orenstein, D. T. Nemeth, F. Ludwig, J. Clarke, P. Merchant, and D. J. Lew, Observation of the Quasiparticle Hall Effect in Superconducting $\mathrm{Y}{\mathrm{Ba}}_{2}{\mathrm{Cu}}_{3}{\mathrm{O}}_{7-\delta}$, Phys. Rev. Lett. 73, 1537 (1994).
  * Shimano  _et al._ [2011] R. Shimano, Y. Ikebe, K. S. Takahashi, M. Kawasaki, N. Nagaosa, and Y. Tokura, Terahertz Faraday rotation induced by an anomalous Hall effect in the itinerant ferromagnet SrRuO3, Europhysics Letters 95, 17002 (2011).

Supplemental Material for “Berry Curvature Spectroscopy from Bloch
Oscillations"

##  I Semiclassical model of electron dynamics

The semiclassical equations of motion for an electron in a two-dimensional
(2D) crystal, occupying an energy band with dispersion $\varepsilon_{n\bm{k}}$
subjected to a uniform electric field $\bm{E}(t)$ are given by [27, 28]

$\displaystyle\hbar\dot{\bm{r}}_{n\bm{k}}$
$\displaystyle=\nabla_{\bm{k}}\varepsilon_{n\bm{k}}-\hbar\dot{\bm{k}}\times\bm{\Omega}_{n\bm{k}},$
(S1a) $\displaystyle\hbar\dot{\bm{k}}$ $\displaystyle=-e\bm{E}(t),$ (S1b)

where the dot stands for the time derivative $d/dt$. Here $e>0$ is the
elementary charge, $n$ is the band index, and
$\bm{\Omega}_{n\bm{k}}=\Omega_{n\bm{k}}\hat{z}$ is the Berry curvature. The
latter is defined as

$\Omega_{n\bm{k}}=i\left(\left\langle\frac{\partial u_{n\bm{k}}}{\partial
k_{x}}\bigg{\rvert}\frac{\partial u_{n\bm{k}}}{\partial
k_{y}}\right\rangle_{\text{cell}}-\left\langle\frac{\partial
u_{n\bm{k}}}{\partial k_{y}}\bigg{\rvert}\frac{\partial u_{n\bm{k}}}{\partial
k_{x}}\right\rangle_{\text{cell}}\right),$ (S2)

where $u_{n\bm{k}}(\bm{r})$ are cell-periodic Bloch functions in periodic
gauge,
$u_{n,\bm{k}+\bm{G}}(\bm{r})=e^{-i\bm{G}\cdot\bm{r}}u_{n\bm{k}}(\bm{r})$ with
$\bm{G}$ a reciprocal lattice vector, and $\langle
u_{n\bm{k}}|u_{m\bm{k}}\rangle_{\text{cell}}=\delta_{nm}$.

In the following, we consider a single band and omit the band index $n$. The
band dispersion and Berry curvature can be expanded as

$\varepsilon_{\bm{k}}=\sum_{\bm{R}}\varepsilon_{\bm{R}}e^{i\bm{k}\cdot\bm{R}},\qquad\Omega_{\bm{k}}=\sum_{\bm{R}}\Omega_{\bm{R}}e^{i\bm{k}\cdot\bm{R}},$
(S3)

where

$\varepsilon_{\bm{R}}=V_{c}\int_{\bm{k}}\varepsilon_{\bm{k}}e^{-i\bm{k}\cdot\bm{R}},\qquad\Omega_{\bm{R}}=V_{c}\int_{\bm{k}}\Omega_{\bm{k}}e^{-i\bm{k}\cdot\bm{R}},\qquad\int_{\bm{k}}=\int_{\text{BZ}}\frac{d^{2}\bm{k}}{\left(2\pi\right)^{2}}.$
(S4)

##  II Boltzmann transport equation

The Boltzmann equation for the distribution function $f(\bm{k},\bm{r},t)$ in
the relaxation-time approximation, is given by

$\frac{\partial f}{\partial
t}+\frac{d\bm{k}}{dt}\cdot\nabla_{\bm{k}}f+\frac{d\bm{r}}{dt}\cdot\nabla_{\bm{r}}f=\frac{f^{0}-f}{\tau},$
(S5)

where $f^{0}(\bm{k})$ is the equilibrium distribution function, i.e.,
$f^{0}(\bm{k})=n_{F}(\varepsilon_{\bm{k}}-\mu)$ with
$n_{F}(z)=1/(e^{z/k_{B}T}+1)$ the Fermi function, with $\mu$ the chemical
potential and $T$ the temperature.

Let us consider a general uniform time-dependent electric field $\bm{E}(t)$.
We are interested in the steady-state solutions (not necessarily static) of

$\frac{\partial f}{\partial
t}-\frac{e}{\hbar}\bm{E}(t)\cdot\nabla_{\bm{k}}f=\frac{f^{0}-f}{\tau}.$ (S6)

In a translational-invariant system,

$f(\bm{k},t)=\sum_{\bm{R}}f_{\bm{R}}(t)\,e^{i\bm{k}\cdot\bm{R}},$ (S7)

where $\bm{R}$ are lattice vectors, and similarly for $f^{0}(\bm{k})$. We then
obtain an ordinary differential equation for each Fourier component
$f_{\bm{R}}(t)$,

$\frac{df_{\bm{R}}}{dt}+\left(\frac{1}{\tau}-\frac{ie}{\hbar}\bm{E}(t)\cdot\bm{R}\right)f_{\bm{R}}(t)=\frac{f^{0}_{\bm{R}}}{\tau},$
(S8)

whose general solution is given by

$f_{\bm{R}}(t)=\frac{f^{0}_{\bm{R}}}{\tau}\int_{t_{0}}^{t}dt^{\prime}\,e^{-\frac{t-t^{\prime}}{\tau}}\exp\left[\frac{ie}{\hbar}\int_{t^{\prime}}^{t}du\,\bm{E}(u)\cdot\bm{R}\right],$
(S9)

with $t_{0}$ an integration constant. In the static limit, i.e., for a time-
independent electric field, we have

$\lim_{\bm{E}(t)\rightarrow\bm{E}}f_{\bm{R}}(t)=\frac{f^{0}_{\bm{R}}}{\tau}\int_{t_{0}}^{t}dt^{\prime}\,e^{-(t-t^{\prime})\left(\frac{1}{\tau}-ie\bm{E}\cdot\bm{R}/\hbar\right)}=f^{0}_{\bm{R}}\,\frac{1-e^{-\frac{t-t_{0}}{\tau}\left(1-ie\tau\bm{E}\cdot\bm{R}/\hbar\right)}}{1-ie\tau\bm{E}\cdot\bm{R}/\hbar}.$
(S10)

The steady-state solution is thus given by

$\displaystyle f_{\bm{R}}(t)$
$\displaystyle=\frac{f^{0}_{\bm{R}}}{\tau}\int_{-\infty}^{t}dt^{\prime}\,e^{-\frac{t-t^{\prime}}{\tau}}\exp\left[\frac{ie}{\hbar}\int_{t^{\prime}}^{t}du\,\bm{E}(u)\cdot\bm{R}\right]$
(S11)
$\displaystyle=f^{0}_{\bm{R}}\int_{0}^{\infty}ds\,e^{-s}\exp\left[\frac{ie}{\hbar}\int_{t-s\tau}^{t}dt^{\prime}\,\bm{E}(t^{\prime})\cdot\bm{R}\right],$
(S12)

where $s=(t-t^{\prime})/\tau$. The exponential factor can be interpreted as
the integrated momentum shift between two scattering events at $t-s\tau$ and
$t$. Going back to momentum space, we have [31]

$f(\bm{k},t)=\sum_{\bm{R}}f_{\bm{R}}(t)e^{i\bm{k}\cdot\bm{R}}=\int_{0}^{\infty}ds\,e^{-s}f^{0}\left(\bm{k}+\frac{e}{\hbar}\int_{t-s\tau}^{t}dt^{\prime}\,\bm{E}(t^{\prime})\cdot\bm{R}\right).$
(S13)

We now consider the following driving field:

$\bm{E}(t)=\bm{E}_{0}+\bm{E}_{1}(t),\qquad\bm{E}_{1}(t)=\bm{\mathcal{E}}_{1}e^{i\omega
t}+\bm{\mathcal{E}}_{1}^{*}e^{-i\omega
t}=2\,\text{Re}\left(\bm{\mathcal{E}}_{1}e^{i\omega t}\right),$ (S14)

where $E_{0}=|\bm{E}_{0}|$ is large compared to $E_{1}=|\bm{E}_{1}|$. In this
case, the Fourier components of the distribution function become

$\displaystyle f_{\bm{R}}(t)$
$\displaystyle=f^{0}_{\bm{R}}\int_{0}^{\infty}ds\,e^{-\left(1-\frac{ie\tau}{\hbar}\,\bm{E}_{0}\cdot\bm{R}\right)s}\exp\left[\frac{ie}{\hbar}\int_{t-s\tau}^{t}dt^{\prime}\,\bm{E}_{1}(t^{\prime})\cdot\bm{R}\right]$
(S15)
$\displaystyle=f^{0}_{\bm{R}}\int_{0}^{\infty}ds\,e^{-\left(1-\frac{ie\tau}{\hbar}\,\bm{E}_{0}\cdot\bm{R}\right)s}\exp\left[\frac{ie}{\hbar}\left(\frac{\bm{\mathcal{E}}_{1}\cdot\bm{R}}{i\omega}\,e^{i\omega
t}\left(1-e^{-i\omega s\tau}\right)+\text{c.c.}\right)\right].$ (S16)

Up to first order in $E_{1}$, we can expand this as

$\displaystyle f_{\bm{R}}(t)$
$\displaystyle\simeq\frac{f^{0}_{\bm{R}}}{1-ie\tau\bm{E}_{0}\cdot\bm{R}/\hbar}+\frac{ief^{0}_{\bm{R}}}{\hbar}\int_{0}^{\infty}ds\,e^{-\left(1-\frac{ie\tau}{\hbar}\,\bm{E}_{0}\cdot\bm{R}\right)s}\left[\frac{\bm{\mathcal{E}}_{1}\cdot\bm{R}}{i\omega}\,e^{i\omega
t}\left(1-e^{-i\omega s\tau}\right)+\text{c.c.}\right]$ (S17)
$\displaystyle=\frac{f^{0}_{\bm{R}}}{1-ie\tau\bm{E}_{0}\cdot\bm{R}/\hbar}\left[1+\frac{e\bm{\mathcal{E}}_{1}\cdot\bm{R}/\hbar}{\omega-e\bm{E}_{0}\cdot\bm{R}/\hbar-\frac{i}{\tau}}\,e^{i\omega
t}-\frac{e\bm{\mathcal{E}}_{1}^{*}\cdot\bm{R}/\hbar}{\omega+e\bm{E}_{0}\cdot\bm{R}/\hbar+\frac{i}{\tau}}\,e^{-i\omega
t}\right].$ (S18)

Defining the frequency-space Fourier components as

$f_{\bm{R}}(t)=\sum_{m=-\infty}^{\infty}f_{\bm{R},m}\,e^{im\omega t},\qquad
f_{\bm{R},m}(\omega)=\frac{\omega}{2\pi}\int_{0}^{2\pi/\omega}dt\,f_{\bm{R}}(t)\,e^{-im\omega
t},$ (S19)

with $f_{\bm{R},-m}=f_{-\bm{R},m}^{*}$, we have, for example,

$\displaystyle f_{\bm{R},0}$
$\displaystyle=\frac{f^{0}_{\bm{R}}}{1-ie\tau\bm{E}_{0}\cdot\bm{R}/\hbar}\left[1-\frac{2|e\bm{\mathcal{E}}_{1}\cdot\bm{R}/\hbar|^{2}}{\left(\omega-e\bm{E}_{0}\cdot\bm{R}/\hbar-\frac{i}{\tau}\right)\left(\omega+e\bm{E}_{0}\cdot\bm{R}/\hbar+\frac{i}{\tau}\right)}+\mathcal{O}(E_{1}^{4})\right],$
(S20) $\displaystyle f_{\bm{R},1}$
$\displaystyle=\frac{f^{0}_{\bm{R}}}{1-ie\tau\bm{E}_{0}\cdot\bm{R}/\hbar}\left[\frac{e\bm{\mathcal{E}}_{1}\cdot\bm{R}/\hbar}{\omega-e\bm{E}_{0}\cdot\bm{R}/\hbar-\frac{i}{\tau}}+\mathcal{O}(E_{1}^{3})\right],$
(S21) $\displaystyle f_{\bm{R},2}$
$\displaystyle=\frac{f^{0}_{\bm{R}}}{1-ie\tau\bm{E}_{0}\cdot\bm{R}/\hbar}\left[\frac{e\bm{\mathcal{E}}_{1}\cdot\bm{R}/\hbar}{\omega-e\bm{E}_{0}\cdot\bm{R}/\hbar-\frac{i}{\tau}}\frac{e\bm{\mathcal{E}}_{1}\cdot\bm{R}/\hbar}{2\omega-e\bm{E}_{0}\cdot\bm{R}/\hbar-\frac{i}{\tau}}+\mathcal{O}(E_{1}^{4})\right].$
(S22)

##  III Dressed optical conductivity

The steady-state current is given by
$\bm{j}(t)=\bm{j}_{\text{Bloch}}(t)+\bm{j}_{\text{geom}}(t)$ with

$\displaystyle\bm{j}_{\text{Bloch}}(t)$
$\displaystyle=-e\int_{\bm{k}}f_{\bm{k}}(t)\nabla_{\bm{k}}\varepsilon_{\bm{k}}=\frac{ie}{\hbar
V_{c}}\sum_{\bm{R}}\bm{R}\,\varepsilon_{-\bm{R}}f_{\bm{R}}(t),$ (S23)
$\displaystyle\bm{j}_{\text{geom}}(t)$
$\displaystyle=-\frac{e^{2}}{\hbar}\left[\bm{E}(t)\times\hat{z}\right]\int_{\bm{k}}f_{\bm{k}}(t)\Omega_{\bm{k}}=\frac{e^{2}}{\hbar
V_{c}}\,\hat{z}\times\sum_{\bm{R}}\Omega_{-\bm{R}}\bm{E}(t)f_{\bm{R}}(t),$
(S24)

with

$V_{c}\int_{\bm{k}}e^{i\bm{k}\cdot\bm{R}}=\delta_{\bm{R},\bm{0}},\qquad\sum_{\bm{R}}e^{i\bm{k}\cdot\bm{R}}=\frac{(2\pi)^{2}}{V_{c}}\,\delta(\bm{k}).$
(S25)

The frequency components of the currents are thus given by

$\displaystyle\bm{j}_{\text{Bloch}}^{(m)}(\omega)$
$\displaystyle=\frac{ie}{\hbar
V_{c}}\sum_{\bm{R}}\bm{R}\,\varepsilon_{-\bm{R}}f_{\bm{R},m},$ (S26)
$\displaystyle\bm{j}_{\text{geom}}^{(m)}(\omega)$
$\displaystyle=\frac{e^{2}}{\hbar
V_{c}}\,\hat{z}\times\sum_{\bm{R}}\Omega_{-\bm{R}}\left(\bm{E}_{0}f_{\bm{R},m}+\bm{\mathcal{E}}_{1}f_{\bm{R},m-1}+\bm{\mathcal{E}}_{1}^{*}f_{\bm{R},m+1}\right),$
(S27)

with $\bm{j}^{(-m)}=\left(\bm{j}^{(m)}\right)^{*}$. For example, the DC
component of the geometric current becomes

$\displaystyle\bm{j}_{\text{geom}}^{(0)}(\omega)$
$\displaystyle=\frac{e^{2}}{\hbar
V_{c}}\,\hat{z}\times\sum_{\bm{R}}\frac{\Omega_{-\bm{R}}f^{0}_{\bm{R}}}{1-ie\tau\bm{E}_{0}\cdot\bm{R}/\hbar}\Bigg{\\{}\left[1-\frac{2|e\bm{\mathcal{E}}_{1}\cdot\bm{R}/\hbar|^{2}}{\left(\omega-e\bm{E}_{0}\cdot\bm{R}/\hbar-\frac{i}{\tau}\right)\left(\omega+e\bm{E}_{0}\cdot\bm{R}/\hbar+\frac{i}{\tau}\right)}\right]\bm{E}_{0}$
(S28)
$\displaystyle+\left[\frac{\left(e\bm{\mathcal{E}}_{1}\cdot\bm{R}/\hbar\right)\bm{\mathcal{E}}_{1}^{*}}{\omega-e\bm{E}_{0}\cdot\bm{R}/\hbar-\frac{i}{\tau}}-\frac{\left(e\bm{\mathcal{E}}_{1}^{*}\cdot\bm{R}/\hbar\right)\bm{\mathcal{E}}_{1}}{\omega+e\bm{E}_{0}\cdot\bm{R}/\hbar+\frac{i}{\tau}}\right]+\mathcal{O}(E_{1}^{4})\Bigg{\\}}.$
(S29)

In lowest order of $|e\bm{\mathcal{E}}_{1}\cdot\bm{R}/\hbar\omega|$, the first
harmonics are given by

$\displaystyle\bm{j}_{\text{Bloch}}^{(1)}(\omega)$
$\displaystyle=\frac{ie}{\hbar
V_{c}}\sum_{\bm{R}}\frac{\bm{R}\varepsilon_{-\bm{R}}f^{0}_{\bm{R}}}{1-ie\tau\bm{E}_{0}\cdot\bm{R}/\hbar}\frac{e\bm{\mathcal{E}}_{1}\cdot\bm{R}/\hbar}{\omega-e\bm{E}_{0}\cdot\bm{R}/\hbar-\frac{i}{\tau}},$
(S30) $\displaystyle\bm{j}_{\text{geom}}^{(1)}(\omega)$
$\displaystyle=\frac{e^{2}}{\hbar
V_{c}}\,\hat{z}\times\sum_{\bm{R}}\frac{\Omega_{-\bm{R}}f^{0}_{\bm{R}}}{1-ie\tau\bm{E}_{0}\cdot\bm{R}/\hbar}\left[\bm{\mathcal{E}}_{1}+\frac{\left(e\bm{\mathcal{E}}_{1}\cdot\bm{R}/\hbar\right)\bm{E}_{0}}{\omega-e\bm{E}_{0}\cdot\bm{R}/\hbar-\frac{i}{\tau}}\right].$
(S31)

We now define the _dressed_ optical conductivity
$\sigma_{ij}(\omega,\bm{E}_{0})$ through
$j_{i}^{(1)}=\sigma_{ij}\mathcal{E}_{1j}$. The current can thus be written as

$j_{i}(t)=j_{i}^{(0)}+2\,\text{Re}\left[\sigma_{ij}(\omega,\bm{E}_{0})\mathcal{E}_{1j}e^{i\omega
t}\right]+\mathcal{O}(E_{1}^{2}).$ (S32)

Making use of $(\bm{a}\times\bm{b})_{i}=\epsilon_{ijk}a^{j}b^{k}$ with
$\epsilon_{i3k}=-\epsilon_{ik}$ the permutation symbol, we find

$\displaystyle\sigma_{ij}(\omega,\bm{E}_{0})$
$\displaystyle=\frac{ie^{2}}{\hbar^{2}V_{c}}\sum_{\bm{R}}\frac{R_{i}R_{j}\varepsilon_{-\bm{R}}f^{0}_{\bm{R}}}{\left(1-ie\tau\bm{E}_{0}\cdot\bm{R}/\hbar\right)\left(\omega-e\bm{E}_{0}\cdot\bm{R}/\hbar-\frac{i}{\tau}\right)}$
(S33) $\displaystyle-\frac{\epsilon_{ik}e^{2}}{\hbar
V_{c}}\sum_{\bm{R}}\frac{\Omega_{-\bm{R}}f^{0}_{\bm{R}}}{1-ie\tau\bm{E}_{0}\cdot\bm{R}/\hbar}\left[\delta_{kj}+\frac{eE_{0k}R_{j}/\hbar}{\omega-e\bm{E}_{0}\cdot\bm{R}/\hbar-\frac{i}{\tau}}\right],$
(S34)

with $\sigma_{ij}(\omega,\bm{E}_{0})=\sigma_{ij}(-\omega,\bm{E}_{0})^{*}$. As
a check, we undress the conductivity:

$\displaystyle\sigma_{ij}(\omega,\bm{0})$
$\displaystyle=\frac{ie^{2}}{\hbar^{2}V_{c}}\sum_{\bm{R}}\frac{R_{i}R_{j}\varepsilon_{-\bm{R}}f^{0}_{\bm{R}}}{\omega-\frac{i}{\tau}}-\epsilon_{ij}\frac{e^{2}}{\hbar}\frac{1}{V_{c}}\sum_{\bm{R}}\Omega_{-\bm{R}}f^{0}_{\bm{R}}$
(S35)
$\displaystyle=\frac{e^{2}}{\hbar^{2}}\int_{\bm{k}}f^{0}_{\bm{k}}\,\frac{\partial_{i}\partial_{j}\varepsilon_{\bm{k}}}{i\omega+\frac{1}{\tau}}-\epsilon_{ij}\frac{e^{2}}{\hbar}\int_{\bm{k}}f^{0}_{\bm{k}}\Omega_{\bm{k}},$
(S36)

with $\partial_{i}=\partial/\partial k_{i}$. This is a well-known result for
the conductivity: the first term is the Drude contribution and the second term
is the anomalous Hall conductivity.

Let us now consider the dressed longitudinal conductivity
$\sigma_{L}=\left(\sigma_{xx}+\sigma_{yy}\right)/2$ which is a scalar and the
dressed Hall conductivity, $\sigma_{H}=\left(\sigma_{xy}-\sigma_{yx}\right)/2$
which is a pseudoscalar. We find

$\displaystyle\sigma_{L}(\omega,\bm{E}_{0})$
$\displaystyle=\frac{ie^{2}}{2\hbar^{2}V_{c}}\sum_{\bm{R}}\frac{R^{2}\varepsilon_{-\bm{R}}f^{0}_{\bm{R}}}{\left(1-ie\tau\bm{E}_{0}\cdot\bm{R}/\hbar\right)\left(\omega-e\bm{E}_{0}\cdot\bm{R}/\hbar-\frac{i}{\tau}\right)},$
(S37) $\displaystyle\sigma_{H}(\omega,\bm{E}_{0})$
$\displaystyle=\frac{1}{2}\epsilon_{ij}\sigma_{ij}=-\frac{e^{2}}{\hbar}\frac{1}{V_{c}}\sum_{\bm{R}}\frac{\Omega_{-\bm{R}}f^{0}_{\bm{R}}}{1-ie\tau\bm{E}_{0}\cdot\bm{R}/\hbar}\left(1+\frac{1}{2}\frac{e\bm{E}_{0}\cdot\bm{R}/\hbar}{\omega-e\bm{E}_{0}\cdot\bm{R}/\hbar-\frac{i}{\tau}}\right),$
(S38)

where we used $\epsilon_{ij}\epsilon_{ik}=\delta_{jk}$. When the static field
is strong, i.e, for $\omega_{B}\tau\gg 1$, where $\omega_{B}=eE_{0}L/\hbar$ is
the Bloch frequency with $L$ the lattice constant, $\sigma_{L}$ and
$\sigma_{H}$ simplify to

$\displaystyle\sigma_{L}(\omega,\bm{E}_{0})$
$\displaystyle=-\frac{e^{2}}{h}\frac{\pi}{\tau\hbar
V_{c}}\sum_{\bm{R}}\frac{R^{2}\varepsilon_{-\bm{R}}f^{0}_{\bm{R}}}{\omega_{\bm{R}}\left(\omega-\omega_{\bm{R}}-\frac{i}{\tau}\right)}=\frac{e^{2}}{h}\frac{\pi}{V_{c}}\sum_{\bm{R}}\left(\frac{-R^{2}\varepsilon_{-\bm{R}}f^{0}_{\bm{R}}}{\hbar\omega_{\bm{R}}}\right)\frac{\left(\omega-\omega_{\bm{R}}\right)\tau+i}{\left(\omega-\omega_{\bm{R}}\right)^{2}\tau^{2}+1},$
(S39) $\displaystyle\sigma_{H}(\omega,\bm{E}_{0})$
$\displaystyle=\frac{e^{2}}{h}\frac{\pi}{i\tau
V_{c}}\sum_{\bm{R}}\frac{\Omega_{-\bm{R}}f^{0}_{\bm{R}}}{\omega-\omega_{\bm{R}}-\frac{i}{\tau}}=\frac{e^{2}}{h}\frac{\pi}{V_{c}}\sum_{\bm{R}}\left(-i\Omega_{-\bm{R}}f^{0}_{\bm{R}}\right)\frac{\left(\omega-\omega_{\bm{R}}\right)\tau+i}{\left(\omega-\omega_{\bm{R}}\right)^{2}\tau^{2}+1},$
(S40)

with $\omega_{\bm{R}}=e\bm{E}_{0}\cdot\bm{R}/\hbar$.