# Graphlets in Multiplex Networks

Tamara Dimitrova  Macedonian Academy of Sciences and Arts, Skopje, Republic of
Macedonia  IKT-Labs, Skopje, Macedonia  These authors contributed equally to
this work  Kristijan Petrovski  Macedonian Academy of Sciences and Arts,
Skopje, Republic of Macedonia  IKT-Labs, Skopje, Macedonia  These authors
contributed equally to this work  Ljupco Kocarev  Macedonian Academy of
Sciences and Arts, Skopje, Republic of Macedonia  Faculty of Computer Science
and Engineering, UKIM, Skopje, Republic of Macedonia  IKT-Labs, Skopje,
Macedonia<EMAIL_ADDRESS>
###### Abstract

We develop graphlet analysis for multiplex networks and discuss how this
analysis can be extended to multilayer and multilevel networks as well as to
graphs with node and/or link categorical attributes. The analysis has been
adapted for two typical examples of multiplexes: economic trade data
represented as a 957-plex network and 75 social networks each represented as a
12-plex network. We show that wedges (open triads) occur more often in
economic trade networks than in social networks, indicating the tendency of a
country to produce/trade of a product in local structure of triads which are
not closed. Moreover, our analysis provides evidence that the countries with
small diversity tend to form correlated triangles. Wedges also appear in the
social networks, however the dominant graphlets in social networks are
triangles (closed triads). If a multiplex structure indicates a strong tie,
the graphlet analysis provides another evidence for the concepts of
strong/weak ties and structural holes. In contrast to Granovetter’s seminal
work on the strength of weak ties, in which it has been documented that the
wedges with only strong ties are absent, here we show that for the analyzed 75
social networks, the wedges with only strong ties are not only present but
also significantly correlated.

###### keywords:

Graphlets, Multiplex Networks, Social Networks, Economy Trade Networks

## Introduction

Networks are ubiquitous. This has triggered a vast amount of research in the
last two decades. When focusing on the local connectivity of subgraphs within
a network, two approaches have been identified, motifs and graphlets. Motifs
are defined as sub-graphs that repeat at frequency higher than in the random
graphs [1, 2]. However, they depend on the choice of the network’s null model.
In contrast, graphlets [3] are induced sub-graphs of a network that appear at
any frequency and hence are independent of the null model. Graphlets have
found numerous applications as building blocks of network analysis in various
disciplines ranging from social science [4, 5] to biology [6, 7]. In social
science, graphlet analysis (known as sub-graph census) is widely adopted in
sociometric studies [4]. Recently, graphlet analysis has also been adapted for
directed networks [8, 9].

Real networks come in various structures including multilayer networks,
multiplex networks, multilevel networks or multiscale networks. Understanding
local connectivity of subgraphs by using graphlets in these networks could,
therefore, lead to improved predictive modes of networks and/or to enhanced
description of their properties. In sociology, the importance of multiplex
networks has been emphasized by many scholars. White, Boorman and Breiger [10]
and Boorman and White [11] treated multiple networks as a foundation of social
structure and argued that the patterning and interweaving of different types
of ties are needed to describe and characterize social structures.
Multiplexity is critical to diverse phenomena such as the mobilization of
social movements [12], the consolidation of political power [13], the
emergence of trust in economic relationships [14], the creation of social
bonds within civic networks [15], and the organization of party coalitions
[16]. Multiplex networks have been studied to understand scientific
collaboration [17], structural logic of intra-organizational networks [18],
formation of ties featuring both an economic and a social component in inter-
organizational networks [19], and formation of relationships among producers
in the multiplex triads [20].

Multilayer/multiplex networks [21, 22] have recently been a subject of
particularly intense research by the network science and physics communities.
Novel structural descriptors [23, 24, 25, 26, 27] and tools from statistical
physics [28, 29] have been developed for studying multilayer networks. By
analyzing multilayer networks, instead of relying on their monolayer
counterparts, scientists have documented evidences of novel features and novel
insights about real systems (see, for example, [30, 31, 32]). Multilayer
networks are fundamental for understanding of dynamical processes on networked
systems, including, for example, spreading processes, such as flows (and
congestion) in transportation networks [33, 34] , and information and disease
spreading in social networks [35, 36, 37, 38, 39].

This paper aims at developing graphlet analysis for multiplex networks. The
study has been motivated by two facts: (1) graphlets are powerful tool for
analyzing single plex/layer networks, and (2) networks with
multiplex/multilayer structures are common in nature and societies. The
central problem in addressing graphlets in a multiplex network comes from the
explosive growth in the number of various edge types with the linear increase
in the number of plexes in the network. We now present an illustrative example
showing graphlet analysis’s challenges and how they are addressed for
multiplex networks.

Figure 1:  Upper-left panel: single plex graph. Nodes of a graphlet are
classified into different orbits. Upper-right panel: multiplex graph with 3
plexes $a,b,$ and $c$. On the same panel the multiplex graph is represented as
edge-labeled graph: each edge is labeled with a single label, an element of
the set $\\{a,b,c,ab,ac,abc\\}$. Left table: triangles and wedge stars. Right
table: wedge paths. The column ‘nodes’ shows the nodes of the corresponding
orbit (triangle, wedge star, or wedge path). The node in the small bracket
(for the first triangle this is the node 1 shown as (1)) indicates the first
node of the orbit. The column ‘edges’ shows the full description of sub-
orbits. While for the single plex graph (upper-left panel) all triangles are
same, for the multiplex graph (upper-right pane) these triangles are different
and are labeled with a single label from the set $E_{t}$. For example the
edges of first triangle are labeled with $ab.ab.b$. The column ‘edges_1’ shows
the first reduction of the number of sub-orbits. Thus, for example, all wedge
stars $ab.ac$, $ab.ab$, $ab.bc$, and $ac.bc$ are represented as $2.2$. The
column ‘edges_2’ represents another reduction in which $ab.ac$, $ab.bc$, and
$ac.bc$ are represented as $2_{x}.2_{y}$, while $ab.ab$ is represented as
$2_{x}.2_{x}$. The number of wedge paths are doubled with paths in which the
other non-central node is a starting node of the path.

Consider a simple undirected (single-plex) graph $G=(V,E)$, see Fig. 1, for
which

$\displaystyle V$ $\displaystyle=$ $\displaystyle\\{1,2,3,4,5,6,7,8\\}$ (1)
$\displaystyle E$ $\displaystyle=$
$\displaystyle\\{12,13,23,28,34,36,38,45,67,68,78\\}.$ (2)

A graphlet $G^{\prime}=(V^{\prime},E^{\prime})$ is an induced subgraph of $G$
(see Materials and Methods). Thus, for example,
$G^{\prime}=(\\{3,4,5\\},\\{34,45\\})$ is a graplet called wedge. By taking
into account the “symmetries” between nodes in a graphlet, the nodes of the
graphlet can be classified into different orbits (see Materials and Methods).
The node 3 of the graphlet $G^{\prime}=(\\{3,4,5\\},\\{34,45\\})$ belongs to
the orbit 1 called “wedge star”, while the node 4 of the graphlet $G^{\prime}$
is classified as “wedge path” and is denoted as orbit 2. The orbit 0 of a
graphlet is node degree. Consider now a 3-plex network shown on Figure 1. We
write $a,b,c$ for network plexes. Let $E_{t}=\\{a,b,c,ab,ac,bc,abc\\}$. A
multiplex graph is a special type of edge-labeled graph in which each edge is
labeled only with a single label – an element of the set $E_{t}$. Thus, for
our example, the multiplex graph can be written as labeled $G=(V,E_{label})$,
see Fig. 1, for which

$\displaystyle V$ $\displaystyle=$ $\displaystyle\\{1,2,3,4,5,6,7,8\\}$ (3)
$\displaystyle E_{label}$ $\displaystyle=$
$\displaystyle\\{12_{ab},13_{ab},23_{b},28_{b},34_{ac},36_{ab},38_{bc},45_{b},67_{abc},68_{c},78_{b}\\},$
(4)

where $e_{\alpha}$ means the edge $e\in E$ is labeled with $\alpha\in E_{t}$ .
Two problems will be addressed when dealing with graphlet analysis of
multiplex networks. The first problem is related to the fact that each orbit
of a multiplex network consists of sub-orbits. Thus, the first orbit, the
orbit 0, consists of 7 sub-orbits: $0_{a}$, $0_{b}$, $0_{c}$, $0_{ab}$,
$0_{ac}$, $0_{bc}$, and $0_{abc}$. To successfully define the sub-orbits for
other orbits (wedges, triangles, and so on) we first introduce two notations:
(1) permutation $[a,b,c]$, stressing the order of edge types as part of the
graphlet and (2) set $\\{a,b,c\\}$ stressing the symmetry the edges have as
part of the graphlet. Thus, for example, the sub-orbit of the orbit wedge star
for the node 4 in the graphlet with vertex set $\\{3,4,5\\}$ is $2_{ac.b}$
which is for simplicity on the Figure 1 written as $ac.b$. For our example,
all sub-orbits of the orbit 2 (wedge star) are elements of the set
$2_{\\{\alpha.\beta\\}}$ where $\alpha,\beta\in E_{t}$; the cardinality of
this set (the set of wedge star sub-orbits) is 28. The sub-orbit of the orbit
wedge path for the node 2 in the graphlet with vertex set $\\{2,8,6\\}$ is
$1_{b.c}$ or simply $b.c$. For our example, all sub-orbits of the orbit 1
(wedge path) are elements of the set $2_{[\alpha.\beta]}$ where
$\alpha,\beta\in E_{t}$; the total number of wedge star sub-orbits is 49.
Section Materials and Methods describes how sub-orbits are defined for a given
orbit and the size of the orbit (the cardinality of the set of sub-orbits
associated with the orbit).

The second problem is the size of the orbit. Even for the simplest graphlet, a
node degree, in a multiplex network the size of orbit 0 (node degree) equals
$2^{d}-1$, where $d$ is the number of plexes. This limits the application of
graphlets for real data even for networks with small number of plexes. In
order to address this problem, we propose two different methods for reducing
the number of graphlets. Both methods are illustrated on Figure 1. The
simplest possible reduction of the $E_{t}=\\{a,b,c,ab,ac,cb,abc\\}$ is to the
set $E_{t}^{{}^{\prime}}=\\{1,2,3\\}$. In another words, an edge of a
multiplex, which has been defined as an element of the set $E_{t}$
representing different relation types, is now (after reduction) defined as an
element of the set $E_{t}^{{}^{\prime}}$ representing the ‘strength’ (plex
count) of the original edge. Thus wedge stars for nodes 3,1,4 and 3,1,6 which
represent two different sub-orbits $ab.ac$ and $ab.ab$, respectively, are now
represented as a single sub-orbit $2.2$. More generally, in this example, all
sub-orbits $\alpha.\beta$ such that $\alpha,\beta\in E_{t}$ and
$|\alpha|=|\beta|=2$ are represented as a single orbit $2.2$ (see Figure 1).
The second reduction takes into account that the sub-orbit $2.2$ is result of
reduction of two subsets:

$\displaystyle A$ $\displaystyle=$
$\displaystyle\left\\{\alpha.\beta:\alpha\neq\beta,\alpha,\beta\in
E_{t},|\alpha|=|\beta|=2\right\\}$ (5) $\displaystyle B$ $\displaystyle=$
$\displaystyle\left\\{\alpha.\alpha:\alpha\in E_{t},|\alpha|=2\right\\}.$ (6)

Both subsets $A$ and $B$, in this reduction, are mapped to two different sub-
orbits $2_{x}.2_{y}$ and $2_{x}.2_{x}$, respectively (see Figure 1). This
reduction results in a multiplex which is called plex-count multiplex with
distinct links inside orbits. Section Materials and Methods provides the full
description of both reductions.

The graphlet analysis for multiplex networks can be extended to graphs with
node and/or link (categorical) attributes as well as to multilayer and/or
multilevel networks as discussed in the supporting information (SI). Further,
economic and social multiplex networks are analyzed using graphlets providing
novel insight of networks’ properties and/or their local structures. In
particular, for economic networks we show that (1) countries produce/trade
products in local structure of triads which are not closed and (2) counties
with small diversity tend to form correlated triangles. For social networks
for which a strong tie is related to the multiplex structure, we provide an
example of social networks for which the wedges with only strong ties are both
present and significantly correlated, in contrast to the Granovetter’s seminal
work on the strength of weak ties, in which it has been shown that the wedges
with only strong ties are absent.

## Materials and Methods

### Multilayer/multiplex networks and graphlets

This paper introduces a method for graphlet analysis of multiplex networks.
This analysis can also be adopted for graphs with node and/or edge
(categorical) attributes. In the SI we explain how the method can be
generalized to graphs with node and/or edge attribute and complex (multiplex,
multilayer, and multilevel) graphs. Next, definitions of multiplex, multilayer
and multilevel networks as well as graphlets for simple graphs are provided.

Multiplex networks – A multiplex network is defined as a $d+1$-tuple
$G=(V,E^{1},\ldots,E^{d})$ where $V$ is the set of nodes and for each
$\alpha\in\\{1,2,\ldots,d\\}$, $E^{\alpha}$ is the set of edges describing the
presence or absence of edges of type $\alpha$ between pairs of nodes. Since a
multiplex network is uniquely defined with the node set $V$ and the edge sets
$E^{1},\ldots,E^{d}$, we write $G(V,E^{1},\ldots,E^{d})$ to denote the
multiplex. The graph $(V,E^{\alpha})$ is also called plex. A $k$-plex network
is a subgraph of the multiplex and is defined as
$G=(V,E^{\alpha_{1}},\ldots,E^{\alpha_{k}})$.

Multilayer networks – A multilayer network is defined as a graph $G=(V,E)$ for
which $V\subseteq V^{1}\times V^{2}\ldots\times V^{d}$ and $E$ is the set of
edges. Typically, $V^{\alpha}=V$, $\alpha=1,\ldots,d$; the elements of $V$ are
called nodes. A layer is a sub-graph $(V^{\alpha},E^{\alpha})$ for which
$E^{\alpha}=\\{ij\in E:i,j\in V^{\alpha}\\}$. Let $A=[a_{ij}^{\alpha\beta}]$,
$i,j=1,\ldots n$ and $\alpha,\beta=1\ldots m$, be the adjacency matrix of the
graph $G=(V,E)$. In general it is assumed that $a_{ij}^{\alpha\beta}\neq 0$
for $\alpha\neq\beta$ and $i\neq j$. In the special case when a network is
represented with an adjacency matrix such that $a_{ij}^{\alpha\beta}=0$ when
both $\alpha\neq\beta$ and $i\neq j$, the network is also called multiplex
network. This definition implies that a multiplex consists of layers (not
plexes). In what follows, we will not consider multilayer approach of
multiplex networks.

Multilevel networks – A multilevel (interconnected) network is a graph
$G=(V,E)$ for which $V=\cup^{d}_{\alpha=1}V^{\alpha}$ and $E$ is the set of
edges. In general, $V^{\alpha}\neq V^{\beta}$ for $\alpha\neq\beta$; each
$V^{\alpha}$ represents a distinct type of nodes. Multilevel networks can be
described with subgraphs $(V^{\alpha},E^{\alpha})$ called levels, where
$E^{\alpha}=\\{ij\in E:i,j\in V^{\alpha}\\}$ and bipartite subgraphs
$(V^{\alpha},V^{\beta},E^{\alpha\beta})$ such that $E^{\alpha\beta}=\\{ij\in
E:i\in V^{\alpha},j\in V^{\beta}\\}$.

Graphs with attributes – Let $G=(V,E)$ be a simple graph, $V_{A}$ be the set
of node attributes, and $E_{A}$ be the set of edge attributes. Nodes are
labeled $i=1,2,\ldots,|V|$, node attributes are labeled
$q=1,2,\ldots,|V_{A}|$, while edge attributes are labeled
$\alpha=1,2,\ldots,|E_{A}|$. We define $|V|\times|V_{A}|$ matrix $D=[d_{iq}]$
as $d_{iq}=1$ if and only if the node $i$ has the attribute $q$, otherwise 0.
We define a sub-graph $G^{\alpha}=(V,E^{\alpha})$ as $E^{\alpha}=\\{e|e\in
E\mbox{ and }e\mbox{ has attribute }\alpha\\}$.

Graphlets – Let $G=(V,E)$ be a graph, where $V$ is a set of nodes and $E$ is a
set of edges. A subgraph $G^{\prime}$ of $G$ is a graph whose set of nodes and
set of edges are subsets of $G$. An induced subgraph of $G$,
$G^{\prime}=(V^{\prime},E^{\prime})$, is a subgraph that consists of a subset
of nodes in $G$ and all of the edges that connect them in $G$, i.e.
$V^{\prime}\subset V$, $E^{\prime}=\\{(u,v):u,v\in V\\},(u,v)\in E^{\prime}$.
The size of a graphlet is the cardinality of its node set. Unless we
explicitly say “induced” in this paper, a subgraph is not necessarily induced.
The top part of Table 1 shows all 2-, 3-, and 4-node undirected graphlets
$G_{i}$, 0 $\leq i\leq 8$. By taking into account the “symmetries” between
nodes in the graphlet $G_{i}$, the nodes of $G_{i}$ are classified into
different automorphism orbits (or just orbits, for brevity), where the nodes
with the same orbit identification are topologically identical. For all
$G_{i}$, $0\leq i\leq 8$, there are 15 orbits, which are also shown as filled
nodes in the top panel of Table 1.

### Graphlets in multiplex graphs

Table 1 depicts all 2-4 node graphlets with their orbits (15 in total) for a
simple graph $G=(V,E)$. For multiplex graphs, these orbits are further
subdivided based on the relation (edge) types, so that each orbit consists of
a number of sub-orbits that are defined by the specific edge types inside the
graphlet. Let ${\cal R}=\\{1,2,\ldots,d\\}$ be the set of relation/edge types.
Let ${\cal P}({\cal R})$ be the power set of the set ${\cal R}$, that is, the
set of all subsets of ${\cal R}$, we define $E_{t}={\cal P}({\cal
R})\setminus\emptyset$. The set $E_{t}$ is the set of combinations of
different edge types with cardinality $|E_{t}|=2^{d}-1$. Multiplex graphs are
treated as having $2^{d}-1$ different types of edges so that each edge of a
multiplex is uniquely defined with an element of the set $E_{t}$. Thus, for
example, for ${\cal R}=\\{1,2,3\\}$, an edge of the multiplex is represented
by an element of the set $\\{1,2,3,12,13,23,123\\}$. To simplify the notation,
the set $\\{a,b,c\\}$ \- the element of the set $E_{t}$ \- will be labeled as
$abc$. Thus, the first orbit – the orbit 0 – see Table 1, consists of 7 sub-
orbits: $0_{1}$, $0_{2}$, $0_{3}$, $0_{12}$, $0_{13}$, $0_{23}$, and
$0_{123}$. To successfully define the orbit class for each orbit we first
define/introduce two notations: (1) permutation $[a,b,c]$, highlighting the
order of the types of edges as part of a graphlet and (2) set $\\{a,b,c\\}$
(can have repeating elements, only their order does not matter), highlighting
the symmetry the edges have as part of a graphlet. Thus for example, in a
multiplex the orbit 1 consists of the orbit class $1_{[a,b]},a,b\in E_{t}$
while the orbit 2 consists of the orbit class $2_{\\{a,b\\}},a,b\in E_{t}$.
For ${\cal R}=\\{1,2\\}$ the sub-orbit instances of orbit 1 are: $1_{1.1}$,
$1_{1.2}$, $1_{1.12}$, $1_{2.1}$, $1_{2.2}$, $1_{2.12}$, $1_{12.1}$,
$1_{12.2}$, and $1_{12.12}$, while the sub-orbit instances of orbit 2 are:
$2_{1.1}$, $2_{1.2}$, $2_{1.12}$, $2_{2.2}$, $2_{2.12}$, and $2_{12.12}$.
Lexical sorting (smaller plex count links first, lexical sort on same plex
count links) is used to order the edges inside every wedge star sub-orbit. In
this way every sub-orbit has an unique representation with which it is further
identified.

A straightforward way to enumerate all possible graphlet sub-orbits is
described in the supplementary material. Table 1 visually shows the first 15
orbits and gives information about their multiplex classes. Each orbit is
uniquely represented with an orbit class, labeled as shown in Table 1 and the
cardinality of this class is called the size of the orbit. To simplify the
notation we write $e^{O_{i}}_{1}.e^{O_{i}}_{2}.\cdots e^{O_{i}}_{k}$ for a
sub-orbit of the given orbit $O_{i}$. When it is clear from the text what
orbit is considered, we will drop the explicit $O_{i}$ in order to keep the
notation uncluttered. Each orbit $O_{i}$ can then viewed as a set of sub-
orbits, $O_{i}=\\{e_{1}.e_{2}.e_{3}..e_{k}\hskip 2.84544pt|\hskip
2.84544pte_{j}\in E_{t}\\}$. Inside the orbit every list of edges
$e_{1}..e_{k}$ is connected in the same way, however they are differentiated
by the specific values of the multiplex links $e_{j}$ coming from the set
$E_{t}$, explained visually in 2.

Figure 2: Sub-orbit breakdown for graphlets up to order 4. The top panel shows
all multiplex $G_{k}$ graphlets, $0\leq k\leq 8$ and their 15 orbits. The
orbit nodes are colored gray. Every edge is identified by a letter. Each blue
colored edge belongs to a permutation group ([]) inside the orbit class, while
each orange edge belongs to a symmetric set group ({}). Edges that belong to
symmetric sets are always lexically sorted (smaller plex count links first,
lexical sort on same size links) before assigning them a specific sub-orbit
instance.

$\begin{array}[H]{lllll}\mbox{Orbit}&\mbox{Orbit Alias}&\mbox{Orbit
class}&&\mbox{Size of orbit}\\\ \hline\cr\\\
0&\mbox{Degree}&{0_{[a]}}&{a\in{E_{t}}}&{|E_{t}|}\\\ 1&\mbox{Wedge
Path}&{1_{[a,b]}}&{a,b\in{E_{t}}}&{|E_{t}|^{2}}\\\ 2&\mbox{Wedge
Star}&{2_{\\{a,b\\}}}&{a,b\in{E_{t}}}&\left({{{|E_{t}|}\choose{2}}}\right)\\\
3&\mbox{Triangle}&{3_{[\\{a,b\\},c]}}&{a,b,c\in{E_{t}}}&{\left({{{|E_{t}|}\choose{2}}}\right)\cdot|E_{t}|}\\\
4&\mbox{Four Path}&{4_{[a,b,c]}}&{a,b,c\in{E_{t}}}&{|E_{t}|^{3}}\\\
5&&{5_{{}[a,[b,c]]}}&{a,b,c\in{E_{t}}}&{|E_{t}|^{3}}\\\ 6&\mbox{Three
Star}&{6_{[a,\\{b,c\\}]}}&{a,b,c\in{E_{t}}}&{|E_{t}|\cdot\left({{{|E_{t}|}\choose{2}}}\right)}\\\
7&&{7_{\\{a,b,c\\}}}&{a,b,c\in{E_{t}}}&{\left({{{|E_{t}|}\choose{3}}}\right)}\\\
8&\mbox{Four
Cycle}&{8_{\\{[a,d],[b,c]\\}}}&{a,b,c\in{E_{t}}}&\left({{{|E_{t}|^{2}}\choose{2}}}\right)\\\
9&\mbox{Tailed
Triangle}&{9_{[a,\\{b,c\\},d]}}&{a,b,c,d\in{E_{t}}}&{|E_{t}|\cdot{{|E_{t}|}\choose{2}}\cdot|E_{t}|}\\\
10&&{10_{[a,\\{b,c\\},d]}}&{a,b,c,d\in{E_{t}}}&{|E_{t}|\cdot{{|E_{t}|}\choose{2}}\cdot|E_{t}|}\\\
11&&{11_{[a,b,c,d]}}&{a,b,c,d\in{E_{t}}}&{|E_{t}|^{4}}\\\ 12&\mbox{Chordal
Cycle}&{12_{[\\{[a,c],[b,d]\\},e]}}&{a,b,c,d\in{E_{t}}}&\left({{{|E_{t}|^{2}}\choose{2}}}\right)\cdot|E_{t}|\\\
13&&{13_{[\\{[a,b],[c,d]\\},e]}}&{a,b,c,d\in{E_{t}}}&\left({{{|E_{t}|^{2}}\choose{2}}}\right)\cdot|E_{t}|\\\
14&\mbox{Four
Clique}&{14_{\\{a,b,c,d,e,f\\}}}&{a,b,c,d,e,f\in{E_{t}}}&{|E_{t}|^{5}}\\\
\end{array}$

Table 1: Sub-orbit breakdown for graphlets up to order 4. Edges that belong to
symmetric sets are always lexically sorted (smaller plex count links first,
lexical sort on same size links) before assigning them a specific sub-orbit
instance. For a set with $n$ elements, the number of $k$-combination with
repetitions is denoted by $\left({{{n}\choose{k}}}\right)$ where
$\left({{{n}\choose{k}}}\right)={{{n+k-1}\choose{k}}}$.

Figure 3: Reducing the number of sub-orbits: plex-count multiplex with
distinct links inside orbits. The first 4 orbits and 36 sub-orbits of a 2-plex
network with ${\cal R}=\\{a,b\\}$. Each plex is identified by a different
color. Each of the nodes with dark gray color is identified with the specific
second reduction sub-orbit.

### Reducing the number of sub-orbits

A multiplex graph is a special type of labeled graph in which each edge is
labeled only with a single label - an element of the set $E_{t}$. The number
of sub-orbits grows exponentially with the linear increase in the number of
plexes in the network. Even for the simplest orbit 0 – the degree – the number
of sub-orbits grows as $2^{d}$. This limits the application of graphlets for
real data even for networks with small number of plexes. In order to address
this problem, we propose two different ways for reducing the number of sub-
orbits.

Let us first consider an example for which ${\cal R}=\\{a,b,c\\}$, so that an
edge of the multiplex is represented by an element of the set
$E_{t}=\\{a,b,c,ab,ac,cb,abc\\}$. In this case, there are 49 wedge-path sub-
orbits. We write $e_{1}.e_{2}$ for a wedge-path sub-orbit, where
$e_{1},e_{2}\in E_{t}$. Thus, for example, $ab.abc$ denotes a sub-orbit with
relation types $a$ and $b$ associated with the first edge and $a,b$ and $c$
with the second edge. The set $O_{wp}$ then contains 49 wedge-path sub-orbits
$O_{wp}=\\{a.a,b.b,c.c,\cdots,abc.ab,abc.ac,abc.bc,abc.abc\\}$. The set of
triangle sub-orbits is similarly large $|O_{tri}|=196$.

For a given orbit, the size of the set of sub-orbits (corresponding to this
orbit) grows exponentially with the linear increase in the number of plexes in
the networks (as shown in the Table 1). In fact, the number of sub-orbits of a
given orbit is polynomial to $|E_{t}|=2^{d}-1$. Therefore, in order to reduce
the size of the set of sub-orbits, one needs to reduce the set $E_{t}$.

The simplest possible reduction of the $E_{t}$ is to the set
$E_{t}^{{}^{\prime}}$ defined as $E_{t}^{{}^{\prime}}=\\{|e|\hskip
2.84544pt|\hskip 2.84544pte\in E_{t}\\}$. Thus, in the above example, the set
$E_{t}=\\{a,b,c,ab,ac,cb,abc\\}$ reduces to the set
$E_{t}^{{}^{\prime}}=\\{1,2,3\\}$. In another words, an edge of a multiplex,
which has been defined as an element of the set $E_{t}$ representing different
relation types, is now (after reduction) defined as an element of the set
$E_{t}^{{}^{\prime}}$ representing the ‘strength’ (plex count) of the original
edge. This reduction leaves us with $d$ possible edge types. The reduced
multiplex is called plex-count multiplex.

For the above example, the reduced wedge-path sub-orbit set is now
$O_{wp}^{{}^{\prime}}=\\{1.1,1.2,1.3,2.1,2.2,2.3,3.1,3.2,3.3\\}$, with
$|O_{wp}^{{}^{\prime}}|=9$, while the triangle sub-orbit set has size
$|O_{tri}^{{}^{\prime}}|=18$. Furthermore, every orbit set has a cardinality
that is $poly(d)$ instead of $poly(2^{d}-1)$. In general, the orbit set is now
defined as
$O_{i}^{{}^{\prime}}=\\{e_{1}^{{}^{\prime}}.e_{2}^{{}^{\prime}}.e_{3}^{{}^{\prime}}..e_{k}^{{}^{\prime}}\hskip
2.84544pt|\hskip 2.84544pte_{j}^{{}^{\prime}}\in E_{t}^{{}^{\prime}}\\}$ where
$e_{i}^{{}^{\prime}}=|e_{i}|$.

The reduction of $E_{t}$ equalizes every $e_{i}$ that has the same plex count.
This greatly reduces the information in the transformed graph. To maintain
some of the original information, we can separate two different links inside a
particular orbit that have the same plex count. Suppose we have two links
$e_{k}=\mathbf{ab},\hskip 2.84544pte_{v}=\mathbf{bc}$, that belong to the
following triangle edge list $\mathbf{ab}.\mathbf{bc}.abc\in O_{tri}$. Clearly
$|e_{k}|=|e_{v}|=2$ however $e_{k}\neq e_{v}$. With the first reduction both
of these links are transformed to $2$, so that $ab.bc.abc\to 2.2.3$. However,
in order to retain the fact that $e_{k}$ and $e_{v}$ are different, we
introduce a new reduction rule $e_{i}^{{}^{\prime\prime}}=d\times
I_{i}+e_{u}^{{}^{\prime}}$, where ${I}_{i}$ represents a distinct id of the
original link $e_{i}$. We assign a different ${I}_{i}$ for every distinct link
that shares the same plex count with another link inside the orbit. With this
mapping, the triangle $ab.bc.abc$ corresponds to $2.5.3$, where
$e_{1}^{{}^{\prime\prime}}=3\times 0+2=2$ and
$e_{2}^{{}^{\prime\prime}}=3\times 1+2=5$, since $d=3$ and
$e_{u}^{{}^{\prime}}=|ab|=|bc|=2$. Specifically for a given orbit list
$e_{1}.e_{2}.\cdots.e_{l}$, the distinct links (of the same plex size) are
given different indices from left to right, starting with 0. This reduction of
the set $E_{t}$ adds some new labels in the set $E_{t}^{{}^{\prime}}$
resulting in the set $E_{t}^{{}^{\prime\prime}}$. For example, the set of sub-
orbits to the wedge-path orbit is now:
$O_{wp}^{{}^{\prime\prime}}=\\{1.1,\mathbf{1.4},1.2,1.3,2.1,2.2,\mathbf{2.5},2.3,3.3\\}$.
Although the size of the set $O_{wp}^{{}^{\prime\prime}}$ is increased
compared with the set $O_{i}^{{}^{\prime}}$, the cardinality of the set
$|O_{i}^{{}^{\prime\prime}}|$ is still polynomial to $d$. We have
$|O_{wp}^{{}^{\prime\prime}}|=11$, and $|O_{tri}^{{}^{\prime\prime}}|=34$ for
$d>2$ and so on. This reduction results in a multiplex which will be called
plex-count multiplex with distinct links inside orbits. Figure 3 shows all
possible sub-orbits of a plex-count multiplex with distinct links inside
orbits with two plexes. The exact sizes of the orbit sets, for both
reductions, are shown in the tables 2 and 3.

$\begin{array}[H]{llll}\mbox{Orbit}&\mbox{Orbit Alias}&\mbox{Orbit
class}&\mbox{Size of orbit}\\\ \hline\cr\\\
0&\mbox{Degree}&{0_{[a_{i}]}}&{d}\\\ 1&\mbox{Wedge
Path}&{1_{[a_{i},a_{j}]}}&{d^{2}}\\\ 2&\mbox{Wedge
Star}&{2_{\\{a_{i},a_{j}\\}}}&{({{d}\choose{2}})}\\\
3&\mbox{Triangle}&{3_{[\\{a_{i},a_{j}\\},a_{k}]}}&{({{d}\choose{2}})*d}\end{array}$
Table 2: Reducing the number of sub-orbits: the size of the orbits 0,1,2, and
3 for a plex-count multiplex network. $a_{i},a_{j},a_{k}\in
E_{t}^{{}^{\prime}}$ represent arbitrary plex-count link.
$\begin{array}[H]{llll}\mbox{Orbit}&\mbox{Orbit Alias}&\mbox{Orbit
class}&\mbox{Size of orbit}\\\ \hline\cr\\\
0&\mbox{Degree}&{0_{[a_{x}]}}&{d}\\\ 1&\mbox{Wedge
Path}&{1_{[a_{x},b_{y}]}}&{d^{2}+(d-1)}\\\ 2&\mbox{Wedge
Star}&{2_{\\{a_{x},b_{y}\\}}}&{({{d}\choose{2}})+(d-1)}\\\
3&\mbox{Triangle}&{3_{[\\{a_{x},b_{y}\\},c_{z}]}}&{\left[({{d}\choose{2}})+(d-1)\right]*d+[(d+1)*(d-1)]}\\\
3*&\mbox{Triangle for
$d>2$}&{3_{[\\{a_{x},b_{y}\\},c_{z}]}}&{\left[({{d}\choose{2}})+(d-1)\right]*d+[(d+1)*(d-1)]+(d-1)}\end{array}$
Table 3: Reducing the number of sub-orbits: the size of the orbits 0, 1, 2,
and 3 for a plex-count multiplex with distinct links inside orbits. $a,b,c\in
E_{t}^{{}^{\prime}}$ represent arbitrary ($a$, $b$ or $c$)-plex count links
while $x,y,z$ are indexes that distinctly identify links of the same plex
size: for two links $a_{x}$ and $b_{y}$, if $a=b$ then $x=y$.

### Graphlet metrics

Here several graphlet metrics for a multiplex network are defined. Given a
node, its graphlet degree of a sub-orbit is the number of times the node is
touched by the sub-orbit. The size of a graphlet is the cardinality of its
node set. Let $G=(V,E^{1},\ldots,E^{D})$ be a multiplex network. We define a
$(n,k)-$signature vector of the node $i$, $SI_{i}(n,k)$, as a vector of
graphlet degrees of the node’s (lexicographically ordered) sub-orbits, for all
graphlets up to the size $n$ of the $k-$plex
$G=(V,E^{\alpha_{1}},\ldots,E^{\alpha_{k}})$.

Let $G$ be a multiplex network with $N$ nodes and $D$ plexes. For large $D$,
computing the signature vector of a node is impossible and we restrict
ourselves, in this case, to sub-multiplexes consisting of $k$ plexes. We first
define a $k$-plex and then construct the vertex $(n,k)-$signature vector of
the $k-$plex for graphlets up to size $n$. In this way, for each vertex and a
given $k-$plex $G(V,E^{\alpha_{1}},\ldots,E^{\alpha_{k}})$, we obtain its
signature vector of length $|SI(n,k)|$. Next, we construct an
$N\times|SI(n,k)|$ matrix whose rows are the $(n,k)-$signature vectors for
each vertex. For a given multiplex network $G$ and its sub-network with plexes
$\alpha_{1},\ldots,\alpha_{k}$, we compute Spearman’s correlation coefficients
between all pairs of columns of the above described matrix and present them in
a $|SI(n,k)|\times|SI(n,k)|$ symmetric matrix which is termed graphlet
correlation matrix of the $k-$plex network
$G(V,E^{\alpha_{1}},\ldots,E^{\alpha_{k}})$. In this way, the network topology
and its local direction patterns, regardless of network size (the number of
vertices) and network volume (the number of edges), are summarized into a
$|SI(n,k)|\times|SI(n,k)|$ matrix for the $k-$plex
$G(V,E^{\alpha_{1}},\ldots,E^{\alpha_{k}})$. Since the Graphlet Correlation
Matrix of a network $G_{1}$ is such that devises a network statistic based on
correlations between the node properties across the multiplex orbits and sub-
orbits, we can examine the network topology of two networks by introducing a
Graphlet Correlation Distance (GCD). Moreover, for two graphs, $G_{1}$ and
$G_{2}$ and their graphlet correlation matrices $GCM_{G_{1}}$ and
$GCM_{G_{2}}$, which are clean of redundancies and encode the information
about the local network topology of a multiplex network that we examine, the
graphlet correlation distance is defined as the Euclidean distance of their
upper triangle values:

$GCD(G_{1},G_{2})=\sqrt{\sum_{i=1}^{d}\sum_{j=i+1}^{d}(GCM_{G_{1}}(i,j)-GCM_{G_{2}}(i,j))^{2}}$
(7)

This metric has been used for a single-plex network by Pržulj et al. [7]. Here
the same distance has been adopted for $k$-plex network, $k\geq 2$,
representing the network topology through the local connectivity.

### Data

Synthetic data – We generate 2000 synthetic multiplex networks, each having
two plexes which are generated from the same algorithm, using four graph
algorithms: Erdos-Renyi (ER), Watts-Strogatz (WS), Barabashi-Albert (BA) and a
modified BA algorithm (PL- powerlaw cluster)[40]. The size of the network
varies between $N=100,200,300,400,500$ and both plexes are generated using the
same algorithm with the same set of parameters. For the PL algorithm the
triangle forming probability is 0.8, while WS has a rewiring probability of
0.01. The other parameters are either the size of the connected component or
the sparsity probability. These are derived as follows:
$p\in\\{0.2,0.35,0.5,0.75,0.8\\},k=\lfloor p*N\rfloor$, where p is the
sparsity probability for ER, while k represents the initial connected
component parameter for BA and PL, and the size of the initial neighbours in
WS. This is fully explained in the SI.

Economic Trade Networks Data – International trade network data from the most
recent available year (2000), provided in [41], are used to construct a
“Multiplex International Trade Network”, in which plexes represent products,
nodes are countries and a link between country $i$ and $j$ in the plex
$\alpha$ exists if at least one of the countries is a significant exporter of
the product $\alpha$ to the other country. The full explanation of how the
network is created is given in the SI. The final network contains of $N=125$
countries (nodes) and $D=957$ products (plexes). From this network we then
focus on a subset of products, or combine products into product categories by
their hierarchical Standard International Trade Classification (SITC) code.

Social Networks Data – Data collected from 75 villages in the region of
Karnatka [42] was used. Each node represents an individual with age ranging
from 18 to 57. The individuals were asked how they interact between each other
in the village, across 12 different aspects of the everyday life: visiting
other’s homes, who they were inviting to visit their’s home, kin, nonrelatives
with whom they socialize, from who they receive medical advice, those from
whom would borrow money and from whom they would lend money, those from which
they would borrow material goods (kerosene, rice, etc.), those to whom they
would lend material goods, giving or getting advice, people with whom they go
to pray (at temple, church or mosque). The data is organized as a multiplex
network with a number of nodes dependent on the village size and $D=12$
plexes.

## Results and Discussions

A typical multiplex network has large number of plexes and, therefore, a full
graphlet analysis is computationally infeasible. Moreover, a network might
contain plexes, which are less significant or contain a small number of links.
Because of this, often one might want to focus on a smaller number of plexes
to be analyzed. Two different strategies have been adopted to address this
problem. First, when the problem in question is such that full graphlet
analysis is needed, we consider for a given multiplex network with $d$ plexes,
the set (or well-defined subset) of all $k-$plex networks
$G(V,E^{\alpha_{1}},\ldots,E^{\alpha_{k}})$ such that
$\alpha_{1}\neq\alpha_{2}\neq\ldots\alpha_{k}$, where $k$ is a small number,
typically $k=2,3,4$. Second, when the problem to be addressed allows sub-orbit
count reduction, we consider both reduced multiplex constructions described in
the previous sections, namely plex-count multiplex and plex-count multiplex
with distinct links inside orbits. In order to compare results for computing
graphlets using full multiplex and two reduced multiplexes, graphlet
correlation matrices (GCMs) are computed for all $k-$plex networks. From these
GCMs we then analyze only those sub-orbit pairs for which significant
correlations (above 0.7) exist in more than 60% of these $k-$plex networks.

Results using all three approaches (full multiplex and two reduced
multiplexes) are comparable and similar to each other for two data sets from
two different domains: economic world trade networks and social networks,
analyzed in more details in this section. For this reason, we only present and
discuss results when the second orbit size reduction is employed. For better
clarity orbits are represented as $O_{i_{x}.j_{y}.k_{z}}$, where $O$ is the
orbit number, $i,j,k$ are the plex counts and $x,y,z$ are the distinction
indices (different if the original plex sets are different). The links are
ordered by the orbit class definition on Table 1 (edges inside symmetric sets
being ordered lexically). As an example the sub-orbit $3_{13.4.25}$ is
rewritten as $3_{2_{x}.1_{x}.2_{y}}$.

#### Economic Trade Networks

The products in the economic trade network are labeled with a hierarchical 4
letter code (SITC Code). The first letters of the code refer to more general
product classes such as products of animal or mineral origin, while the full
code refers to more specific products such as skimmed milk or pork. The full
network contains info for 957 specific products, which are arranged as a
multiplex (957-plex) network.

In order to find more general conclusions we focus on 2 letter product class
pairs. From each of these product classes we then uniformly chose up to 2000
individual products and create 2-plex networks. For example, if the chosen
product classes are paper and furniture products, then from each of them we
can create 2-plex networks using (craft paper, leather furniture) or
(newspaper rolls, wood furniture) etc.

Figure 4:  Histograms containing normalized frequencies of 22 (second
reduction) sub-orbits from 2-plex economic trade networks (a) and social
networks (b). 4950 2-plex social networks and 15390 2-plex economic trade
networks were analyzed. Detailed explanation how the 2-plex networks were
constructed is provided in the main text and in the SI.

Figure 4a shows the histogram (normalized frequencies) of the sub-orbits in
the economic network. Dominant graphlets are wedges, reflecting how economic
networks are built: as trading networks between two countries. Moreover, in
economic trade networks wedge paths and wedge stars are (almost) equally
represented. We now examine the correlations between graphlet sub-orbits. We
randomly chose 100 product class pairs. For each class pair, we select up to
2000 2-plex networks constructed from individual products belonging to the two
classes respectively. For every product 2-plex network, we compute its GCM and
remember the significant correlations ($>0.7$). Subsequently for each class
pair we retain strong correlations appearing in more than 60% of the
individual networks. To generalize about economic networks we than examine
correlations appearing in more than 90 % of the 100 random product class
pairs.

Figure 5:  (a) Full international trade network of 2 plexes: diary and meat
products. The green and orange links represent trade of dairy and meat
respectively, while blue links represent both dairy and meat trade. Small
induced subgraph from United Kingdom’s trade network is further visualized for
better understanding for highlighting different structural patterns: (b)
strong preferences of trading with neighbors of node’s traders, (c) trade hubs
of the network, (d) triangle single product relations.

The significant correlations that exist in a majority of these pairs provide a
more general overview of graphlet correlations in economic trade networks. The
full correlation tables are presented in the SI. Part of the correlations that
emerged and are more interesting and nontrivial to our understanding are shown
on Figure 5 and described thoroughly below. For these networks, only positive
correlations were found, which in the case of wedges and triangles follow an
interesting pattern of behavior which also emerged when we employed the same
method, for 3 and 4-plex networks. This can be reviewed in the figures
provided in the SI, since here we describe in details only the 2-plex
correlations and their meaning and importance.

The correlation $0_{2_{x}}$ and $2_{2_{x}.1_{x}}$ appears in all 100 randomly
chosen economic 2-plexes. This implies that there are trade relations from a
certain node, some of which are strong (trading both products), while others
are weak (trading one product). Another similar strong correlation is between
$2_{2_{x}.2_{x}}$ and $2_{1_{x}.2_{x}}$ (appearing in 98 % of observed
networks). From these correlations we can infer that certain countries(nodes)
act as a type of network hub, ‘specializing’ in two specific products (links),
which they trade interchangeably with different countries. Furthermore from
the frequency histogram on the same networks, it is clear that strong
relations are not common (which makes sense as we choose random two products),
however wedge star relations of the type $2_{2_{x}.1_{x}}$ appear more
frequently.

We also observe triangle trade relations in the economic networks. The sub-
orbits $0_{1_{x}}$ and $3_{1_{x}.1_{x}.1_{x}}$ are strongly correlated in 96 %
of observed networks. This implies that countries form trade deals with a
specific product. This is supported also by the histogram data, which shows
that single product triangles appear more frequently than other triangles in
the networks. Furthermore we observe correlations between $0_{2_{x}}$ and
$3_{2_{x}.1_{x}.1_{y}]}$ or $3_{2_{x}.1_{x}.1_{x}}$, which appear in 97% of
sampled networks. Therefore, when a country trades one or two products with
two different countries, these countries also form trade deals with each other
(with a single product), meaning that one of the countries in the triangle is
a ’stronger’ trader, trading in two products, while the other two countries
mostly trade with a single product. This construction is furthermore supported
with the correlation between $3_{2_{x}.1_{x}.1_{x}}$ and
$3_{2_{x}.1_{x}.1_{y}}$, appearing in 95 % of networks, suggesting that there
are ’trade triangles’ that are composed of a strong trade link $2_{x}$, (which
are rarer) and single product trade links $1_{i}$ composed of any of the two
products. Moreover, the formation of these deals can be supported with the
following correlations: $2_{2_{x}.1_{x}}$ and $3_{2_{x}.1_{x}.1_{x}}$, with 93
% of observed networks; $2_{2_{x}.1_{x}}$ and $3_{2_{x}.1_{x}.1_{y}}$, with 93
% of observed networks; and $2_{1_{x}.1_{x}}$ and $3_{1_{x}.1_{x}.1_{x}}$,
with 92 % of observed networks. This implies that some of the wedge trade
relations are closed into triangles with a single product trade link. The
triangle relations are represented in Fig. 5 b, where we can see that United
Kingdom acts as a hub for trading meat and dairy (thus having a number of
trade relations), and it forms triangle $3_{1_{x}.2_{x}.1_{x}}$ with Austria
and New Zealand, or in Fig. 5 d we have the triangle $3_{1_{x}.1_{x}.1_{x}}$
between Chile, Italy and the trading hub Netherlands.

The concept of economic complexity has been recently introduced [43, 44, 45,
46, 47] with the aim to reflect the amount of knowledge that is embedded in
the productive structure of an economy. Capability-driven economic
competitiveness has been analyzed using three methods: methods for
reflections, fitness-complexity method, and modified fitness-complexity
method. Two simple measures have been introduced both related to degrees: the
first is country degree (in the bipartite network) and is called diversity and
the second is the product degree (in the bipartite network) and is called
ubiquity. Here graphlet analysis provides another view of the economic
competitiveness. For those countries with large diversity (computed as in
[43], for example) we found that sub-orbits $2_{2_{x}.2_{x}}$ and
$2_{1_{x}.2_{x}}$ are correlated as well as sub-orbits $1_{2_{x}.2_{x}}$ and
$1_{1_{x}.2_{x}}$ but not $2_{2_{x}.2_{x}}$ and $2_{1_{x}.1_{x}}$ or other
pairs for wedges. This implies that the countries with large diversity are
structurally (locally) well described with wedges and have both double and
single plex links. On the other hand, the fact that correlations between
degrees and triangles are found among sub-orbits $0_{2_{x}}$ and
$3_{2_{x}.1_{x}.1_{y}}$ or $0_{2_{x}}$ with $3_{2_{x}.1_{x}.1_{x}}$,
$0_{1_{x}}$ and $3_{1_{x}.1_{x}.1_{x}}$, and not between other degree-triangle
pairs, in particular not between $0_{2_{x}}$ and $3_{2_{x}.2_{x}.2_{x}}$
provide evidence that that the countries with small diversity tend to form
correlated triangles.

#### Social Networks

Social network data is organized as 75 (the number of villages) multiplex
networks with a number of nodes dependent on the village size and $D=12$
plexes. Each 12-plex network is further separated into all different
combinations of 2-plexes and the graphlet analysis, similar to one performed
for economic data, is carried out for social data as well. We run our analysis
for all possible relation pairs, 66 in total. For each relation pair we create
2-plex networks for each village, then we extract the strong correlations that
appear in a majority of villages $(>60\%)$. These correlations are assumed to
be representative of the specific relation pair. By finding correlations that
appear in a majority of relation pairs $(>80\%)$ we aim to find more general
correlations that appear in social networks.

Figure 6:  (a) 2-plex social network from one of the villages (village 28):
give advice (orange links) and friends relations (green links). Blue links
represent both, giving advice and friend relations. Focus, again, for better
visualization is given for a smaller induced subgraph. (b) Highlights the
strong tie cliques in the network; (c) structural holes of the network, and
(d) links between the structural holes and the outreach node of a clique.

Figure 4b depicts the histogram (normalized frequencies) of the sub-orbits in
the social network. Dominant graphlets are sub-orbits: degree $0_{2_{x}}$,
wedges with at least one $2_{x}$ link, that is $1_{1_{x}.2_{x}}$,
$1_{2_{x}.1_{x}}$, $1_{2_{x}.2_{x}}$, and $2_{1_{x}.2_{x}}$,
$2_{2_{x}.2_{x}}$, and the triangle $3_{2_{x}.2_{x}.2_{x}}$). A significant
correlation (100% of all social relation combinations, in the majority of
villages) was found between degree sub-orbit $0_{2_{x}}$ and triangle sub-
orbit $3_{2_{x}.2_{x}.2_{x}}$ (shown in Fig. 6b). In the seminal work
Granovetter [48] suggested the strength of dyadic ties to be the tool linking
micro and macro levels of sociological theory. He showed that dyadic ties are
related to larger structures by implementing the following principle: the
stronger the tie between two individuals, the larger the proportion of
individuals to whom they will both be tied. The impact of this principle on
diffusion of influence and information, mobility opportunity, and community
organization is well documented. This principle has been supported by
providing evidence that the triads in which two ties are strong and the third
is absent are unlikely to occur. Following Granovetter here we suggest
multiplexity as a way of indicating a strong tie. Thus, the tie $2_{x}$ (in
which both plexes are presented) is called strong tie, while ties $1_{i}$ are
weak ties. We also call $3_{2_{x}.2_{x}.2_{x}}$ a strong triangle. The
graphlet analysis shows that triangles other than $3_{2_{x}.2_{x}.2_{x}}$ are
unlikely to occur; moreover, the occurrence of the strong triangles is highly
correlated with the occurrence of the strong ties, supporting Granovetter’s
principle. However, in our case, this support is direct. The occurrence of
wedges $1_{2_{x}.2_{x}}$ and $2_{2_{x}.2_{x}}$ are significant (in contrast to
[48] in which such wedges are unlikely to occur). Moreover, we found
significant correlation in 89% of 2-plexes between triads in which two ties
are strong and the third is absent. In another words, we found that the wedge
stars $2_{2_{x}.2_{x}}$ and the wedge paths $1_{2_{x}.2_{x}}$ are strongly
correlated.

Another concept related to weak ties is the concept of structural holes [49],
introduced to explain the origin of differences in social capital. An
individual holds certain positional advantages/disadvantages from how she/he
is embedded in neighborhoods or other social structures. A structural hole
represents the gap between two individuals who have complementary sources of
information. A simple measure of structural holes in a network is the bridge
count. According to Granovetter, no strong tie is a bridge [48]. Several
significant graphlet correlations support the concept of structural holes.
Thus, significant correlations ($>0.7$) are found in 100% of the tested 2-plex
networks between wedge star $2_{1_{x}.2_{x}}$ with $1_{1_{x}.2_{x}}$, also
$1_{1_{x}.1_{x}}$ with $1_{1_{x}.2_{x}}$ in 95% of networks. When looking at
the data we found that these structures appear often around the same person.
Furthermore we observe the following correlations between $0_{1_{x}}$ and
these structures: $2_{1_{x}.2_{x}}$, $1_{1_{x}.2_{x}}$ with 100% of relations
pairs, also $1_{1_{x}.1_{x}}$ with 97%, and $2_{1_{x}.1_{x}}$ with 94% of
relation pairs. If we examine the connections that these people have within
the network, we can conclude that they behave as structural holes, meaning
they connect different social cliques (households). This can be viewed on
Figure 6c which depicts several examples of individuals that are mediators
between two or more households as part of different cliques. A correlation
that appears in 83% of networks is between wedge paths $1_{1_{x}.1_{y}}$ with
$1_{1_{x}.2_{x}}$. This correlation again signifies that the network has
individuals that serve as a connection between two cliques. Furthermore, these
graphlet structures can also detect individuals inside the clique that are
linked to outside, having broader communication reach. This can be observed on
Figure 6d where only a few people inside a clique are connected to an outside
person. For the full list of correlations we refer the reader to the SI.

We remark that these conclusions are based on the analysis of a single data
set consisting of 75 social networks. Therefore, more thoughtful social
network analysis for reaching more general conclusion is needed, which is,
however, beyond the scope of this paper and will be provided in a future
study.

Figure 7:  Synthetic networks: 3D Visualization using graphlet correlations.
Purple points: Erdos-Renyi (ER) network, gray points: Watts-Strogatz (WS)
network, blue points: Barabashi-Albert (BA) network, and orange points:
modified BA (PL- power-law cluster) network. (a), (b), and (c) flattened
networks of 2, 3 and 4 plexes; (d), (e), and (f) 2-plex, 3-plex and 4-plex
networks, respectively.  Figure 8:  3D visualization using graphlet
correlations. Green points are the economic 2-plex networks while social
2-plexes are shown with cyan points. (a) For better visualization, 2% of
2-plex flattened networks are shown; (b) full 2-plex networks; (c) social and
economic 2-plex networks positioned together with the synthetic 2-plex
networks.

## Conclusions

Graphlets are a powerful tool for analyzing local network structure. Multiplex
networks, multilayer networks, and networks with node and/or link
(categorical) attributes are pervasive and graphlet analysis developed here
can further enhance our understanding of complex networks. Graphlets provide
discriminatory property for different type networks. Even a simple graphlet
histogram plot of economic and social networks, see Fig. 4, shows the
differences between these structures and provides evidences on how the
networks are built. Wedges occur more often in economic networks rather than
social networks, indicating the tendency of a country to produce/trade of a
product in local structure of triads which are not closed (that is, wedges not
triangles). Wedges (open triads) also appear in the social networks, however
the dominant graphlets in social networks are triangles (closed triads). If
the multiplex is the indicator of the strong tie, the graphlet analysis
provides another evidence for the concepts of strong/weak ties and structural
holes. In contrast to the work of Granovetter[48], however, in our work
related to a single data set consisting of 75 social networks, wedges with
only strong ties are not only present but they are strongly correlated.

Graphlets can also provide clustering. Graphlet correlation matrix of a given
network is represented as point in some multidimensional space, which is then
visualized in 3D space, by using multidimensional scaling. This has been
demonstrated with synthetic networks in Fig. 7. Several conclusions can be
drawn from the figure: (1) different graphs are clearly separated, (2) same
graphs with different sparsity are also distinguishable, and (3) aggregating a
multiplex network in a single-plex network provides less information on the
network structure. Graphlet correlation matrices for economic and social
networks are visualized in Figure 8. We again notice a good separation of the
economic and social networks. This is present both in the flattened and the
2-plex representation. However, due to the additional plex information, 2-plex
graph clouds have better resolution. The final part of Figure 8 shows the
synthetic and real networks in the same 3d space, and again there exists a
strong separation among the different graph types. The full description on how
graphlets can be used for clustering will be provided in a separate
manuscript.

## References

  * [1] Milo, R. _et al._ Network motifs: simple building blocks of complex networks.  _Science_ 298, 824–827 (2002).
  * [2] Sporns, O. & Kötter, R.  Motifs in brain networks.  _PLoS Biol_ 2, e369 (2004).
  * [3] Pržulj, N., Corneil, D. G. & Jurisica, I.  Modeling interactome: scale-free or geometric?  _Bioinformatics_ 20, 3508–3515 (2004).
  * [4] Holland, P. W. & Leinhardt, S.  Local structure in social networks.  _Sociological methodology_ 7, 1–45 (1976).
  * [5] Faust, K.  A puzzle concerning triads in social networks: Graph constraints and the triad census.  _Social Networks_ 32, 221–233 (2010).
  * [6] Pržulj, N.  Biological network comparison using graphlet degree distribution.  _Bioinformatics_ 23, e177–e183 (2007).
  * [7] Yaveroğlu, Ö. N. _et al._ Revealing the hidden language of complex networks.  _Scientific reports_ 4, 4547 (2014).
  * [8] Sarajlić, A., Malod-Dognin, N., Yaveroğlu, Ö. N. & Pržulj, N.  Graphlet-based characterization of directed networks.  _Scientific reports_ 6 (2016).
  * [9] Trpevski, I., Dimitrova, T., Boshkovski, T., Stikov, N. & Kocarev, L.  Graphlet characteristics in directed networks.  _Scientific Reports_ 6 (2016).
  * [10] Boorman, S. A. & White, H. C.  Social structure from multiple networks. ii. role structures.  _American journal of sociology_ 1384–1446 (1976).
  * [11] White, H. C., Boorman, S. A. & Breiger, R. L.  Social structure from multiple networks. i. blockmodels of roles and positions.  _American journal of sociology_ 730–780 (1976).
  * [12] Gould, R. V.  Multiple networks and mobilization in the paris commune, 1871.  _American Sociological Review_ 716–729 (1991).
  * [13] Padgett, J. F. & Ansell, C. K.  Robust action and the rise of the medici, 1400-1434.  _American journal of sociology_ 1259–1319 (1993).
  * [14] Granovetter, M.  Economic action and social structure: The problem of embeddedness.  _American journal of sociology_ 91, 481–510 (1985).
  * [15] Baldassarri, D. & Diani, M.  The integrative power of civic networks 1.  _American Journal of Sociology_ 113, 735–780 (2007).
  * [16] Grossmann, M. & Dominguez, C. B.  Party coalitions and interest group networks.  _American Politics Research_ 37, 767–800 (2009).
  * [17] Maggioni, M. A., Breschi, S. & Panzarasa, P.  Multiplexity, growth mechanisms and structural variety in scientific collaboration networks.  _Industry and Innovation_ 20, 185–194 (2013).
  * [18] Rank, O. N., Robins, G. L. & Pattison, P. E.  Structural logic of intraorganizational networks.  _Organization Science_ 21, 745–764 (2010).
  * [19] Ferriani, S., Fonti, F. & Corrado, R.  The social and economic bases of network multiplexity: Exploring the emergence of multiplex ties.  _Strategic Organization_ 11, 7–34 (2013).
  * [20] Shipilov, A. V. & Li, S. X.  The missing link: The effect of customers on the formation of relationships among producers in the multiplex triads.  _Organization Science_ 23, 472–491 (2012).
  * [21] De Domenico, M. _et al._ Mathematical formulation of multilayer networks.  _Physical Review X_ 3, 041022 (2013).
  * [22] Kivelä, M. _et al._ Multilayer networks.  _Journal of complex networks_ 2, 203–271 (2014).
  * [23] Solá, L. _et al._ Eigenvector centrality of nodes in multiplex networks.  _Chaos: An Interdisciplinary Journal of Nonlinear Science_ 23, 033131 (2013).
  * [24] Halu, A., Mondragón, R. J., Panzarasa, P. & Bianconi, G.  Multiplex pagerank.  _PloS one_ 8, e78293 (2013).
  * [25] Battiston, F., Nicosia, V. & Latora, V.  Structural measures for multiplex networks.  _Physical Review E_ 89, 032804 (2014).
  * [26] De Domenico, M., Solé-Ribalta, A., Omodei, E., Gómez, S. & Arenas, A.  Ranking in interconnected multilayer networks reveals versatile nodes.  _Nature communications_ 6 (2015).
  * [27] Cozzo, E. _et al._ Structure of triadic relations in multiplex networks.  _New Journal of Physics_ 17, 073029 (2015).
  * [28] Bianconi, G.  Statistical mechanics of multiplex networks: Entropy and overlap.  _Physical Review E_ 87, 062806 (2013).
  * [29] Menichetti, G., Remondini, D., Panzarasa, P., Mondragón, R. J. & Bianconi, G.  Weighted multiplex networks.  _PloS one_ 9, e97857 (2014).
  * [30] Cardillo, A. _et al._ Emergence of network features from multiplexity.  _Scientific Reports_ 3, 1344 (2013).
  * [31] Braun, U. _et al._ Dynamic reconfiguration of frontal brain networks during executive cognition in humans.  _Proceedings of the National Academy of Sciences_ 112, 11678–11683 (2015).
  * [32] Filiposka, S., Gajduk, A., Dimitrova, T. & Kocarev, L.  Bridging online and offline social networks: Multiplex analysis.  _Physica A: Statistical Mechanics and its Applications_ 471, 825–836 (2017).
  * [33] Morris, R. G. & Barthelemy, M.  Transport on coupled spatial networks.  _Physical review letters_ 109, 128703 (2012).
  * [34] Solé-Ribalta, A., Gómez, S. & Arenas, A.  Congestion induced by the structure of multiplex networks.  _Physical review letters_ 116, 108701 (2016).
  * [35] Wang, Z., Andrews, M. A., Wu, Z.-X., Wang, L. & Bauch, C. T.  Coupled disease–behavior dynamics on complex networks: A review.  _Physics of life reviews_ 15, 1–29 (2015).
  * [36] Funk, S. _et al._ Nine challenges in incorporating the dynamics of behaviour in infectious diseases models.  _Epidemics_ 10, 21–25 (2015).
  * [37] Granell, C., Gómez, S. & Arenas, A.  Dynamical interplay between awareness and epidemic spreading in multiplex networks.  _Physical review letters_ 111, 128701 (2013).
  * [38] Sanz, J., Xia, C.-Y., Meloni, S. & Moreno, Y.  Dynamics of interacting diseases.  _Physical Review X_ 4, 041005 (2014).
  * [39] Lima, A., De Domenico, M., Pejovic, V. & Musolesi, M.  Disease containment strategies based on mobility and information dissemination.  _Scientific reports_ 5, 10650 (2015).
  * [40] Holme, P. & Kim, B. J.  Growing scale-free networks with tunable clustering.  _Physical review E_ 65, 026107 (2002).
  * [41] Feenstra, R. & Lipsey, R.  Nber-united nations trade data 1962-2000 (2000).
  * [42] Banerjee, A., Chandrasekhar, A. G., Duflo, E. & Jackson, M. O.  The diffusion of microfinance.  _Science_ 341, 1236498 (2013).
  * [43] Hidalgo, C. A., Klinger, B., Barabási, A.-L. & Hausmann, R.  The product space conditions the development of nations.  _Science_ 317, 482–487 (2007).
  * [44] Hidalgo, C. A. & Hausmann, R.  The building blocks of economic complexity.  _proceedings of the national academy of sciences_ 106, 10570–10575 (2009).
  * [45] Tacchella, A., Cristelli, M., Caldarelli, G., Gabrielli, A. & Pietronero, L.  A new metrics for countries’ fitness and products’ complexity.  _Scientific reports_ 2, 723 (2012).
  * [46] Cristelli, M., Gabrielli, A., Tacchella, A., Caldarelli, G. & Pietronero, L.  Measuring the intangibles: A metrics for the economic complexity of countries and products.  _PloS one_ 8, e70726 (2013).
  * [47] Stojkoski, V., Utkovski, Z. & Kocarev, L.  The impact of services on economic complexity: Service sophistication as route for economic growth.  _PloS one_ 11, e0161633 (2016).
  * [48] Granovetter, M. S.  The strength of weak ties.  _American journal of sociology_ 1360–1380 (1973).
  * [49] Holes, S.  The social structure of competition (1992).

## Author contributions statement

All authors conceived the experiments. T.D. and K.P. conducted the
experiments. All authors analyzed the results, wrote and reviewed the
manuscript.

## Additional information

Competing financial interests The authors declare no conflicts of interest.


# Effect of Holstein phonons on the optical conductivity of gapped graphene

Kh. Jahanbani  Institute for Advanced Studies in Basic Sciences (IASBS),
Zanjan, 45195-1159, Iran  School of Physics, Institute for Fundamental
Sciences, (IPM) 19395-5531 Tehran, Iran  Reza Asgari 111Corresponding author:
Tel: +98 21 22280692; fax: +98 21 22280415.
E-mail address<EMAIL_ADDRESS> School of Physics, Institute for
Fundamental Sciences, (IPM) 19395-5531 Tehran, Iran

###### Abstract

We study the optical conductivity of a doped graphene when a sublattice
symmetry breaking is occurred in the presence of the electron-phonon
interaction. Our study is based on the Kubo formula that is established upon
the retarded self-energy. We report new features of both the real and
imaginary parts of the quasiparticle self-energy in the presence of a gap
opening. We find an analytical expression for the renormalized Fermi velocity
of massive Dirac Fermions over broad ranges of electron densities, gap values
and the electron-phonon coupling constants. Finally we conclude that the
inclusion of the renormalized Fermi energy and the band gap effects are indeed
crucial to get reasonable feature for the optical conductivity.

###### pacs:

78.67.-n, 71.10.Ay, 73.25.+i, 72.80.-r

##  1 Introduction

There is a considerable interest in understanding the effects on properties of
particle due to the interactions with environment, for instance the coupling
of electrons to lattice vibrations or electron-phonon coupling. The electron-
phonon coupling plays an essential role in the theory of high temperature
superconductivity and they exist in other material such as nanotubes, C60
molecules and other fullerenes rmp . Also it is important to consider the
electron-phonon coupling in transport properties.

Graphene, a single layer of carbon atoms, Geim  is disputable the first true
two-dimensional lattices. Graphene is thermodynamically stable and there is
indeed ripple structures on graphene sheets. Lattice displacements due to the
ripple structures are symmetric with respect to their close carbon atoms and
couple to the carrier densities. The electrons moving through the sheet are
coupled to the out-of-plane phonons and therefore the electron-phonon coupling
plays an important role in the transport properties akturk ; basko ; park .
The coupling of electrons to out-of-plane optical phonons can be modeled by a
Holstein type coupling Holstein . In this model the coupling of electrons to
dispersionless optical phonons is essentially local. The electron-phonon
coupling has been carefully examined and has been shown to give rise to Kohn
anomalies in the phonon dispersion at edge points in the Brillouin zone where
the phonons can be studied by Raman spectroscopy piscanec1 ; piscanec2 ;
pisana . An alternative strategy for the electron-phonon coupling measurement
is based on the analysis of the $G$-peak linewidths and its broadening.

The optical conductivity is one of the most useful tools to investigate the
basic properties of materials. Both the excitation spectrum of materials such
gaps, phonons and interband transitions and the scattering mechanisms leave
their distinct traces in transport. It was shown that the infrared
conductivity of graphene is basically independent of the frequency peres2 ;
peresijmp ; peresprb78 ; gusynin  and experimentally confirmed this manner li
; nair . The effect of electron-phonon interaction in gapless graphene has
been discussed by several authors Stauber ; calandra ; tse ; peres ; stauperes
directed towards understanding this effect on the optical conductivity.

The energy spectrum of the Dirac electrons in a graphene layer that
epitaxially grown on a SiC substrate has been measured by Zhou et al. Zhou
and they observed an energy gap of about 200 meV opened up in the electronic
spectrum. They attributed the opening up of the gap is due to the breaking of
the $A$ and $B$ sublattices symmetry novoselov . The optical response of a
gapped graphene is of important for an understanding of optoelectronic
devices. Moreover, the optical spectroscopy can be used for measurements of
the magnitude of the energy gap.

In this paper we consider the sublattice symmetry breaking mechanism for a gap
opening in a pristine doped graphene sheet and study the impact of the
electron-phonon coupling on the electronic conductivity of the electron-doped
gapped graphene using Kubo formula at zero-temperature. We show that the
renormalized velocity is suppressed due to the electron-phonon interaction.
There is a shift in the chemical potential and we show that the interacting
chemical potential is less than the noninteracting one due to the electron-
phonon coupling. The optical conductivity is affected by Pauli blocking below
twice value of the renormalized interacting chemical potential and gap values.

##  2 Model Hamiltonian and theory

We consider the simplest form of Hamiltonian that describes the interaction of
electron with an optical phonon mode, called the Holstein model. The honeycomb
lattice can be consider in terms of two triangular sublattices $A$ and $B$. We
consider electrons in $\pi$-orbital of carbon atoms by using the tight-binding
Hamiltonian in addition to the effect of the electron-phonon coupling due to
localized Holstein phonons and a gap opening procedure due to sublattice
symmetry breaking alireza . The total Hamiltonian in momentum space can be
expressed as

$\displaystyle H$ $\displaystyle=$
$\displaystyle-t\sum_{k,\sigma}[\phi(k)a_{{\bf k},\sigma}^{\dagger}b_{{\bf
k},\sigma}+h.c]$ (1) $\displaystyle+$ $\displaystyle
D\sum_{p,k,\sigma}\chi_{0}[a^{\dagger}_{{\bf p},\sigma}a_{{\bf
p+k},\sigma}+b^{\dagger}_{{\bf p},\sigma}b_{{\bf p+k},\sigma}](c_{\bf
k}+c_{\bf-k}^{\dagger})$ $\displaystyle+$
$\displaystyle\sum_{k}\omega_{0}c_{\bf k}^{\dagger}c_{\bf
k}+\Delta\sum_{k}[a_{{\bf k},\sigma}^{\dagger}a_{{\bf k},_{\sigma}}-b_{{\bf
k},\sigma}^{\dagger}b_{{\bf k},\sigma}]$ $\displaystyle-$
$\displaystyle\mu_{0}\sum_{k}[a_{{\bf k},\sigma}^{\dagger}a_{{\bf
k},\sigma}+b_{{\bf k},\sigma}^{\dagger}b_{{\bf k},\sigma}]$

where $a_{{\bf k},\sigma}$ or $b_{{\bf k},\sigma}$ is the fermion annihilation
operator in $k-$space on sublattice $A$ or $B$, respectively and $t$ is the
nearest neighbor hopping parameter castro . The band gap, $2\Delta$ has a
nonzero value as a result of breaks the symmetry between sublattices, A and B.
We consider that the noninteracting chemical potential, $\mu_{0}$ be larger
than the gap value representing the electron-doped system. The electron-phonon
coupling is determined by $D$ and furthermore $\omega_{n}$ denotes the
fermionic Matsubara frequency. Moreover, $c_{\bf k}$ is the annihilation
phonon operator. $\omega_{0}$ is the frequency of the out of plane vibrations
of the optical phonon and $\chi_{0}=\sqrt{\frac{\hbar}{2MN\omega_{0}}}$ with M
is ion’s mass and N denotes the number of unit cells.
$\phi(k)=\sum_{\delta}\exp^{-i{\bf\delta}\cdot{\bf k}}$ with $\delta$ being
the vectors connecting the three nearest neighbors on the honeycomb lattice
Mahan . $\phi(k)$ reduces to $\hbar v_{\rm F}k/t$ in the Dirac cone
approximation castro .

The matrix element of noninteracting Green’s function with the gap of the
electronic spectrum is determined by following expression

$\displaystyle G_{\alpha\beta}^{0}(k,i\omega_{n})=\frac{1}{2}\sum_{\lambda=\pm
1}\left(\delta_{\alpha,\beta}+\frac{\lambda\Upsilon_{\alpha,\beta}}{\xi_{{}_{\Delta}}(k)}\right)\frac{1}{i\hbar\omega_{n}+\mu_{0}-\lambda\xi_{{}_{\Delta}(k)}},$

in which $\alpha,\beta=A,B$ and we have defined parameters
$\Upsilon_{AA}=-\Upsilon_{BB}=\Delta$,
$\Upsilon_{AB}=\Upsilon^{*}_{BA}=-t\phi(k)$. The quasiparticle excitation
energy is $\xi_{{}_{\Delta}}(k)=\sqrt{t^{2}|\phi(k)|^{2}+\Delta^{2}}$. Note
that at zero-temperature $\mu_{0}=\xi_{{}_{\Delta}}(k_{\rm F})$ with $k_{\rm
F}$ is the Fermi momentum of charge carriers.

An exact evaluation of the self-energy is only possible in some special cases.
The matrix elements of the self-energy calculated to the lowest order in the
electron-phonon interaction and is defined as

$\Sigma_{\alpha,\beta}(i\omega_{n},p)=-k_{B}T\sum_{k,\nu}D^{2}\chi_{0}^{2}D^{(0)}(k,i\nu)G^{(0)}_{\alpha,\beta}(p-q,i\omega_{n}-i\nu)$
(2)

where $G^{(0)}_{\alpha,\beta}$ and $D^{(0)}$ are the zero-order electron and
phonon Green’s functions, respectively Mahan ; Jonson . In Holstein phonons,
$\chi_{0}$ and $D^{(0)}(k,i\nu)$ are momentum independent and thus the phonon
propagator is simplified by

$D^{(0)}(k,i\nu)=-2\omega_{0}/(\nu^{2}+\omega_{0}^{2}).$ (3)

We restrict our calculations to the lowest order self-energy that is
sufficient if Migdal’s theorem, states that vertex corrections in the
electron-phonon interaction can be neglected if the typical phonon frequencies
are sufficiently smaller than the electronic energy scale, is valid.
Therefore, we can neglect the vertex corrections since the self-energy is
$k$-independent. Using the contour integration, we can perform the summation
over the bosonic frequency in the expression of the self-energy and finally
the self-energy yields as

$\displaystyle\Sigma_{AA}(i\omega_{n})$ $\displaystyle=$
$\displaystyle\frac{D^{2}\chi_{0}^{2}}{2}\sum_{k,\lambda=\pm
1}\left(1+\frac{\lambda\Delta}{\xi_{{}_{\Delta}(k)}}\right)\times$ (4)
$\displaystyle\left\\{\frac{N_{0}+n_{F}(\lambda\xi_{{}_{\Delta}(k)}-\mu_{0})}{i\hbar\omega_{n}+\hbar\omega_{0}-\lambda\xi_{{}_{\Delta}(k)}+\mu_{0}}\\!+\\!\frac{N_{0}+1-n_{F}(\lambda\xi_{{}_{\Delta}(k)}-\mu_{0})}{i\hbar\omega_{n}-\hbar\omega_{0}-\lambda\xi_{{}_{\Delta}(k)}+\mu_{0}}\right\\},$

where $N_{0}=1/(e^{\hbar\omega_{0}/k_{B}T}-1)$ and $n_{F}(x)$ denotes the
Fermi-Dirac distribution function. To calculate $\Sigma_{BB}(i\omega_{n})$,
the gap value $\Delta$ might be replaced by $-\Delta$ in Eq. 4. It should be
noted that $\Sigma_{AB}(i\omega_{n})=\Sigma_{BA}(i\omega_{n})=0$ in the Dirac
cone approximation. The explicit expression of the self-energy will be
computed in the following.

###  2.1 Finite doping with a gap opening

We consider the low excited electron energy where the noninteracting electron
spectrum energy is given by $\sqrt{(\hbar v_{\rm F}k)^{2}+\Delta^{2}}$ alireza
. To evaluate the zero-temperature retarded self-energy evaluated at the Fermi
surface, we integrate Eq. 4 over $k$ and then decompose the results into
$\Sigma^{j}(\omega)=\Sigma^{j}_{0}(\omega)+\Delta\Sigma^{j}(\omega)$ where

$\displaystyle\Re e\Sigma^{j}_{0}(\omega)=\frac{\hbar
A_{c}}{2\pi}(g\omega_{0})^{2}\\{-\frac{\omega_{j}}{v_{F}^{2}}\ln|\frac{\Delta^{2}+\hbar^{2}v_{F}^{2}k_{c}^{2}}{(\hbar\omega+\mu_{0})^{2}-(\hbar\omega_{0}+\Delta)^{2}}|+\frac{\omega_{0}}{v_{F}^{2}}\ln|\frac{\omega_{+}+\omega_{0}}{\omega_{-}-\omega_{0}}|\\}$
(5) $\displaystyle\Im m\Sigma^{j}_{0}(\omega)$ $\displaystyle=$
$\displaystyle-\pi\frac{\hbar
A_{c}}{2\pi}(g\omega_{0})^{2}\\{-\frac{\omega_{j}+\omega_{0}}{v_{F}^{2}}\Theta(-\omega_{+}-\omega_{0})\Theta(\hbar\omega+\mu_{0}+\hbar\omega_{0}+\sqrt{\hbar^{2}v_{F}^{2}k_{c}^{2}+\Delta^{2}})$
$\displaystyle+$
$\displaystyle\frac{\omega_{j}-\omega_{0}}{v_{F}^{2}}\Theta(\omega_{-}-\omega_{0})\Theta(-\hbar\omega-\mu_{0}+\hbar\omega_{0}+\sqrt{\hbar^{2}v_{F}^{2}k_{c}^{2}+\Delta^{2}})\\}$

here $\omega_{j}=\omega+(\mu_{0}+j\Delta)/\hbar$ with $j=+1(-1)$ refers to
sublattice $A$($B$). $k_{c}$ is the ultraviolet cut-off momentum Stauber  and
finally the coupling constant $g=\sqrt{N}D\chi_{0}/\omega_{0}$ being the order
of unity. The area of the unit cell is $A_{c}=a^{2}3\sqrt{3}/2$ with
$a=1.42$Å. The extra terms take the following form as

$\displaystyle\Re e\Delta\Sigma^{j}(\omega)=\frac{\hbar
A_{c}}{2\pi}(g\omega_{0})^{2}\\{-\frac{\omega_{j}}{v_{F}^{2}}\ln|\frac{(\omega+\omega_{0})(\omega_{-}-\omega_{0})}{(\omega-\omega_{0})(\omega_{-}+\omega_{0})}|-\frac{\omega_{0}}{v_{F}^{2}}\ln|\frac{\omega^{2}-\omega_{0}^{2}}{(\omega_{-})^{2}-\omega_{0}^{2}}|\\}$
(6) $\displaystyle\Im m\Delta\Sigma^{j}(\omega)=-\pi\frac{\hbar
A_{c}}{2\pi}(g\omega_{0})^{2}\\{$
$\displaystyle\\!\\!\\!\\!\\!\\!\frac{\omega_{j}+\omega_{0}}{v_{F}^{2}}~{}\Theta(\omega_{-}+\omega_{0})\Theta(-\omega-\omega_{0})$
(7) $\displaystyle-$
$\displaystyle\frac{\omega_{j}-\omega_{0}}{v_{F}^{2}}\Theta(-\omega+\omega_{0})~{}\Theta(\omega_{-}-\omega_{0})\\}$

If $\Delta=0$, the self-energy reduces to massless Dirac graphene which
addressed in Ref Stauber . Therefore, we have generalized the retarded self-
energy expression to gapped graphene. Once the retarded self-energy is
obtained, the quasiparticle properties of system due to the interaction of the
electron-phonon can be calculated. The renormalized electronic spectrum is
given by the Dyson equation as $E_{\bf k}=\xi_{\Delta}(k)+\Re e\Sigma(E_{\bf
k})$. Notice that according to the Dyson equation, we might distinguish the
noninteracting chemical potential from the chemical potential of the
interacting system due to the fact that $\Re e\Sigma(\omega)$ is not vanished
for doped graphene when $\omega$ tends to zero. We thus have

$\mu=\mu_{0}+\Re e\Sigma(\omega)|_{\omega=0}~{}.$ (8)

The renormalized velocity, on the other hand, is given by

$\frac{v^{\star}}{v_{\rm F}}=\frac{\hbar v_{\rm F}k/\xi_{\Delta}(k)+(v_{\rm
F}\hbar)^{-1}\partial_{k}\Re
e\Sigma(k,\omega)}{1-\hbar^{-1}\partial_{\omega}\Re
e\Sigma(k,\omega)}|_{\omega=0,k=k_{F}}$ (9)

within the Dyson scheme Mahan . The self-energy is independent of the
momentum, accordingly its $k$-derivative is zero. Consequently, the
renormalized velocity is obtained analytically

$\displaystyle\frac{v_{\rm F}}{v^{*}(1+(\Delta/\hbar v_{\rm F}k_{\rm
F})^{1/2})}$ $\displaystyle=$ $\displaystyle
1+\left(\frac{g{\omega}_{0}}{v_{\rm
F}{k}_{c}}\right)^{2}\\{\ln|\frac{({\Delta}^{2}+({\hbar v_{\rm
F}k_{c}})^{2})}{(({\mu_{0}}+{\hbar\omega}_{0})^{2}-{\Delta}^{2})}|$ (10)
$\displaystyle-$
$\displaystyle({\mu_{0}}+j{\Delta}+\hbar{\omega}_{0})\frac{2({\mu_{0}}+\hbar{\omega}_{0})}{({\mu_{0}}+\hbar{\omega}_{0})^{2}-{\Delta}^{2}}+2\frac{{\mu_{0}}+j{\Delta}}{\hbar{\omega}_{0}}\\}.$

###  2.2 Optical Conductivity

The optical conductivity can be calculated from the Kubo formalism. To this
end, we need to obtain the current operator which is a composition of the
paramagnetic and diamagnetic terms, i.e.
$j_{\alpha}=j^{P}_{\alpha}+j^{D}_{\alpha\beta}{\emph{A}}_{\beta}$. We do need
to modify the hopping parameter in the presence of an electromagnetic field
Stauber  and then expand it up to the second order in the vector potential
$\overrightarrow{\emph{A}}(t)$. The current operator expressions do not change
in the presence of the gap value and therefore by assuming that the electric
field is in the direction of $x$-axis, we have

$j_{x}^{P}=-i\xi\sum_{\sigma,k}[(\phi(k)-3)a_{\sigma}^{\dagger}(k)b_{\sigma}(k)-(\phi^{\star}(k)-3)a_{\sigma}(k)b_{\sigma}^{\dagger}(k)]$
(11)

where $\xi=tea/\hbar$ and then the Kubo formula for conductivity is given by

$\sigma_{xx}(\omega)=\frac{<j_{x}^{D}>}{iA_{s}(\omega+i\eta)}+\frac{\Lambda_{xx}(\omega+i\eta)}{i\hbar
A_{s}(\omega+i\eta)}$ (12)

where $A_{s}$ is the area of sample and
$\Lambda_{xx}(i\omega_{n})=\int_{0}^{\hbar/k_{B}T}d\tau
e^{i\omega_{n}\tau}<T_{\tau}j_{x}^{P}(\tau)j_{x}^{P}(0)>$ Mahan . We have
ignored vertex corrections in the Kubo formula since we worked in nearly
highly electron doped graphene for which the Dirac cone approximation is
applicable. It was shown that the vertex corrections is essential for the low
density carriers of the DC conductivity of graphene. cappelluti  After a
lengthy but straightforward algebra, we find

$\displaystyle\Im m\Lambda_{xx}(\omega)\\!$ $\displaystyle=$
$\displaystyle\\!\\!\xi^{2}\frac{A_{s}}{8\pi}\\!\int_{0}^{k_{c}}\\!kdk\int_{-\infty}^{\infty}\frac{d\epsilon}{2\pi}(n_{F}(\epsilon+\omega)-n_{F}(\epsilon))$
(13) $\displaystyle\times$
$\displaystyle\left\\{(2t^{2}|\phi(k)|^{4})A_{AB}(k,\epsilon)A_{AB}(k,\epsilon+\omega)\right.$
$\displaystyle+$
$\displaystyle(9-|\phi(k)|^{2})~{}[A_{AA}(k,\epsilon)A_{BB}(k,\epsilon+\omega)$
$\displaystyle+$ $\displaystyle
A_{BB}(k,\epsilon)A_{AA}(k,\epsilon+\omega)]\\}$

where the spectral functions are the imaginary part of Green’s function which
take the following forms:

$\displaystyle A_{\alpha,\beta}=-2\Im
m\\{\frac{\Phi_{\alpha,\beta}}{{(\Omega_{+}-\Sigma_{BB}(i\omega_{n}))(\Omega_{-}-\Sigma_{AA}(i\omega_{n}))-t^{2}|\phi(k)|^{2}}}\\}.$

Here

$\displaystyle\Omega_{\pm}=i\hbar\omega_{n}+\mu\pm\Delta,~{}~{}~{}~{}~{}~{}~{}~{}~{}\Phi_{AA(BB)}=\Omega_{+(-)}-\Sigma_{BB(AA)}(i\omega_{n})$

and $\Phi_{AB}=\Phi_{BA}=1$. The integral over $k$ in Eq. 13 can be performed
analytically and accordingly one dimensional integral will be needed to be
calculated numerically. Note that the interacting chemical potential is used
instead of the noninteracting one because of the nonzero value of
$\Sigma_{j}(0)$ . It should be noted that by setting $\Delta=0$, the optical
conductivity results are different with the results given in Ref. Stauber  due
to the fact that we have implemented the interacting Fermi energy in the
formalism.

##  3 Numerical Results

We have considered the system with the phonon energy being
$\hbar\omega_{0}=0.2~{}eV$ peres . Although the order of coupling constant is
unity, we consider a larger value to seek its effect better. We have found
that the value of the quasiparticle properties for sublattices $A$ and $B$ are
different at most about $0.8\%$ due to the gap opening. We will then present
only the results of the sublattice $A$.

In Fig. 1, we have shown the results of the real and imaginary parts of the
retarded self-energy for the electron-doped system, ($\mu>\Delta$) at
$n=5\times 10^{12}$cm-2. $\Im m\Sigma(\omega)$ vanishes in
$|\omega|<\omega_{0}$ at which point it jumps up to a finite value because
only then can a quasiparticle decay by boson emission. It drops towards the
zero for $\omega<-\omega_{0}$ and then increases linearly showing a marginal
type physics which happens in the Coulomb electron interactions in undoped
graphene polini . Notice that $\Im m\Sigma$ is not symmetric with respect to
change of the sign of frequency. In addition, $\Im m\Sigma(\omega)$ vanishes
when $|\hbar\omega+\hbar\omega_{0}+\mu_{0}|<\Delta$ due to the effect of the
gap opening and $\Im m\Sigma$ tends to zero at $-\hbar\omega_{0}-\mu_{0}$ for
$\Delta=0$. These behaviors can be determined explicitly from expressions
given by Eqs. 2.1 and 7.

In Fig. 1b we can see logarithmic type singularities dogan  at
$\omega=\pm\omega_{0}$ and $\omega=-\omega_{0}-(\mu_{0}-\Delta)/\hbar$ for the
results of $\Re e\Sigma$. The extra singular behavior is due to the gap
effect. It should be noted that the singularity at $\omega=\pm\omega_{0}$
would be washed out if a momentum dependence of phonon spectra is used. In
addition, there is a cancelation of the logarithmic singularity at
$\omega=-\omega_{0}-(\mu_{0}+\Delta)/\hbar$. The logarithmic singularity can
be determined to the argument of the logarithm in Eqs. 5 and 6. We have
obtained an expression for the interacting density of states too through the
spectral function. The singularities manner lead to kink structures in the
interacting electronic density of states. In the results, there are three kink
structures in the interacting density of states where one of them is
associated to the gap. The kink structures would affect to physical quantities
and transport properties through the interacting electronic density of states.

The renormalized velocity as function of the densities, gap values and the
coupling constants are shown in Fig. 2. The renormalized velocity is
suppressed due to the electron-phonon interaction and the gap values too. We
have found a nonmonotonic behavior of $v^{*}$ with respect to the electron
density when the gap value increases and results are shown in Fig. 2b. At
small gap values, $v^{*}$ decreases with increasing density however it changes
behavior at large gap values and behaves like conventional two-dimensional
electron systems. Therefore, we expect that the electron-phonon interaction
renormalized the electronic quantities at the Fermi surface by a factor
$v^{*}/v_{\rm F}$.

The optical conductivity scaled by $\sigma_{0}=e^{2}/4\hbar$ as a function of
energy for different values of (a) the coupling constants and (b) the gap
values are shown in Fig. 3. First of all, $\sigma$ tends to a minimum value at
$\omega_{0}$. Moreover it basically increases around $\omega>\omega_{0}$ due
to the contribution of the Holstein phonon sideband. In the case of
noninteracting electron-phonon system, $\sigma$ has a sharp structure, step
function manner, at $2\mu_{0}$ due to interband transitions and the
conductivity increases by a factor of two, $\sigma=2\sigma_{0}$ at
$\omega=2\Delta$ and finally at higher frequencies decreases and approaches to
$\sigma_{0}$ gusynin . By switching interaction on, the chemical potential
becomes weaker and consequently the position of the sharp structure changes to
$2\mu$ which is smaller than $2\mu_{0}$. This behavior is clearly shown in the
Fig. 3 which did not consider in results discussed in Ref. Stauber . At $g=0$,
the conductivity is larger than $\sigma_{0}$ about $2\mu$ and then tends to
$\sigma_{0}$ in gapped graphene. However, $\sigma$ always remains smaller than
$\sigma_{0}$ in gapless graphene. The gap dependence on the optical
conductivity is shown in Fig. 3b. First, the gap opening makes the chemical
potential bigger therefore the sharp structure in the $\sigma$ tends to larger
$\omega$ values. Second, the scattering mechanism increases by increasing the
electron densities and then the optical conductivity changes and becomes
smaller.

Another point of interest for experiments is the density dependence ( in units
of 1012 cm-2) of the optical conductivity ( Fig. 4) as a function of frequency
at $2\Delta=0.2$ eV. Note that the noninteracting chemical potential values
associated to the electron densities used in Fig. 4 are
$\mu_{0}=0.154,0.279,0.382$ and $0.831$ eV, respectively with giving
$\Delta=0.1$ eV. The optical conductivity increases by increasing the electron
density around $\omega_{0}$ however $\sigma$ decreases faster by increasing
the density at high frequency. The sharp structure of the optical conductance
tends to higher frequency by increasing the electron density. The sharp
position occurs at $2\mu$ which is always smaller than $2\mu_{0}$ for the same
system.

##  4 Conclusion

we have calculated the optical conductivity of gapped graphene, including the
effect of the lowest order self-energy diagram due to the electron-phonon
interaction by Holstein Hamiltonian. We have reported an extra logarithmic
singular behavior associated to gap value in the real part of the self-energy.
We have found the density, gap value and the electron-phonon coupling
dependence of the renormalized velocity and the interacting chemical
potential. The optical conductivity is affected by these physical quantities
and Pauli blocking below twice value of the renormalized chemical potential
and the gap values. We conclude that the inclusion of the renormalized Fermi
energy and the band gap affects are indeed crucial to get reasonable feature
for the optical conductivity. The gap dependence of the optical conductivity
would be verified by experiments.

##  5 acknowledgments

R. A thank S.G. Sharapov for stimulating discussion. We are grateful A.
Qaiumzadeh for useful comments. We thank Centro de Ciencias de Benasque, Spain
where this work was completed.

Note added- In final stage of preparing this manuscript, we became aware of a
related work for gapless graphene carbotte .

## References

  * (1) O. Gunnarsson, Rev. Mod. Phys. 69 (1997) 575 and references therein.
  * (2) K. S. Novoselov, A. K. Geim, S. V. Morozov, et.al., Science, 306 (2004) 666 .
  * (3) Akin Akturk and Neil Goldman, Journal of Applied physics, 103 (2008) 053702 .
  * (4) D. M. Basko and I. L. Aleiner, Phys. Rev. B 77 (2008) 041409 (R) .
  * (5) Cheol-Hwang Park, Feliciaon Giustino, Marivin L. Cohen and Steven G. Louie, Nano Letters 8 (2008) 4229 .
  * (6) T. Holstein, Ann. Phys. (N.Y.) 8, 325 (1959); 8 (1959) 343 .
  * (7) S. Piscanec et al., Phys. Rev. Lett. 93 (2004) 185503 .
  * (8) S. Piscanec et al., Phys. Rev. B 75 (2007) 035427 .
  * (9) S. Pisana et al. Nature Mater. 6 (2007) 198 .
  * (10) N. M. R. Peres, F. Guinea and A. H. Castro Neto, Phys. Rev. B 73(2006) 125411
  * (11) N. M. R. Peres, and T. Stauber, Int. J. Mod. Phys. B, 16 (2008) 2529 .
  * (12) T. Stauber, N. M. R. Peres, and A. K. Geim, Phys. Rev. B, 78 (2008) 085432 .
  * (13) V. P. Gusynin, S. G. Sharapov and J. P. Carbotte, Phys. Rev. Lett. 96(2006) 256802 .
  * (14) Z. Q. Li, et.al., Nat. Phys., 4(2008) 532 .
  * (15) R. R. Nair et. al., Science, 320 (2008) 1308 .
  * (16) T. Stauber, and N. M. R. Peres, J. Phys.: Condens. Matter 20 (2008) 055002.
  * (17) M. Calandra and F. Mauri, Phys. Rev. B 76 (2007) 205411 .
  * (18) W.-K. Tse and Das Sarma, Phys. Rev. Lett 99 (2007) 236802 .
  * (19) N. M. R. Peres, T. Stauber and A. H. Castro Neto, EPL 84 (2008) 38002 .
  * (20) T. Stauber, and N. M. R. Peres, Phys. Rev. B, 78(2008) 085418 .
  * (21) S.Y. Zhou et. al., Nature Mater., 76 (2007) 770 .
  * (22) K. S. Novoselov, et.al., nature 438(2005) 197 .
  * (23) E. Cappelluti and L. Benfatto, Phys. Rev. B 79 (2009) 035419 .
  * (24) A. Qaiumzadeh and R. Asgari, Phys. Rev. B 79 (2009) 075414 .
  * (25) A. H. Castro Neto et al., Rev. Mod. Phys., 81 (2009) 109 .
  * (26) Gerald D. Mahan, Many-Particle physics, (Plenum Peress, New York 1990) 2nd edition.
  * (27) M. Jonson, and G. D. Mahan, Phys. Rev. B, 21 (1980) 4223 .
  * (28) F. Dogan and F. Marsiglio, Phys. Rev. B 68 (2003) 165102 .
  * (29) M. Polini et al., Phys. Rev. B 77 (2008) 081411(R) .
  * (30) J. P. Carbotte, E. J. Nicol and S. G. Sharapov, arXiv:0908.2608 .

Figure 1: (Color online) Imaginary (a) and real (b) parts of the self-energy
as a function of energy evaluated at Fermi energy for different gap values at
the coupling constant $g=3.0$ and density $n=5.0\times 10^{12}$ cm-2.

Figure 2: (Color online) Renormalized electron velocity for (a) the different
values of coupling constants at $n=5.0\times 10^{12}$ cm-2, (b) the different
value of density (in units of $10^{12}$ cm-2) at $g=3.0$.

Figure 3: (Color online) Optical conductivity as a function of energy for (a)
the different values of coupling constant at $\Delta=0.1~{}eV$ and (b) the
different value of $\Delta$ at $g=3.0$. We consider $n=1\times
10^{13}~{}cm^{-2}$.

Figure 4: (Color online) Optical conductivity as a function of energy for
different values of density ( in unites of $10^{12}$ cm-2) at $g=3.0$ and
$\Delta=0.1~{}eV$