# Link overlap influences opinion dynamics on multiplex networks of Ashkin-
Teller spins

Cook Hyun Kim  CCSS, CTP and Department of Physics and Astronomy, Seoul
National University, Seoul 08826, Korea  Minjae Jo  CCSS, CTP and Department
of Physics and Astronomy, Seoul National University, Seoul 08826, Korea  J. S.
Lee  School of Physics, Korea Institute for Advanced Study, Seoul 02455, Korea
G. Bianconi  School of Mathematical Sciences, Queen Mary University of London,
E1 4GF, London, United Kingdom  Alan Turing Institute, The British Library,
NW1 2DB, London, United Kingdom  B. Kahng  Center for Complex Systems, KI of
Grid Modernization, Korea Institute of Energy Technology, Naju, Jeonnam 58217,
Korea

###### Abstract

Consider a multiplex network formed by two layers indicating social
interactions: the first layer is a friendship network and the second layer is
a network of business relations. In this duplex network each pair of
individuals can be connected in different ways: they can be connected by a
friendship but not connected by a business relation, they can be connected by
a business relation without being friends, or they can be simultaneously
friends and in a business relation. In the latter case we say that the links
in different layers overlap. These three types of connections are called
multilinks and the multidegree indicates the sum of multilinks of a given type
that are incident to a given node. Previous opinion models on multilayer
networks have mostly neglected the effect of link overlap. Here we show that
link overlap can have important effects in the formation of a majority
opinion. Indeed, the formation of a majority opinion can be significantly
influenced by the statistical properties of multilinks, and in particular by
the multidegree distribution. To quantitatively address this problem, we study
a simple spin model, called the Ashkin-Teller model including 2-body and
4-body interactions between nodes in different layers. Here we fully
investigate the rich phase diagram of this model which includes a large
variety of phase transitions. Indeed the phase diagram or the model displays
continuous, discontinuous, and hybrid phase transitions, and successive jumps
of the order parameters within the Baxter phase.

##  I INTRODUCTION

Over the past two decades, network theory NS ; Doro_book ; Newman_book ; jkps
; Doro_crit  has provided the pivotal framework for characterizing the
interplay between graph structures and dynamics of complex systems. Recently,
multilayer networks Bianconi2018multilayer ; Buldyrev2010 ; PhysicsReports ;
Kivela ; Goh_review  are attracting considerable scientific interest. These
network of networks are able to integrate information on various types of
links characterizing complex systems where interactions have different nature
and connotation. Therefore, they provide a useful perspective for analyzing
complex social, transportation, or biological systems Thurner ; Weighted ;
Boccaletti ; Bullmore2009 ; Makse  etc. Multilayer networks not only have rich
correlated structures PRE ; Vito_corr ; Goh ; Raissa  that encode more
information than a single layer, but also contain various dynamical processes
that are strongly affected by the multiplexity of the network. These dynamical
processes include percolation Buldyrev2010 ; Goh ; Baxter2012 ; Cellai1 ;
Algorithm ; BD1 ; Cellai2 , diffusion Arenas1 ; Arenas2 , epidemic spreading
Boguna_epidemics ; Cozzo_epidemics ; Arenas_aware , and game theory game1 ;
game2  etc.

Multiplex networks are a special class of multilayer network consisting of a
set of nodes connected by $M$ different types of links. Each network
consisting of a given type of link interaction forms one of the $M$ layers of
a multiplex network.

Most social networks are multiplex. In fact, social ties have different
connotations possibly indicating friends, colleagues, acquaintances and family
relations, etc. Moreover, in the modern society, online social interactions
can occur between different online social networks such as Twitter, Facebook
and LinkedIn etc. The vast majority of data on multiplex social networks
display a significant link overlap Weighted ; PRE ; Thurner . This property
indicates that a significant fraction of pair of nodes can be connected at the
same time by more than one type of interaction. For example, it might occur
that a colleague is also a friend or that a two individuals might be connected
at the same time in Facebook and Twitter.

The opinion dynamics on social multiplex networks have been investigated
recently using spin models such as the voter models Masuda ; Marina1 ;
Marina_Vito ; Chmiel1 ; Chmiel2 , election models Election , and Hamiltonian
spin systems Vito_spin . The observed dynamics on social multiplex networks
cannot be reduced to the dynamics on a single social aggregated network that
treats all the interactions of the multiplex network on an equal footing. In
adaptive voter models, an absorbing and shattered fragmentation transition
Marina1 ; Marina_Vito  occurs in which one layer can be fragmented into two
clusters each one reaching consensus on a different opinion, whereas the other
layer remains connected in one cluster. In election models, the competing
campaigns of two parties can give rise to election outcomes in which both
parties have a large electorate Election . Additionally, the party investing
more in building a connected network of supporters is more likely to win the
election Election . In studies of the opinion dynamics on multiplex networks,
where different opinion can be spread across different layers, an important
question is whether each node maintains coherent behavior, that is, has a
similar opinion in all the layers. A spin opinion model displaying a
coherence–incoherence transition was numerically investigated recently
Vito_spin . Spins are coupled within each layer to represent the interaction
between one node and its neighbors on a given topic, and also across layers to
represent the tendency of each node to take a coherent opinion on all the
topics.

Another spin model, which illustrates opinion dynamics in social networks due
to the influence of interdependence between different social communities, is
the Ashkin-Teller model Ashkin1943 . It was studied on scale-free (SF) network
in which the degree distribution follows power law and an analytical approach
revealed that a rich phase diagram including the critical end point was
obtained AT . It was considered on a duplex network with identical topology;
however, the most realistic multiplex network Buldyrev2010  can be the case in
which the layers of the bilayer network are distinct and the amount of overlap
is tunable.

Here our goal is to investigate to what extent link overlap affects the
opinion dynamics defined on multiplex networks and whether link overlap favors
coherent opinions.

We consider a duplex network formed by two layers where a two-state opinion
dynamics takes place. For example, one could consider a voting model for the
city council and for the national parliament. For each vote, nodes can be
influenced by a different set of nodes. In the previous example, the first
layer indicates the network influencing the city council vote, the second
layer indicates the network determining the national vote. The link overlap
has a clear effect on this opinion dynamics by coupling the two layers. In
fact if two nodes are connected in both layers it is natural to assume that
the simultaneous alignment of the opinions in both layers must be favored by
the dynamics. This considerations allow us to model the opinion dynamics in
presence of link overlap, with a spin Hamiltonian model that is a variation of
the Ashkin-Teller (AT) model Ashkin1943 ; AT  that we call $g$-AT model.

The model contains two species of Ising spins, the $s$-spin and $\sigma$-spin,
with each species of spin located on a single layer of the duplex network. The
duplex network is a maximum entropy duplex network with given multidegree
distribution PRE  and as such it is very suitable to modulate the role of
overlapping multilinks. In particular we here assume that non-overlapping
multilinks and overlapping multilinks have a SF multidegree distribution
characterized by a different power-law exponent. Here we provide a complete
analytical mean-field solution to this model and we reveal the complex phase
diagram of the model. We show that favoring the simultaneous alignment of the
opinions of nodes connected in both layers provides a simple mechanism to
generate coherence of opinions.

This paper is organized as follows: We introduce the Hamiltonian of the $g$-AT
model and the duplex network topology under study in Sec II. In Sec. III we
derive the free energy density using the mean-field approximation, and then
self-consistency equations for the order parameters by minimizing the free
energy density. Next, from these self-consistency equations, we obtain the
susceptibilities. In Secs. IV and V, we obtain rich phase diagrams in which
different phases in the parameter space are delimited by lines indicating
phase transitions (PTs) of different order. Note that the phase diagrams are
richer than those of the original AT model on SF networks AT , because the
links are classified into two types: non-overlapping and overlapping links.
Finally, we summarize the results in Sec. VI.

Figure 1: (Panel (a)) The $g$-AT model on a duplex network: two species
($s_{i}$, $\sigma_{i}$) of Ising spins describe respectively the opinion of
node $i$ in layer 1 and in layer 2. Each pair of nodes of the duplex network
can be connected by a different type of multilink: multilinks $(1,1)$ connect
pair of nodes in both layer 1 and layer 2; multilinks $(1,0)$ and $(0,1)$
connect pair of nodes only in layer 1 and only in layer 2 respectively.
Therefore multilinks $(1,1)$ describe overlapping links while multilinks
$(1,0)$ and $(0,1)$ describe non overlapping links. The model can be also
interpreted as a model on a colored network in which nodes are associated pair
of spin and the interactions between each pair of nodes can be distinguished
in mutlilinks $(1,1),(1,0)$ and $(0,1)$ (panel (b)).

Figure 2: The $g$-AT model is an Hamiltonian model combining two-body and
four-body interactions. The four body interactions characterizes the
interactions between the spins of species $s_{i}$, and of species $\sigma_{i}$
connected by a multilink $(1,1)$. The two-body interactions characterize the
coupling between spins of a given species (either the spins $s_{i}$ or the
spins $\sigma_{i}$) connected by either a multilink $(1,0)$ or $(0,1)$.

##  II MODEL AND FORMALISM

We consider a duplex network formed by $N$ nodes $i\in\\{1,2,\ldots,N\\}$.
Every pair of nodes $(i,j)$ of the duplex network can be connected in
multiplex ways. In order to indicate these different type of connections we
use multilinks introduced in Ref. PRE . In particular we say that a pair of
nodes $(i,j)$ is connected by a multilink $(1,0)$ if they are only connected
in layer 1, they are connected by a multilink $(0,1)$ if they are only
connected in layer 2, and they are connected by a multilink $(1,1)$ if they
are connected in both layers. Every pair of nodes can be connected only by one
type of multilink, alternatively they can be unconnected in both layers (see
Figure 1). We call multiadjacency matrices the matrices of elements
$A^{(1,0)}_{ij}$, $A^{(0,1)}_{ij}$, and $A^{(1,1)}_{ij}$ indicating whether or
not the pair of nodes $(i,j)$ is connected by a multilink $(1,0)$, a multilink
$(0,1)$, and a multilink $(1,1)$, respectively. This general duplex network
topology includes link overlap captured by the multilinks $(1,1)$. The
presence of such multilinks has been observed in a variety of social networks
Weighted ; Thurner . Here and in the following we indicate with multidegrees
$k_{i}^{(1,0)},k_{i}^{(0,1)}$, and $k_{i}^{(1,1)}$ the numbers of multilinks
incident to the node $i$, i.e.,

$\displaystyle k_{i}^{(1,0)}$ $\displaystyle=$
$\displaystyle\sum_{j=1}^{N}A_{ij}^{(1,0)},$ $\displaystyle k_{i}^{(1,0)}$
$\displaystyle=$ $\displaystyle\sum_{j=1}^{N}A_{ij}^{(1,0)},$ $\displaystyle
k_{i}^{(1,1)}$ $\displaystyle=$ $\displaystyle\sum_{j=1}^{N}A_{ij}^{(1,1)}.$
(1)

On such a duplex network, we consider the $g$-AT model that describes opinion
dynamics and takes into account the role that link overlap has on this
dynamics. We consider two species of Ising spins $s_{i}$ and $\sigma_{i}$
associated with the dynamics on layer 1 and layer 2, respectively. The two
spins take values $s_{i}\in\\{-1,1\\}$ and $\sigma\in\\{-1,1\\}$. These spin
variables are interacting via 2-body interactions and 4-body interactions (see
Figure 2). In particular, for each multilink $(1,0)$ connecting node $i$ to
node $j$, we have a 2-body Ising interaction between the spins $s_{i}$ and
$s_{j}$ with coupling constant $J_{2}$. Similarly, for each multilink $(0,1)$
connecting node $i$ to node $j$, we have a 2-body Ising interaction between
the spins $\sigma_{i}$ and $\sigma_{j}$ with coupling constant $J_{2}$. For
each multilink $(1,1)$ instead we consider a combination of 2-body and 4-body
interactions. The 2-body interactions tend to align spins associated to the
same layer with coupling constant $J_{2}$. The 4-body interactions couples
instead the four spins $s_{i}$, $s_{j}$, $\sigma_{i}$, and $\sigma_{j}$, and
is modulated by a coupling constant $J_{4}$. In particular, the Hamiltonian of
the $g$-AT model without an external magnetic field is expressed as the sum of
three terms,

$\displaystyle\mathcal{H}_{o}=\mathcal{H}^{(1,0)}+\mathcal{H}^{(0,1)}+\mathcal{H}^{(1,1)}$
(2)

where

$\displaystyle\mathcal{H}^{(1,0)}$ $\displaystyle=$ $\displaystyle-
J_{2}\sum_{\langle i,j\rangle}A_{ij}^{(1,0)}s_{i}s_{j}$
$\displaystyle\mathcal{H}^{(0,1)}$ $\displaystyle=$ $\displaystyle-
J_{2}\sum_{\langle i,j\rangle}A_{ij}^{(0,1)}\sigma_{i}\sigma_{j}$
$\displaystyle\mathcal{H}^{(1,1)}$ $\displaystyle=$ $\displaystyle-
J_{4}\sum_{\langle i,j\rangle}A_{ij}^{(1,1)}s_{i}s_{j}\sigma_{i}\sigma_{j}$
(3) $\displaystyle-J_{2}\sum_{\langle
i,j\rangle}A_{ij}^{(1,0)}s_{i}s_{j}-J_{2}\sum_{\langle
i,j\rangle}A_{ij}^{(0,1)}\sigma_{i}\sigma_{j}$

with the pairs of connected nodes $\langle i,j\rangle$. Alternatively, the
Hamiltonian $\mathcal{H}_{o}$ of the $g$-AT model without an external magnetic
field can be expressed more concisely as

$\displaystyle-\beta\mathcal{H}_{o}=$ $\displaystyle K_{2}\sum_{\langle
i,j\rangle}\mathbf{s}_{i}^{\mathbf{T}}\textbf{A}_{ij}\mathbf{s}_{j}\,,$ (4)

where $\mathbf{s}_{i}=(s_{i},\sigma_{i},s_{i}\sigma_{i})^{\mathbf{T}}$, and
the matrix $\textbf{A}_{ij}$ is given by

${\bf A}_{ij}=\left(\begin{array}[]{ccc}A^{(1,0)}_{ij}+A^{(1,1)}_{ij}&0&0\\\
0&A^{(0,1)}_{ij}+A^{(1,1)}_{ij}&0\\\ 0&0&xA^{(1,1)}_{ij}\end{array}\right),$
(5)

where $x\equiv J_{4}/J_{2}$. Moreover, $\beta=1/k_{\rm B}T$, where $k_{\rm B}$
is the Boltzmann constant, $T$ is the temperature, and $K_{2}\equiv\beta
J_{2}$ with coupling constant $J_{2}$. For later discussion, we define
similarly $K_{4}\equiv\beta J_{4}$.

Here, we investigate the critical properties of this model on a maximum
entropy duplex network model with given multidegree distribution PRE . In
order to distinguish between multilinks $(1,1)$ which imply link overlap and
the other multilinks $(1,0)$ and $(0,1)$ which do not, we assume for
simplicity that each node $i$ of the multiplex network has the same
multidegree $(1,0)$ and multidegree $(0,1)$, and we indicate the multidegree
of non-overlapping multilinks and of overlapping multilinks as

$\displaystyle k_{i}^{(1,0)}$ $\displaystyle=$ $\displaystyle
k_{i}^{(0,1)}=k_{\textrm{n},i}\,,$ $\displaystyle k_{i}^{(1,1)}$
$\displaystyle=$ $\displaystyle k_{\textrm{o},i}.$ (6)

where the subscript n of $k_{\textrm{n},i}$ indicates and the subscript o of
$k_{\textrm{o},i}$ indicates non-overlap and overlapping multilinks,
respectively. We assume that the degree distributions corresponding to
overlapping and non-overlapping multilinks are power-law functions with
exponents $\lambda_{\textrm{o}}$ and $\lambda_{\textrm{n}}$, respectively. The
degree distribution is shortly written as

$\displaystyle P_{d}(k_{a})\sim k_{a}^{-\lambda_{a}},$ (7)

where $a\in\\{\textrm{o},\textrm{n}\\}$.

In the considered ensemble of duplex networks PRE  a pair of nodes $(i,j)$ is
connected by $(1,0)$ multilinks with probability $p^{(1,0)}_{ij}$, by $(0,1)$
multilinks with probability $p^{(0,1)}_{ij}$, and by $(1,1)$ multilinks with
probability $p^{(1,1)}_{ij}$, where we have

$\displaystyle p^{(1,0)}_{ij}$ $\displaystyle=$
$\displaystyle\frac{k_{i}^{(1,0)}k_{j}^{(1,0)}}{\left\langle{k^{(1,0)}}\right\rangle
N},$ $\displaystyle p^{(0,1)}_{ij}$ $\displaystyle=$
$\displaystyle\frac{k_{i}^{(0,1)}k_{j}^{(0,1)}}{\left\langle{k^{(0,1)}}\right\rangle
N},$ $\displaystyle p^{(1,1)}_{ij}$ $\displaystyle=$
$\displaystyle\frac{k_{i}^{(1,1)}k_{j}^{(1,1)}}{\left\langle{k^{(1,1)}}\right\rangle
N}$ (8)

with $\langle k^{(1,0)}\rangle$, $\langle k^{(0,1)}\rangle$, and $\langle
k^{(1,1)}\rangle$ being the average multidegrees. Indeed these marginal
probabilities are obtained in the maximum entropy ensemble with given
multidegree distribution as long as the degree distribution display the
structural cutoff. Here, we consider the thermodynamic limit ($N\to\infty$)
and power-law exponents greater than $3$, so that the effect of structural
cutoff can be ignored.

The phase diagram of this model will be affected by the topology of multiplex
network and the strength of the interlayer interaction. This can be studied as
a function of three parameters, $\lambda_{\textrm{n}}$,
$\lambda_{\textrm{o}}$, and $x\equiv J_{4}/J_{2}$. The ratio $x$ quantifies
the degree strength of 4-body interaction with respect to the strength of
2-body interactions for between nodes linked by multilinks $(1,1)$.

The original AT model AT  comprises two species of Ising spins, $s_{i}$ and
$\sigma_{i}$, locating at each node $i$ on a monolayer network. The original
AT-model can be thus recovered as a limit case of the $g$-AT model in absence
of non-overlapping multilinks and when $x=1$ (i.e. $J_{4}=J_{2}$). Indeed in
this limit we recover the Hamiltonian for the original AT model given by

$\displaystyle\mathcal{H}=-J_{2}\sum_{\langle
ij\rangle}(s_{i}s_{j}+\sigma_{i}\sigma_{j})-J_{2}\sum_{\langle
ij\rangle}s_{i}s_{j}\sigma_{i}\sigma_{j},$ (9)

that can be rewritten in the form of the 4-state Potts model as

$\displaystyle\mathcal{H}=-4J_{2}\sum_{\langle
ij\rangle}\left(\delta(q_{i},q_{j})-1/4\right),$ (10)

where $q_{i}$ is a Potts spin with value $0,1,2$ or $3$ at node $i$ and
$\delta(q_{i},q_{j})=1$ for $q_{i}=q_{j}$, and zero otherwise Kadanoff . Since
the non-overlapping multilinks are absent, the phase diagram of the original
model is a function of a single power-law exponent $\lambda$ of the degree
distribution. Clearly this power-law exponent correspond to the power-law
exponent $\lambda_{\textrm{o}}$ of overlapping links of the $g$-AT model.

##  III Mean-field solution

To obtain the Landau free energy, we calculate the Hamiltonian in Eq. (4) by
the mean-field approximation. We first take the local order parameters
$\mathbf{m}_{i}=(m^{s}_{i},m^{\sigma}_{i},m^{s\sigma}_{i})^{\mathbf{T}}$,
whose components are defined as $m^{s}_{i}=\langle s_{i}\rangle$,
$m^{\sigma}_{i}=\langle\sigma_{i}\rangle$, and $m^{s\sigma}_{i}=\langle
s_{i}\sigma_{i}\rangle$. Here $\langle\cdots\rangle$ is the ensemble average
of a given quantity. Next, we expand each spin variable with respect to the
respective local order parameter as $\mathbf{s}_{i}=(m^{s}_{i}+\delta
m^{s}_{i},m^{\sigma}_{i}+\delta m^{\sigma}_{i},m^{s\sigma}_{i}+\delta
m^{s\sigma}_{i})^{\mathbf{T}}$. We can neglect the higher-order terms in
$\delta m^{s}_{i}$, $\delta m^{\sigma}_{i}$, and $\delta m^{s\sigma}_{i}$
because the magnitude of these terms is very small compared to that of the
local order parameter. The mean-field Hamiltonian $\mathcal{H}_{\textrm{mf}}$
can be written as

$\displaystyle-\beta\mathcal{H}_{\textrm{mf}}\simeq$ $\displaystyle-
K_{2}\sum_{i,j}\mathbf{m}_{i}^{\mathbf{T}}\textbf{A}_{ij}\mathbf{m}_{j}+K_{2}\sum_{i,j}\mathbf{m}_{i}^{\mathbf{T}}\textbf{A}_{ij}(\mathbf{s}_{j}+\mathbf{\sigma}_{j})\,.$
(11)

Then, we obtain the mean-field Landau free energy $\mathcal{F}$, which is
given by

$\displaystyle\beta\mathcal{F}$ $\displaystyle=-\ln Z$ (12)
$\displaystyle=-\ln\sum_{\\{s_{i},\sigma_{i}\\}}e^{-\beta\mathcal{H}_{\rm
mf}}\simeq-\sum_{i}\ln
Z_{i}+K_{2}\sum_{i,j}\mathbf{m}_{i}^{\mathbf{T}}\textbf{A}_{ij}\mathbf{m}_{j},$
(13)

where

$\displaystyle
Z_{i}=4\left[\mathcal{C}_{i}(s)\mathcal{C}_{i}(\sigma)\mathcal{C}_{i}(s\sigma)+\mathcal{S}_{i}(s)\mathcal{S}_{i}(\sigma)\mathcal{S}_{i}(s\sigma)\right],$
(14)

with

$\displaystyle\mathcal{C}_{i}(s)\equiv\cosh\left(\sum_{j\in\text{nn}(i)}K_{2}m^{s}_{j}\right),\quad\mathcal{S}_{i}(s)\equiv\sinh\left(\sum_{j\in\text{nn}(i)}K_{2}m^{s}_{j}\right)\,.$
(15)

Here $\sum_{j\in\text{nn}(i)}$ indicates that the summation runs over all the
nearest neighbors $j$ of node $i$ for each of the three types of links.

Next, we use the annealed approximation to perform the summation:

$\displaystyle\left.\begin{array}[]{ll}\sum_{\langle
i,j\rangle}A^{(1,0)}_{ij}\mathcal{A}_{ij}\rightarrow\frac{1}{2}\sum_{i,j}p_{ij}^{(1,0)}\mathcal{A}_{ij}\,,\\\
\\\ \sum_{\langle
i,j\rangle}A^{(0,1)}_{ij}\mathcal{A}_{ij}\rightarrow\frac{1}{2}\sum_{i,j}p_{ij}^{(0,1)}\mathcal{A}_{ij}\,\textrm{
and}\\\ \\\ \sum_{\langle
i,j\rangle}A^{(1,1)}_{ij}\mathcal{A}_{ij}\rightarrow\frac{1}{2}\sum_{i,j}p_{ij}^{(1,1)}\mathcal{A}_{ij}\,,\end{array}\right.$
(21)

where $\mathcal{A}_{ij}$ is a given function of $i$ and $j$ and
$p_{ij}^{(1,0)},p_{ij}^{(0,1)}$ and $p_{ij}^{(1,1)}$ are defined in Eq. (8).

We define a global order magnetization for $s$ spin :

$\displaystyle m^{(1,0)}_{s}=\dfrac{\sum_{i}k_{i}^{(1,0)}m^{s}_{i}}{N\langle
k^{(1,0)}\rangle}\textrm{ and
}m^{(1,1)}_{s}=\dfrac{\sum_{i}k_{i}^{(1,1)}m^{s}_{i}}{N\langle
k^{(1,1)}\rangle},$ (22)

where $m^{s}_{i}$ is the local order parameter for $s$ spin. We introduce
global order parameters for $\sigma$ and $s\sigma$ spins similarly. Then, we
set that $M\equiv m^{(1,1)}_{s\sigma}$.

Since the considered duplex network ensemble has the same multidegree
distribution of the non-overlapping multilinks we can set

$\displaystyle m_{s}^{(1,0)}=m_{\sigma}^{(0,1)}\equiv
m_{\textrm{n}}\,,\qquad\textrm{ }m_{s}^{(1,1)}=m_{\sigma}^{(1,1)}\equiv
m_{\textrm{o}}\,.$ (23)

The three order parameters are now denoted as $m_{\textrm{o}}$,
$m_{\textrm{n}}$, and $M$-magnetization, respectively. Applying the annealed
approximation, we rewrite the free energy density
($f\equiv\beta\mathcal{F}/N$) in terms of the order parameters
$m_{\textrm{n}}$, $m_{\textrm{o}}$, and $M$. The free energy density $f$ is
given by

$\displaystyle f\simeq K_{2}m_{\textrm{n}}^{2}\langle
k_{\textrm{n}}\rangle+K_{2}m_{\textrm{o}}^{2}\langle
k_{\textrm{o}}\rangle+\dfrac{1}{2}K_{4}M^{2}\langle k_{\textrm{o}}\rangle$
(24) $\displaystyle~{}-2\int_{k_{\rm min}^{\textrm{n}}}^{\infty}\int_{k_{\rm
min}^{\textrm{o}}}^{\infty}dk_{\textrm{n}}dk_{\textrm{o}}P_{d}(k_{\textrm{n}})P_{d}(k_{\textrm{o}})\ln\left[\cosh\left(K_{2}(m_{\textrm{n}}k_{\textrm{n}}+m_{\textrm{o}}k_{\textrm{o}})\right)\right]$
(25) $\displaystyle~{}-\int_{k_{\rm
min}^{\textrm{o}}}^{\infty}dk_{\textrm{o}}P_{d}(k_{\textrm{o}})\ln\left[\cosh\left(K_{4}Mk_{\textrm{o}}\right)\right]-\mathcal{B}_{1}\,,$
(26)

where $K_{4}=\beta J_{4}$ with coupling constant $J_{4}$ and

$\displaystyle\mathcal{B}_{1}$ $\displaystyle=\int_{k_{\rm
min}^{\textrm{n}}}^{\infty}\int_{k_{\rm
min}^{\textrm{o}}}^{\infty}dk_{\textrm{n}}dk_{\textrm{o}}P_{d}(k_{\textrm{n}})P_{d}(k_{\textrm{o}})\ln\left(1+\mathcal{T}_{2}^{2}\mathcal{T}_{4}\right)\,,$
(27)

with

$\displaystyle\mathcal{T}_{2}\equiv\tanh\left(K_{2}(m_{\textrm{n}}k_{\textrm{n}}+m_{\textrm{o}}k_{\textrm{o}})\right),\,\mathcal{T}_{4}\equiv\tanh\left(K_{4}Mk_{\textrm{o}}\right)\,.$
(28)

Minimizing the free energy density $f$, ${\partial f}/{\partial m_{a}}=0$ and
${\partial f}/{\partial M}=0$, we obtain the following self-consistency
relations:

$\displaystyle m_{a}\langle k_{a}\rangle$ $\displaystyle=\int_{k_{\rm
min}^{\textrm{n}}}^{\infty}\int_{k_{\rm
min}^{\textrm{o}}}^{\infty}dk_{\textrm{n}}dk_{\textrm{o}}P_{d}(k_{\textrm{n}})P_{d}(k_{\textrm{o}})\dfrac{\mathcal{T}_{2}\left(1+\mathcal{T}_{4}\right)}{1+\mathcal{T}_{2}^{2}\mathcal{T}_{4}}k_{a}$
(29)

where $a\in\\{\textrm{o},\textrm{n}\\}$, and

$\displaystyle M\langle k_{\textrm{o}}\rangle$ $\displaystyle=\int_{k_{\rm
min}^{\textrm{n}}}^{\infty}\int_{k_{\rm
min}^{\textrm{o}}}^{\infty}dk_{\textrm{n}}dk_{\textrm{o}}P_{d}(k_{\textrm{n}})P_{d}(k_{\textrm{o}})\dfrac{\mathcal{T}_{4}+\mathcal{T}_{2}^{2}}{1+\mathcal{T}_{2}^{2}\mathcal{T}_{4}}k_{\textrm{o}}.$
(30)

The self-consistency relations (29) and (30) admit three solutions,
corresponding to the paramagnetic phase ($m_{a}=0,M=0$), the Baxter phase
($m_{a}>0,M>0$), and the $\langle\sigma s\rangle$ phase ($m_{a}=0,M>0$).

To obtain the susceptibility, we also consider a Hamiltonian including an
external magnetic field, given by

$\displaystyle-\beta\mathcal{H}$
$\displaystyle=-\beta\mathcal{H}_{o}+\sum_{i}\mathbf{H}_{i}^{\mathbf{T}}\mathbf{s}_{i}\,,$
(31)

where
$\textbf{H}_{i}=(k_{\textrm{o},i}H_{\textrm{o}}+k_{\textrm{n},i}H_{\textrm{n}},k_{\textrm{o},i}H_{\textrm{o}}+k_{\textrm{n},i}H_{\textrm{n}},k_{\textrm{o},i}H_{4})^{\mathbf{T}}$.
$H_{a}$ is the external magnetic field applied to $s$ and $\sigma$ spins in
proportion to the multidegree $k_{\textrm{a}}$ and $H_{4}$ is another external
magnetic field applied to $s\sigma$ spins in proportion to degree
$k_{\textrm{o}}$. Minimizing the free energy density, we obtain the self-
consistency equations for magnetizations with respect to external magnetic
fields:

$-\dfrac{\partial f}{\partial H_{a}}=m_{a}\langle
k_{a}\rangle,\quad-\dfrac{\partial f}{\partial H_{4}}=M\langle
k_{\textrm{o}}\rangle\,.$ (32)

These self-consistency equations can be obtained by substituting
$K_{2}m_{a}k_{a}$ and $K_{4}Mk_{\textrm{o}}$ with $(K_{2}m_{a}+H_{a})k_{a}$
and $(K_{4}M+H_{4})k_{\textrm{o}}$, respectively, in Eqs. (26), (27), (29),
and (30) (see Appendix A). The susceptibilities are calculated using the
following relations:

$\chi_{a}\equiv\dfrac{\partial m_{a}}{\partial H_{a}}\Big{|}_{H_{a},H_{4}\to
0},\quad\chi_{M}\equiv\dfrac{\partial M}{\partial
H_{4}}\Big{|}_{H_{a},H_{4}\to 0}\,.$ (33)

Using the above relations, the susceptibilities $\chi$ can be obtained as
follows:

$\displaystyle\chi_{a}$
$\displaystyle=\frac{\mathcal{A}_{aa}+\mathcal{A}_{a\bar{a}}K_{2}\partial
m_{\bar{a}}/\partial H_{a}+\mathcal{A}_{aM}K_{4}\partial M/\partial
H_{a}}{\langle k_{a}\rangle-K_{2}\mathcal{A}_{aa}}\,,$ (34)
$\displaystyle\chi_{M}$
$\displaystyle=\frac{\mathcal{A}_{MM}+\mathcal{A}_{M\textrm{o}}K_{2}\partial
m_{\textrm{o}}/\partial H_{4}+\mathcal{A}_{M\textrm{n}}K_{2}\partial
m_{\textrm{n}}/\partial H_{4}}{\langle k_{\textrm{o}}\rangle-
K_{4}\mathcal{A}_{MM}}\,,$ (35)

where $a\in\\{\textrm{n},\textrm{o}\\}$ and
$\bar{a}\in\\{\textrm{n},\textrm{o}\\}$ with $a$ different from $\bar{a}$.
Here the $\mathcal{A}$ terms are obtained as follows

$\displaystyle\mathcal{A}_{aa}$ $\displaystyle=\dfrac{\partial
m_{a\textrm{I}}}{\partial H_{a}},\quad\,\mathcal{A}_{a\bar{a}}=\dfrac{\partial
m_{a\textrm{I}}}{\partial H_{\bar{a}}},\quad\,\mathcal{A}_{aM}=\dfrac{\partial
m_{a\textrm{I}}}{\partial H_{4}},\,$ (36) $\displaystyle\mathcal{A}_{Ma}$
$\displaystyle=\dfrac{\partial M_{\textrm{I}}}{\partial
H_{a}},\quad\hbox{and}\quad\mathcal{A}_{MM}=\dfrac{\partial
M_{\textrm{I}}}{\partial H_{4}}.$ (37)

In Appendix B, we provide the extensive formulas for $\mathcal{A}$s in the
limit $H_{a}\to 0$ and $H_{4}\to 0$, where $m_{a\textrm{I}}$ and
$M_{\textrm{I}}$ are presented in intergral form in Eqs. (61) and Eqs. (62).

##  IV Phase diagram I: $J_{4}/J_{2}$-DEPENDENCE

###  IV.1 Phases of the model

The $g$-AT model has one of the three phases, paramagnetic phase, Baxter
phase, and $\langle\sigma s\rangle$ phase, in equlibrium state, depending on
$\lambda_{\textrm{n}}$, $\lambda_{\textrm{o}}$, $x=J_{4}/J_{2}$, and $T$.

  * –

The paramagnetic phase is characterized by the order parameters
$m_{a}=\langle{s}\rangle=\left\langle{\sigma}\right\rangle=0,M=\langle{\sigma
s}\rangle=0$. This is the characteristic phase found in the high temperature
region, where the stochastic element of the dynamics is dominant. This phase
corresponds to an equilibrium configuration in which there is no majority
opinion in either layer ($m_{a}=0$), and each node has a random and
uncorrelated opinion in the two different layers ($M=0$). Therefore this is
the phase entirely dominated by noise.

  * –

The Baxter phase is characterized by the order parameters
$m_{a}=\langle{s}\rangle=\langle{\sigma}\rangle>0$, and $M=\langle{\sigma
s}\rangle>0$. This is the phase in which we observe the formation of a
majority opinion which is the same in both layers ($m_{a}>0$). Therefore, each
node has coherent opinions in the two distinct layers ($M>0$).

  * –

The $\langle{\sigma s}\rangle$ phase (Coherent phase) is characterized by the
order parameters $m_{a}=\langle{s}\rangle=\langle{\sigma}\rangle=0$, and
$M=\langle{\sigma s}\rangle>0$. This phase occurs for high temperature and
when $x=J_{4}/J_{2}$ is sufficiently high, in which the 4-body interactions
are stronger than the 2-body interactions; therefore each single node of the
multiplex network tends to have the same opinion in both layers but these
opinions are not yet aligned with the opinion of their neighbours. As a
consequence there is no yet formation of a majority opinion in each layer
(i.e. $m_{a}=0$). Note that the term of $\langle{\sigma s}\rangle$ phase
originates from the original paper in physics Ashkin1943 . To impose a meaning
on the phase in the perspective of opinion formation, we call it Coherent
phase hereafter.

###  IV.2 Classification of critical points and regions of the phase diagram

In the $g$-AT model the transitions between the phases Para, Baxter and
Coherent occurs as a function of the temperature $T$ and very diverse and rich
critical phenomena are observed. Indeed the PTs can be continuous,
discontinuous, hybrid and in general we can observe more than one PT as the
temperature $T$ is lowered, while the other parameters are kept unchanged.

To be concrete we discuss here an exemplar phase diagram of the $g$-AT model
in the parameter space [$x,T^{-1}$] (see Fig. 3). This phase digram is
obtained for the power-law exponents $\lambda_{\textrm{n}}=3.53$ and
$\lambda_{\textrm{o}}=3.90$. For this value of the power-law exponents, the
phase diagram is similar to that of the original AT model AT  in the range
$\lambda_{c}<\lambda_{\textrm{o}}<4$, where $\lambda_{c}\approx 3.503$
indicates the tricritical point (TP) of the original AT model AT . In
particular we recall that in the original AT-model for
$\lambda_{\textrm{o}}>\lambda_{c}$, the PT is of the first-order; otherwise,
it is of the second-order AT ; Potts ; Potts_mendes . Here we will describe in
detail this phase diagram while the dependence of the phase diagram on the
power-law exponents $\lambda_{\textrm{n}}$ and $\lambda_{\textrm{o}}$ will be
treated in the next section.

In the phase diagram shown in Fig. 3, the three phases of the dynamics are
denoted by Para ($m_{a}=0,M=0$), Baxter phase ($m_{a}>0,M>0$), and
$\langle\sigma s\rangle$ phase ($m_{a}=0,M>0$). Dotted and solid lines
represent discontinuous and continuous PTs, respectively. The critical
temperature $T_{s}$ denotes the temperature at which a second-order PT occurs
from the Para phase to the Baxter phase. Note that $T_{s}$ is independent on
$x=J_{4}/J_{2}$ for $x<x_{e}$.

The phase diagram has characteristic points denoted as GZs and CEs. We
indicate with GZ a point at which the jump size (gap) of the order parameter
becomes zero at each side of the dotted curve. We indicate with CE a critical
endpoint, locating at the end of a continuous PT line, at which a line of
first-order PT and a line of discontinuity of the order parameter merge. We
will show that a mixed-order (or hybrid) transition occurs at these CE points.
There are two GZs and three CEs in Fig. 3. Their $x$ positions are asymmetric.

Figure 3: Schematic phase diagram of the $g$-AT model for a given set of
$\lambda_{\textrm{n}}=3.53$ and $\lambda_{\textrm{o}}=3.90$. Solid and dotted
curves represent continuous and discontinuous PTs, respectively. This phase
diagram is mostly similar to the one of the original AT model AT .

In the phase diagram shown in Fig. 3 we distinguish six regions based on the
ratio $x=J_{4}/J_{2}$:

  * –

In region (i), as the temperature is lowered the system undergoes a continuous
PT at $T_{s}$ from the Para to the Baxter phase. Therefore as the noise is
reduced the system goes continuously from a phase with the absence of any
order, to a state with a clear majority opinion which is the same in both
layers. This transition is denoted as (i)-type PT.

  * –

In region (ii), a continuous PT occurs at $T_{s}$ between the Para and the
Baxter phase. As $T$ is lowered further, a discontinuous jump of the order
parameters $m_{a}$ and $M$ occurs subsequently at $T_{f}$, in which we observe
a discontinuity in $m_{a}$ and $M$ between two non-zero values. This indicates
that at $T_{f}$, there is a discontinuous increment in the fraction of nodes
adopting the majority opinion. This transition is denoted as (ii)-type PT.

  * –

In region (iii), a discontinuous PT occurs at $T_{f}$ between the Para and the
Baxter phase. This implies that a majority opinion is formed abruptly in both
layers. This transition is denoted as (iii)-type PT.

  * –

In region (iv), a continuous PT occurs between the Para and the Coherent phase
at $T_{s,M}$. As $T$ is decreased further, a discontinuous PT occurs at
$T_{f}$ from the Coherent phase to the Baxter phase. This implies that as the
noise is reduced, at temperatures below the first continuous transition each
single node tend to adopt a coherent opinion in both layers, and then when the
temperature is further reduced a majority opinion is reached abruptly in both
layers. These transitions are denoted as (iv)-type PT.

  * –

In region (v), two continuous PTs occur successively: between the Para and the
Coherent phase at $T_{s,M}$ and between the Coherent and the Baxter phase at
$T_{s}^{\prime}$, respectively. Then as the temperature is decreased further,
the order parameters $m_{a}$ and $M$ jumps at $T_{f}$ from one finite value to
another. These transitions are denoted as (v)-type PT.

  * –

In region (vi), two continuous PTs occur between the Para and the Coherent
phase at $T_{s,M}$ and between the Coherent and the Baxter phase at
$T_{s}^{\prime}$. These transitions are denoted as (vi)-type PT.

###  IV.3 Free-energy landscape for the $x$-dependence of phase transitions

In this paragraph we will discuss the critical behavior of the $g$-AT model as
a function of the parameter $x=J_{4}/J_{2}$.

The phase diagram of the $g$-AT model can be treated separately for $x<x_{M}$
and $x>x_{M}$, where $x_{M}$ indicates the characteristic ratio between
$J_{4}$ and $J_{2}$. For $x<x_{M}$, the $J_{2}$ interactions are dominant, and
$O(m_{a})\gg O(M)$ near the transition temperature. On the other hand, for
$x>x_{M}$, $J_{4}$ interactions (interlayer interaction) become dominant and
$O(m_{a})\ll O(M)$ near the transition temperature. Thus, the Coherent phase
can emerge. For the original AT model, $x_{M}=1$; however, for the $g$-AT
model, $x_{M}$ depends on $\lambda_{\textrm{n}}$ and $\lambda_{\textrm{o}}$.
$x_{M}$ locates between $x_{e}$ and $x_{e,M}$ in Fig. 3. Explicitly formula to
derive $x_{M}$ will be presented in Eq. (84).

The free-energy landscape determines the location and type of PTs with respect
to $x$. Here we provide the discussion of the main results obtained by
investigating the properties of the free-energy density illustrated by in
Figs. 4 and 5 near $T_{s}$. We refer the interested reader to the exact
formula of the free-energy density $f$ given in Appendix E.

Figure 4: (a)–(d) Plot of the order parameters $m_{a}$ and $M$ as a function
of $T/J_{2}$, and (e)–(h) Plot of the free energy density landscape as a
function of $m_{\textrm{n}}$ for $\lambda_{\textrm{n}}=3.53$ and
$\lambda_{\textrm{o}}=3.90$ and various interlayer interaction ratios:
$x=1.30$ for (a) and (e); $x=1.40$ for (b) and (f); $x=1.62$ for (c) and (g);
and $x=1.80$ for (d) and (h). The transition types are second-order in region
(i) (a) and (e); successive continuous-discontinuous in region (ii) for (b)
and (f); mixed-order at CE1 for (c) and (g); and discontinuous transition in
region (iii) for (d) and (h).

####  IV.3.1 Case $0<x<x_{M}$

In this paragraph we describe the critical behavior of the $g$-AT model for
$x<x_{M}$ including regions (i), (ii) and (iii), and a point CE1. In region
(i) a continuous PT occurs between the Para and the Baxter phase. Therefore as
the temperature is lowered, both $m_{a}$ and $M$ increase continuously for
$T<T_{s}$.

Using Eqs. (87) and (90), we can obtain the critical behavior of the order
parameters for $T<T_{s}$

$\displaystyle m_{a}$ $\displaystyle\sim\,(T_{s}-T)^{\beta_{m}}$ with
$\displaystyle\beta_{m}$ $\displaystyle=\dfrac{1}{\lambda_{\textrm{min}}-3},$
(38) $\displaystyle M$ $\displaystyle\sim\,(T_{s}-T)^{\beta_{M}}$ with
$\displaystyle\beta_{M}$
$\displaystyle=\dfrac{\lambda_{\textrm{o}}-2}{\lambda_{\textrm{min}}-3},$ (39)

where
$\lambda_{\textrm{min}}=\textrm{min}(\lambda_{\textrm{o}},\lambda_{\textrm{n}})$.
The specific heat scales as

$\displaystyle C\,\sim(T_{s}-T)^{-\alpha}$ with
$\displaystyle\alpha=\dfrac{\lambda_{\textrm{min}}-5}{\lambda_{\textrm{min}}-3}.$
(40)

The susceptibility diverges as (see Appendix E for the derivation)

$\displaystyle\chi_{a}\sim\left\\{\begin{array}[]{cccc}(T_{s}-T)^{-\gamma^{-}}&\quad\hbox{with}\quad\gamma^{-}=1&\hbox{for}&T<T_{s},\cr\cr(T-T_{s})^{-\gamma^{+}}&\quad\hbox{with}\quad\gamma^{+}=1&\hbox{for}&T>T_{s}.\end{array}\right.$
(44)

Secondly we observe that as $x$ is increased but still remains less than
$x_{M}$, a jump arises in the order parameters $m_{a}$ and $M$ in region (ii),
observed for $x_{g}<x<x_{e}$ in Fig. 3. We observe that the system undergoes a
continuous second order transition at $T_{s}$ between the Para and the Baxter
phase characterized by the same critical exponents listed above. Moreover as
the temperature $T$ is lowered further the system undergoes a sudden increase
of the order parameter at $T_{f}$. Indeed at $T_{f}$ the free energy density
$f$ displays a global minimum at a finite $m_{a}$, leading to the abrupt
change of the order parameters (see Fig. 4(b) and (f)).

For $x\to x_{e}$, the temperature $T_{s}$ becomes equal to $T_{f}$. Therefore,
two global minima of $f(m_{a})$ occur at $m_{a}=0$ and $m_{a}>0$,
simultaneously. At this point, the second-order and the first-order transition
lines merge. Therefore, the critical behavior appears, together with the jump
of the order parameters $m_{a}$ and $M$ as illustrated in Fig. 4 (c) and (g).
This type of PT is referred to as a mixed-order (or hybird) transition and
this point is named as critical endpoint. However, the susceptibility
$\chi_{m}$ diverges at $T_{s}^{+}$ as it appears in a continuous PT.

In region (iii), for $x_{e}<x<x_{M}$, a single discontinuous transition occurs
at $T_{f}$ between the Para and the Baxter phase. The critical bahavior in
this region is illustrated in Fig. 4(d) and (h).

Figure 5: (a)–(d) Plot of the order parameters $m_{a}$ and $M$ as a function
of $T/J_{2}$. (e)–(h) Plot of the free energy density landscape as a function
of $m_{\textrm{n}}$ for $\lambda_{\textrm{n}}=3.53$ and
$\lambda_{\textrm{o}}=3.90$ and various interlayer interaction ratios:
$x=2.10$ for (a) and (e); $x=2.30$ for (b) and (f); $x=2.45$ for (c) and (g);
and $x=2.60$ for (d) and (h). The transition types are first-order in regime
(iii) for (a) and (e); successive continuous-discontinuous in regime (iv) for
(b) and (f); successive continuous-discontinuous in regime (v) for (c) and
(g); and continuous transition in regime (vi) for (d) and (h).

####  IV.3.2 Case $x\geq x_{M}$

Here, we consider the critical behavior of the $g$-AT model for $x\geq x_{M}$
including the regions (iii), (iv), (v) and (vi) and two critical endpoints
(CE2, CE3).

At $x=x_{M}$ in region (iii), the free-energy densities $f(m_{a})$ and $f(M)$
develop a local minimum at $m_{a}^{*}$ larger than 0 and at a temperature
$T_{f}$ larger than $T_{s,M}$ (Fig. 5(a) and (e)). A discontinuous transition
of $m_{a}$ and $M$ occurs between the Para and the Baxter phase at $T_{f}$.
This phenomenology remains unchanged for $x_{M}<x<x_{e,M}$.

At $x=x_{e,M}$, the system is at the boundary between the regions (iii) and
(iv), and we find the critical endpoint denoted as CE2 in Fig. 3. At CE2, when
the temperature $T$ approaches $T_{s,M}$ from below, i.e., when we explore the
critical behavior for $T\to T_{s,M}^{-}$, we observe a discontinuity in the
value of $M$ from a non-zero value to zero, and thus the susceptibility does
not diverge. However, when $T$ approaches $T_{s,M}$ from above, i.e., when we
explore the critical behavior for $T\to T_{s,M}^{+}$ even if $M$ jumps
suddenly and shows a behavior reminiscent of a first-order transition, the
susceptibility $\chi_{M}$ diverges. This is due to the fact that CE2 is the
endpoint of a line of second-order PTs between the Para and the Coherent
phase. Thus, the magnetization $M$ exhibits the properties of a mixed-order
transition.

In region (iv), for $x_{e,M}<x<x_{e}^{\prime}$ we observe a continuous second
order PT between the Para and the Coherent phase occurs at $T_{s,M}$. This is
due to the behavior of the free-energy densities $f(m)$ and $f(M)$ which
display a global minimum at $m=0$ and $M=0$ for $T>T_{s,M}$, while for
$T<T_{s,M}$, a global minimum of $f(M)$ appears at finite $M>0$. The value of
this global minimum of $f(M)$ increases continuously as $T$ is lowered. As $T$
is decreased further and reaches $T_{f}$, new global minima of $f(m_{a})$ and
$f(M)$ appear at certain finite $m_{a}$ and $M$. Thus, a first-order
transition occurs and both order parameters $m_{a}$ and $M$ display a
discontinuous jump. The critical behavior of the model in region (iv) is shown
in Fig. 5(b) and (f).

At the boundary between the regions (iv) and (v), for $x=x_{e}^{\prime}$, we
observe the CE3, where the magnetization $m_{a}$ changes discontinuously from
$0$ to a finite value at ${T_{s}^{\prime}}^{-}$; however, the susceptibility
$\chi_{m}$ diverges at ${T_{s}^{\prime}}^{+}$. Thus, a mixed-order PT occurs
at the CE3.

In region (v), for $x_{e}^{\prime}<x<x_{g}^{\prime}$ as the temperature is
gradually lowered we observe first a send-order PT between the Para and the
Coherent phase at $T_{s,M}$, then we observe another second-order PT between
the Coherent and the Baxter phase at $T_{s}^{\prime}$. In addition to these
two PTs we observe a jump of the order parameter $m_{a}$ and $M$ from non-zero
values. This discontinuity can be obtained by studying the free-energy
densities $f(m_{a})$ and $f(M)$. Indeed when $T>T_{s,M}$, the global minima of
$f(m_{a})$ and $f(M)$ remain at $m_{a}=0$ and $M=0$. For
$T_{s}^{\prime}<T<T_{s,M}$, the global minimum of $f(M)$ occurs at a finite
$M$, which increases continuously as $T$ is lowered gradually. Correspondly,
in this same range of temperatures, the global minimum of $f(m_{a})$ remains
at still $m_{a}=0$. As $T$ gets below $T_{s}^{\prime}$, a global minimum of
$f(m_{a})$ emerges at a finite $m_{a}>0$ in a gradual way. Thus, $m_{a}$ is
finite, and a second-order PT occurs at $T_{s}^{\prime}$. For this same range
of temperature the global minimum of $f(M)$ is achieved at an increasingly
larger value of $M$. When $T$ reaches $T_{f}$, new global minima of $f(m_{a})$
and $f(M)$ emerge at finite $m_{a}$ and $M$, which this minima being separated
from the respective value of the free-energy minima obtained for $T_{f}^{+}$.
Thus, a discontinuity occurs for the order parameter at $T_{f}$. These
behaviors are schematically shown in Fig. 5 (c) and (g).

In region (vi) corresponding to high values of $x$, or $x\to\infty$ two second
order PTs are observed. The first PT occurs between the Para and the Coherent
phase at $T_{s,M}$ and the second PT between the Coherent and the Baxter phase
occurs at $T_{s}^{\prime}$. These behaviors close to these two PTs are
schematically shown in Fig. 5 (d) and (h).

Using Eqs. (100) and (104), we can obtain the following critical behaviors for
$m_{a}$ and $M$:

$\displaystyle m_{a}$ $\displaystyle\sim(T_{s}^{\prime}-T)^{\beta_{m}}$ with
$\displaystyle~{}\beta_{m}$
$\displaystyle=\dfrac{1}{\lambda_{\textrm{min}}-3},$ (45) $\displaystyle M$
$\displaystyle\sim(T_{s,M}-T)^{\beta_{M}}$ with $\displaystyle~{}\beta_{M}$
$\displaystyle=\dfrac{1}{\lambda_{\textrm{o}}-3},$ (46)

where
$\lambda_{\textrm{min}}=\textrm{min}(\lambda_{\textrm{o}},\lambda_{\textrm{n}})$.
Using these results, we obtain the specific heats, which scale as

$\displaystyle C_{m}$ $\displaystyle\sim(T_{s}^{\prime}-T)^{-\alpha_{m}}$ with
$\displaystyle\alpha_{m}$
$\displaystyle=\dfrac{\lambda_{\textrm{min}}-5}{\lambda_{\textrm{min}}-3}\,,$
(47) $\displaystyle C_{M}$ $\displaystyle\sim(T_{s,M}-T)^{-\alpha_{M}}$ with
$\displaystyle\alpha_{M}$
$\displaystyle=\dfrac{\lambda_{\textrm{o}}-5}{\lambda_{\textrm{o}}-3}\,.$ (48)

The susceptibilities behave as follows:

$\displaystyle\chi_{m}\sim\left\\{\begin{array}[]{cccc}(T-T_{s}^{\prime})^{-\gamma_{m}^{+}}&~{}\quad\hbox{with}\quad\gamma_{m}^{+}=1\,,\\\
(T_{s}^{\prime}-T)^{-\gamma_{m}^{-}}&~{}\quad\hbox{with}\quad\gamma_{m}^{-}=1\,,\end{array}\right.$
(51)

and

$\displaystyle\chi_{M}\sim\left\\{\begin{array}[]{cccc}(T-T_{s,M})^{-\gamma_{M}^{+}}&\quad\hbox{with}\quad\gamma_{M}^{+}=1\,,\\\
(T_{s,M}-T)^{-\gamma_{M}^{-}}&\quad\hbox{with}\quad\gamma_{M}^{-}=1\,.\end{array}\right.$
(54)

Detailed derivations of $\chi_{M}$ and $\chi_{m}$ near $T_{s,M}$ and
$T_{s}^{\prime}$, respectively, are given in Appendix E.

Table 1: Critical exponents for $3<(\lambda_{\textrm{n}},\lambda_{\textrm{o}})<4$: Here, $\alpha$ is the exponent of the specific heat, $\beta_{m}$ ($\beta_{M}$) is the exponent of the magnetization $m_{a}$ ($M$) at zero external magnetic field, and $\gamma_{m}$ ($\gamma_{M}$) is the exponent of the susceptibility for $m_{a}$ ($M$)-magnetization near the transition temperature.  Range of $x$ | $\alpha_{m}$ | $\alpha_{M}$ | $\beta_{m}$ | $\beta_{M}$ | $\gamma_{m\pm}$ | $\gamma_{M\pm}$
---|---|---|---|---|---|---
$x=0$ | $\frac{\lambda_{\textrm{mim}}-5}{\lambda_{\textrm{min}}-3}$ | - | $\frac{1}{\lambda_{\textrm{min}}-3}$ | - | $1$ | -
$0<x<x_{M}$ | $\frac{\lambda_{\textrm{min}}-5}{\lambda_{\textrm{min}}-3}$ | $\frac{\lambda_{\textrm{min}}-5}{\lambda_{\textrm{min}}-3}$ | $\frac{1}{\lambda_{\textrm{min}}-3}$ | $\frac{\lambda_{\textrm{o}}-2}{\lambda_{\textrm{o}}-3}$ | $1$ | $0$
$x=x_{M}~{}(\lambda_{\textrm{n}}>\lambda_{\textrm{o}})$ | $\frac{\lambda_{\textrm{o}}-5}{\lambda_{\textrm{o}}-3}$ | $\frac{\lambda_{\textrm{o}}-5}{\lambda_{\textrm{o}}-3}$ | $\frac{1}{\lambda_{\textrm{o}}-3}$ | $\frac{1}{\lambda_{\textrm{o}}-3}$ | $1$ | $1$
$x=x_{M}~{}(\lambda_{\textrm{n}}<\lambda_{\textrm{o}})$ | $\frac{\lambda_{\textrm{n}}-5}{\lambda_{\textrm{n}}-3}$ | $\frac{\lambda_{\textrm{o}}-5}{\lambda_{\textrm{o}}-3}$ | $\frac{1}{\lambda_{\textrm{n}}-3}$ | $\frac{1}{\lambda_{\textrm{o}}-3}$ | $1$ | $1$
$x>x_{M}$ | $\frac{\lambda_{\textrm{min}}-5}{\lambda_{\textrm{min}}-3}$ | $\frac{\lambda_{\textrm{o}}-5}{\lambda_{\textrm{o}}-3}$ | $\frac{1}{\lambda_{\textrm{min}}-3}$ | $\frac{1}{\lambda_{\textrm{o}}-3}$ | $1$ | $1$

###  IV.4 Anomalous Scaling Relations

The critical exponents of the continuous transition are listed in Table 1 for
all ranges of $x$.

The scaling relation for $m_{a}$ satisfies the conventional relation:

$\displaystyle\alpha_{m}+2\beta_{m}+\gamma_{m}=2\,.$ (55)

By contrast, the scaling relation for $M$ shows an unusual behavior for
$x<x_{M}$. The scaling relation for $M$ does not hold for $x<x_{M}$ as

$\displaystyle\alpha_{M}+2\beta_{M}+\gamma_{M}=\left\\{\begin{array}[]{cc}3&~{}\text{for
}\lambda_{\textrm{n}}>\lambda_{\textrm{o}},\\\ \\\
3+2\dfrac{\lambda_{\textrm{o}}-\lambda_{\textrm{n}}}{\lambda_{\textrm{n}}-3}&~{}\text{for
}\lambda_{\textrm{n}}<\lambda_{\textrm{o}}.\end{array}\right.$ (59)

Note that for the original AT model, the scaling relation for $M$ is written
as $\alpha+2\beta_{M}+\gamma_{M}=3$. This relation can be confirmed by setting
$\lambda_{\textrm{n}}=\lambda_{\textrm{o}}$ in the second equation of (59).

##  V Phase diagram II: $\lambda_{a}$-dependence

###  V.1 General remarks

The $g$-AT model may be regarded as a combination of the original AT model on
the network of overlapping links and two independent Ising models on the
respective network of non-overlapping links. In order to fully appreciate the
general phase diagram of the $g$-AT model, let us recall three important
results revealing the interplay between network structure and spin models, the
Ising, Potts, and AT models.

The Ising model on a single SF network with power-law exponent $\lambda$
exhibits a second-order PT at a finite temperature
$T_{c}\propto\langle{k^{2}}\rangle/\langle{k}\rangle$ for $\lambda>3$ within
the annealed approximation Ising2 . Thus, as $\lambda$ increases, $T_{c}$
decreases. Since the magnetization corresponds to the formation of a majority
opinion, this implies that the larger the branching ratio
$\langle{k^{2}}\rangle/\langle{k}\rangle$ of the network is, the easier is for
the network to display a majority opinion. In a single SF network, as the
power-law exponent of degree distributions $\lambda\to 3^{+}$ the branching
ratio of the network increases as
$\langle{k^{2}}\rangle/\langle{k}\rangle\simeq(\lambda-2)/(\lambda-3)$.
Therefore this implies that tuning the power-law exponent $\lambda$, the
network undergoes a topological change that affects the dynamics of spin
model, in particular can modify the value of its critical temperature.

Consequently, we expect that the general phase diagram of the $g$-AT model
will display a significant dependence on the pair of power-law exponents
$(\lambda_{\textrm{n}},\lambda_{\textrm{o}})$. In particular the relative
value of $\lambda_{\textrm{n}}$ with respect to $\lambda_{\textrm{o}}$ allows
to tune the relative influence of non-overlapping multilinks with respect to
overlapping multilinks. We have already seen that $x=J_{4}/J_{2}$ modifies the
phase diagram as it modulates the strength of the $4$-body interactions
(mediated by overlapping multilinks) and the strength of $2$-body interactions
(mediated by non-overlapping multilinks). We expect that the phase diagram
depends on not only $x$ but also the power-law exponents
$(\lambda_{\textrm{n}},\lambda_{\textrm{o}})$ significantly.

Let us recall that the original AT model can be recast in the Potts model with
four states when we set $x=J_{4}/J_{2}=1$ which display a tricritical PT when
the power-law exponent $\lambda_{\textrm{o}}=\lambda_{c}\approx 3.503$. This
implies that for $\lambda_{\textrm{o}}>\lambda_{c}$, the four state Potts
model displays a first-order PT; otherwise, it displays a second-order PT
Potts . From this observation we conclude that the phase diagram of the $g$-AT
model is expected to be more rich around the values
$\lambda_{\textrm{n}}\simeq\lambda_{c}$ and
$\lambda_{\textrm{o}}\simeq\lambda_{c}$.

Figure 6: Schematic phase diagrams of the $g$-AT model in the parameter space
[$\lambda_{\textrm{n}}$, $\lambda_{\textrm{o}}$] for (a) $x\approx 0$, (b)
$x=x_{M}$, and (c) $x\gg x_{M}$. The notations of the phases (i)$-$(vii) are
the same as the ones presented in Figs. 3 and 10.

Figure 7: Schematic phase diagrams of the $g$-AT model in the parameter space
[$x,T^{-1}$] for various values $\lambda_{\textrm{n}}$ and
$\lambda_{\textrm{o}}$. ($\lambda_{\textrm{n}}$, $\lambda_{\textrm{o}}$) = (a)
(3.53, 3.51), (b) (3.53, 3.90), (c) (3.53, 3.49), (d) (3.35, 3.90), (e) (3.53,
3.30), and (f) (3.30, 3.90).

###  V.2 Additional regions of the general phase diagram

In Sec. IV we have described the $x$-dependence of the phase diagram of the
$g$-AT model for a choice of power-law exponents
$(\lambda_{\textrm{n}},\lambda_{\textrm{o}})=(3.53,3.90)$. This phase diagram
displays the PTs of type (i)-(vi), whose implications for opinion dynamics has
been discussed in Sec. IV.

When we consider all possible values of the power-law exponents
$(\lambda_{\textrm{n}},\lambda_{\textrm{o}})$ we observe five more
characteristic regions and lines, denoted as (vii)–(xi) in Figs. 6 (b):

  * –

In region (vii) of Fig. 6 and Fig. 7 (a), two discontinuous behaviors occur
successively as the temperature is lowered. This region can be divided into
four subregions. These subregions are similiar to (ii)-(v) regions in Fig. 7
(b) but a discontinuous jump line additionally exists in a lower temperature
region. The discontinuous jump line originates from the correlations between
$m_{\textrm{n}}$ of one layer and $m_{\textrm{n}}$ or $m_{\textrm{o}}$ of the
other layer, whereas the other discontinuous jump lines in the interval
$[x_{g2},x_{e}]$ and $[x_{e}^{\prime},x_{g1}^{\prime}]$ in a higher
temperature region originates from the correlation between the same
$m_{\textrm{o}}$s but on different layers.

  * 1.

In region (vii)1 ($x_{g2}<x<x_{e}$), a continuous PT between the Para and the
Baxter phase occurs at $T_{s}$ and two discontinuous jumps of the order
parameters $m_{a}$ and $M$ occur successively as the temperature is lowered.
As the noise (temperature) is reduced, the opinion dynamics exhibits first a
continuous PT in which a majority opinion is formed in both layers, and
spreads abruptly twice over a finite fraction of nodes in the multiplex
network.

  * 2.

In region (vii)2 ($x_{e}<x<x_{e,M}$), as the temperature is lowered, a
discontinuous PT between the Para and the Baxter phase occurs firstly and
subsequently a discontinuous jump of the order parameters $m_{a}$ and $M$
occurs in the same Baxter phase. Therefore as the noise is reduced, the
opinion dynamics has first a discontinuous PT in which a majority opinion is
formed in both layers, and then we observe an additional jump in the
magnetization as the majority opinion gets adopted by a larger fraction of
nodes of the multiplex network.

  * 3.

In region (vii)3 ($x_{e,M}<x<x_{e}^{\prime}$), a continuous PT occurs between
the Para and the Coherent phase at $T_{s,M}$. As the temperature is decreased
further, a discontinuous PT occurs at $T_{f}$ from the Coherent phase to the
Baxter phase and a discontinuous jump of the order parameters $m_{a}$ and $M$
occurs at $T_{f}^{\prime}$ successively. This implies that as the noise is
reduced, at temperatures lower than $T_{s,M}$, each single node prefers to
adopt a coherent opinion in both layers. As the temperature is further
reduced, we observe a discontinuous PT in which majority opinion is formed in
both layers, and then the majority opinion spreads abruptly over a finite
fraction of nodes in the multiplex network.

  * 4.

In region (vii)4 ($x_{e}^{\prime}<x<x_{g2}^{\prime}$), two continuous PTs
occur successively: between the Para and the Coherent phase at $T_{s,M}$ and
between the Coherent and the Baxter phase at $T_{s}^{\prime}$. Then as the
temperature is decreased further two discontinuous jumps of the order
parameter occur at $T_{f}$ and $T_{f}^{\prime}$, respectively. Overall, as the
noise is reduced, the opinion dynamics exhibits two types of continuous PTs
successively in which the coherent and the majority opinion are formed,
respectively. When the noise is decreased further, the majority opinion
spreads abruptly twice over a finite fraction of nodes in the multiplex
network.

  * –

In regions or points (viii) and (x) of Figs. 6 and Fig. 7 (c) and (d),
respectively, a continuous PT between the Para and the Baxter phase and a
discontinuous jump of the order parameters $m_{a}$ and $M$ occur successively
as the temperature is lowered. The critical behavior at $x_{M}$ differs from
that in region (ii) as we observe at a tricritical point. The transition point
at $x_{M}$ in (viii) and (x) acts as a branching point of the critical line to
the two critical lines of the Coherent phase. Overall, as the noise is reduced
the opinion dynamics exhibits first a continuous PT in which a majority
opinion is formed in both layers, and subsequently the majority opinion
spreads abruptly over a large fraction of nodes in the multiplex network.

  * –

In regions (ix) and (xi) of Fig. 6 and in Fig. 7 (e) and (f) at $x_{M}$,
respectively, a continuous PT between the Para and the Baxter phase is
observed as the temperature is lowered. The critical behavior at $x_{M}$
differs from that in (i) as we observe at a tricritical point. The transition
point at $x_{M}$ in (ix) and (xi) acts as a branching point of the critical
line to the two critical lines of the Coherent phase. As the noise is reduced,
the opinion dynamics has continuous PT in which a majority opinion is formed
in both layers.

In order to explore the dependence of this rich phase diagram on the exponents
$\lambda_{\textrm{n}}$ and $\lambda_{\textrm{o}}$ around $x=x_{M}$, we plot
the phase diagram in the space [$x,T^{-1}$] for various values of the degree
pairs $(\lambda_{\textrm{n}},\lambda_{\textrm{o}})$ (see Fig. 7). We find that
if a first-order PT occurs at $x=x_{M}$, the overall phase diagram is close to
the phase diagram discussed in Sec. IV (see Fig. 3 and Fig. 7 (b)). If $x_{M}$
demarks the boundary between type-(ii) and type-(v) PTs, then the phase
diagram is similar to Fig. 7 (c) and (d), respectively. When
$\lambda_{\textrm{n}}<\lambda_{\textrm{o}}\approx 3.90$ (see Fig. 7 (b), (d),
and (f)), as $\lambda_{\textrm{n}}$ is decreased, a discontinuous transition
curve shrinks and moves left and upward as shown in Figs. 7 (b)$\to$(f).
Moreover, when $\lambda_{\textrm{n}}$ is slightly larger than
$\lambda_{\textrm{o}}\approx\lambda_{c}^{+}$ (see Fig. 7 (a)), double
discontinuous transition curves appear in the phase diagram, where two
discontinuous PTs occur successively as $T$ is decreased.

###  V.3 Free-energy landscape for the $\lambda$-dependence of phase
transitions at $x\approx x_{M}$

Here, we will investigate the free-energy landscape of the $g$-AT model for
the (vii)–(xi) types of PTs. Around $x\approx x_{M}$, phase and PT type are
determined by the free energy density presented in Appendix D.

In region (vii)1, PT type is determined by the free energy density given as
Eq. (87) because $x<x_{M}$. Note that the higher-order term (89) of Eq. (87)
is negative. As $x\to x_{M}$, the terms with $B_{\textrm{o}}$ of Eq. (89) and
with $D_{m}$ of Eq. (87) become comparable in their magnitudes to the terms
with $C_{3}$ and $D_{\textrm{o}}$ in Eq. (106), respectively. Thus, these
terms with $D_{m}$ and $B_{\textrm{o}}$ play a similar role to the terms with
$C_{\textrm{o}}$ and $C_{3}$. For $T>T_{s}$, the terms with $D_{m}$ and
$B_{\textrm{o}}$ are not large in magnitude, so that the global minimum of
$f(m_{a})$ remains at $m_{a}=0$ and $M=0$, and thus a continuous PT occurs at
$T_{s}$. However, when $T$ is lowered further, the term with $D_{m}$ increases
and becomes comparable to the leading order terms. Then a discontinuous jump
of the order parameters occurs at $T_{f}$. As $T$ is lowered further, another
negative term with $B_{\textrm{o}}$ term increases, another jump of the order
parameters occurs at $T_{f}^{\prime}$. Hence in the region (vii)1, as $T$ is
lowered from $T_{s}^{+}$, a second-order PT occurs firstly and then two
discontinuous jumps occur successively in the Baxter phase.

In region (vii)2, the term with $C_{\textrm{o}}$, induced by the correlation
between $m_{\textrm{o}}$s on different layers, becomes negative when
$\lambda_{\textrm{o}}>\lambda_{c}$, thus produces a discontinuous PT at
$T_{f}$ higher than $T_{s}$. Another negative term with $C_{3}$, induced by
the correlations between $m_{\textrm{n}}$ of one layer and $m_{\textrm{n}}$ or
$m_{\textrm{o}}$ of the other layer, becomes larger as $T$ is decreased, and
thus a discontinuous jump of the order parameters $m_{a}$ and $M$ in the
Baxter phase occurs at $T_{f}^{\prime}<T_{s}$. $f(m_{a})$ and $f(M)$ develop a
global minimum at a temperature $T_{f}$, leading to a discontinuous PT between
the Para and the Baxter phase. As $T$ is further lowered from $T_{f}$, the
global minimum position of $m_{a}$ and $M$ increases continuously until a
certain temperature $T_{f}^{\prime}$. When $T$ reaches $T_{f}^{\prime}$,
another global minimum of $f(m_{a})$ and $f(M)$ emerge at another finite
$m_{a}$ and $M$, which lead to the jumps in the order parameters $m_{a}$ and
$M$ in the Baxter phase. The order parameters and free energy landscape in
this region are depicted in Fig. 8 (a) and (b), respectively.

In region (vii)3, PT type is determined by the free energy density Eq. (100),
because $x>x_{M}$. Note that the higher-order term (102) of Eq. (100) is
negative. As $x\to x_{M}$, the terms with $B_{\textrm{o}}^{\prime}$ of Eq.
(102) and with $D_{M}$ of Eq. (100) become comparable in their magnitudes to
the terms with $C_{3}$ and $D_{\textrm{o}}$ in Eq. (106), respectively. Thus,
these terms with $D_{M}$ and $B_{\textrm{o}}^{\prime}$ play a similar role to
the terms with $C_{\textrm{o}}$ and $C_{3}$. For $T>T_{s,M}$, the terms with
$D_{M}$ and $B_{\textrm{o}}^{\prime}$ are small in magnitude, so that a global
minimum of $f(m_{a})$ and $f(M)$ remain at $m_{a}=0$ and $M=0$, respectively.
As $T$ is lowered from $T_{s,M}^{+}$, a second-order PT for $M$ occurs from
the Para to the Coherent phase at $T_{s,M}$ and a global minimum of $f(M)$
increases continuously. When $T$ is decreased further to $T_{f}$, the term
with $D_{M}$ increases in its magnitude, a new global minimum of $f(m_{a})$
and $f(M)$ appears far from values at $T_{f}^{+}$, respectively, a
discontinuous PT between the Coherent and the Baxter phase appears at $T_{f}$.
When $T$ is lowered further, the term with $B_{\textrm{o}}^{\prime}$ becomes
large, a discontinuous jump of the order parameters occurs at $T_{f}^{\prime}$
in the Baxter phase. Hence, in the region (vii)3, as $T$ is lowered from
$T_{s,M}^{+}$, a continuous PT from the Para to the Coherent phase occurs
firstly at $T_{s,M}$ and then a discontinuous PT between the Coherent and the
Baxter phase occurs at $T_{f}$ and then a discontinuous jump occurs at
$T_{f}^{\prime}$ successively.

In region (vii)4, PT type is investigated through Eq. (100), because
$x>x_{M}$. When $T>T_{s,M}$, the terms with $D_{M}$ and
$B_{\textrm{o}}^{\prime}$ are too small, and the global minimum of $f(m_{a})$
and $f(M)$ remain at $m_{a}=0$ and $M=0$, respectively. For
$T_{s}^{\prime}<T<T_{s,M}$, a second-order PT for $M$ from the Para to the
Coherent phase occurs at $T_{s,M}$, and the global minimum of $f(M)$ grows
continuously as $T$ is lowered from $T_{s,M}^{+}$. Meanwhile, the global
minimum of $f(m_{a})$ remains at still $m_{a}=0$. As $T$ is lowered across
$T_{s}^{\prime}$, a second-order PT for $m_{a}$ from the Coherent to the
Baxter phase occurs at $T_{s}^{\prime}$. When $T$ is lowered and reaches
$T_{f}$ and $T_{f}^{\prime}(<T_{f})$, the global minimum of $f(m_{a})$ and
$f(M)$ jump discontinuously from the previous positions at $T_{f}^{+}$ and
$T_{f}^{{}^{\prime}+}$, respectively. Hence, in the region (vii)4, as $T$ is
decreased from $T_{s,M}^{+}$, a continuous PT from the Para to the Coherent
phase occurs at $T_{s,M}$ and a continuous PT from the Coherent to the Baxter
phase occurs at $T_{s}^{\prime}$, and then two discontinuous jump of the order
parameters in the Baxter phase occur at $T_{f}$ and $T_{f}^{\prime}$,
successively.

In region (viii), $f(m_{a})$ and $f(M)$ still display a global minimum at
$m=0$ and $M=0$ for $T>T_{s}$ and a continuous PT between the Para and the
Baxter phase at $T=T_{s}$. The values of the critical exponents for this PT
are listed in Table I for the case $x\approx x_{M}$ and
$\lambda_{\textrm{n}}>\lambda_{\textrm{o}}$. However, as $T$ is further
lowered below a certain temperature $T_{f}$, a global minimum emerges at a
non-zero value of the magnetization $m_{a}>0$ and $M>0$, and a discontinuous
transition occurs. Thus, as $T$ is decreased from $T_{s}^{+}$, a continuous
transition occurs firstly at $T_{s}$, followed by a discontinuous transition
at $T_{f}$. The order parameters and free energy density landscape are
depicted in Fig. 9 (a) and (c), respectively.

In region (ix), $f(m_{a})$ and $f(M)$ have a global minimum at $m_{a}=0$ and
$M=0$, respectively, for $T\geq T_{s}$ while for $T=T_{s}^{-}$ a global
minimum emerges continuously at a non-zero value of the magnetization.
Therefore at $T=T_{s}$ we observe a continuous PT between the Para and the
Baxter phase. The values of the critical exponents for this PT are listed in
Table I for the case $x\approx x_{M}$ and
$\lambda_{\textrm{n}}>\lambda_{\textrm{o}}$. The order parameters and free
energy density landscape are illustrated in Fig. 9 (b) and (d), respectively.
This continuous transitions is similar to the (i)-type PTs; however, the
critical behavior for $M$ of this continuous transition differs from that of
the (i)-type PT, and thus we denote this type of a continuous transition as
the (ix)-type PT to distinguish this from the (i)-type PTs.

Figure 8: For the (vii)-type of PT, schematic plots of (a) the order
parameters $m_{a}$ and $M$ as a function of $T$ and (b) the free energy
density landscape as a function of $m_{\textrm{o}}$ for various $T$s. The
exponents of degree distributions are taken as ($\lambda_{\textrm{n}}$,
$\lambda_{\textrm{o}}$) = (3.53, 3.51).

Figure 9: For the (viii)- and (ix)-type of PTs, schematic plots of (a) and (b)
the order parameters $m_{a}$ and $M$ as a function of $T$, respectively. (c)
and (d) schematic plots of the free energy density landscape as a function of
$m_{\textrm{o}}$ for the temperatures around $T_{f}$. The exponents of degree
distributions for (a) and (c) are taken as ($\lambda_{\textrm{n}}$,
$\lambda_{\textrm{o}}$) = (3.53, 3.30) and for (b) and (d) are taken as (3.53,
3.49).

Figure 10: (a) and (b) Plots of the order parameters $m_{a}$ and $M$ at
$x=x_{M}$ as a function of $T/J_{2}$. (c) and (d) Plots of the free energy
density landscape as a function of $m_{\textrm{n}}$ for various $T$. The
exponents of degree distributions taken for (a) and (c) are
($\lambda_{\textrm{n}}$, $\lambda_{\textrm{o}}$) = (3.30, 3.55); for (b) and
(d) are (3.43, 3.55).

In region (x), $f(m_{a})$ and $f(M)$ have a global minimum at $m_{a}=0$ and
$M=0$, respectively. $f(m_{a})$ remains at $m_{a}=M=0$ for $T>T_{s}$. When $T$
is lower than $T_{s}$, a continuous transition occurs. The values of the
critical exponents for this PT are listed in Table I for the case $x\approx
x_{M}$ and $\lambda_{\textrm{n}}<\lambda_{\textrm{o}}$. As $T$ is further
lowered, the order parameter gradually increases. When $T$ reaches $T_{f}$,
the order parameter jumps by a finite amount and a new global minimum of
$f(m_{a})$ occurs at a finite $m_{a}$. Thus, as $T$ is decreased from
$T_{s}^{+}$, a continuous PT occurs at $T_{s}$ firstly and then a
discontinuous jump of the order parameter occurs at $T_{f}$ as shown in Fig.
10(a) and (b).

In region (xi), $f(m_{a})$ and $f(M)$ have a global minimum at $m_{a}=0$ and
$M=0$, respectively. They remain at $m_{a}=M=0$ for $T>T_{s}$. At $T_{s}$,
$m_{a}$ and $M$ exhibit continuous PTs. The values of the critical exponents
for this PT are listed in Table I for the case $x\approx x_{M}$ and
$\lambda_{\textrm{n}}<\lambda_{\textrm{o}}$. When $T$ is decreased from
$T_{s}$, a global minimum occurs at finite $m_{a}$ and $M$. These behaviors
are schematically shown in Fig. 10 (b) and (d). Note that the $\beta_{M}$ and
$\chi_{M}$ of the (xi)-type PTs are different from those of the (i)-type PTs
(Eq. (39) and $\gamma_{M}=0$).

Figure 11: Schematic contour lines of $x_{e,M}$ in the parameter space
[$\lambda_{\textrm{n}}$, $\lambda_{\textrm{o}}$]. The rightmost (light yellow)
region contains the domains (viii), (ix), (x) and (xi) as shown in Fig. 6 (b).
In this region, a continuous PT appears at $T_{s}$ at $x_{M}$, and thus
$x_{e,M}=x_{M}$, the Coherent phase appears at $x_{M}$. The regions denoted as
\scriptsize1⃝$-$\scriptsize4⃝ correspond to the regions (iii) and (vii)
denoted in Fig. 6 (b). In this region, a discontinuous PT appears at $x_{M}$,
and thus $x_{e,M}>x_{M}$. The contour lines represent in term of the ratio
$x_{e,M}/x_{M}$. In \scriptsize1⃝, $1<x_{e,M}/x_{M}<1.05$; in \scriptsize2⃝,
$1.05<x_{e,M}/x_{M}<1.10$; in \scriptsize3⃝, $1.10<x_{e,M}/x_{M}<1.15$; and in
\scriptsize4⃝, $1.15<x_{e,M}/x_{M}$.

###  V.4 $\lambda$-dependence of the Coherent phase

When $x\to 0$, a second-order PT occurs from the Baxter to the Para phase,
whereas when $x\gg x_{M}$, a first-order PT occurs from the Baxter to the
Coherent phase (i.e., $\langle\sigma s\rangle$), followed by another PT occurs
from the Coherent to the Para phase. The Coherent phase appears for $x\geq
x_{e,M}(\lambda_{a})$. In Fig. 11, we display the $\lambda_{a}$-dependence of
$x_{e,M}$ in unit of $x_{M}$. If $x_{e,M}=x_{M}^{+}$, then the Coherent phase
appears in the range of $x>x_{M}$ as shown in Fig. 7 (a), (b), (d), and (e),
which is denoted as (ix), (viii), (xi), and (x) (light yellow) in Fig. 11,
respectively. In the regions (iii) and (vii), when a discontinuous PT occurs
at $x_{M}$, then $x_{e,M}>x_{M}$, and the Coherent phase appears in range of
$x>x_{e,M}$ as shown in Fig. 7 (c) and (f), respectively. The contour lines
between different regions with different circle numbers represent different
ratios $x_{e,M}/x_{M}$.

We find that the transition point $x_{e,M}$ is delayed as both
$\lambda_{\textrm{n}}$ and $\lambda_{\textrm{o}}$ are increased. This is
caused by the following reasons: if $\lambda_{\textrm{n}}$ and
$\lambda_{\textrm{o}}$ are large, then the branching ratios of non-ovelapping
and overlapping links, respectively, become small. Thus, a larger value of
$x=J_{4}/J_{2}$ is needed, i.e., the strength of 4-body interaction needs to
be reinforced to form a Coherence phase.

##  VI Conclusion

To investigate the effect of link overlap on the opinion dynamics defined on a
multiplex network, we studied the so-called $g$-AT model, a spin model in
thermal equilibrium systems. The $g$-AT model describes the dynamics of two
species of Ising spins, namely the $s$ and $\sigma$ spins, each of which is
located on a single layer of the duplex network under consideration. Here, the
spin model is defined on duplex networks with an SF multidegree distribution,
which facilitates tuning of the effect of overlapping links with respect to
non-overlapping links. In particular, we distinguish between multilinks
$(1,1)$ that characterize overlapping links and multilinks $(1,0)$ and $(0,1)$
that do not. We assume that the multidegrees $k^{(1,1)}=k_{\textrm{o}}$ and
$k^{(1,0)}=k^{(0,1)}=k_{\textrm{n}}$ follow power-law distributions associated
with the tunable power-law exponents $\lambda_{\textrm{n}}$ and
$\lambda_{\textrm{o}}$. This system is illustrated in Fig. 1.

Pairs of $s$-spins (pairs of $\sigma$-spins) connected by overlapping and non-
overlapping links interact through a 2-body interaction of strength $J_{2}$.
Four spins comprising two $s$-spins and two $\sigma$-spins connected by
overlapping links interact through a 4-body interaction with strength $J_{4}$
(see Fig. 2). The ratio $x\equiv J_{4}/J_{2}$ is a control parameter that can
alter the critical properties of the model, and the system is assumed to be in
thermal contact with a heat reservoir at temperature $T$. Here, $T$ represents
the diversity of each individual opinion in a social community. Thus, there
exist four control parameters, namely $\lambda_{\textrm{n}}$,
$\lambda_{\textrm{o}}$, $x$, and $T$. By applying the Landau–Ginzburg theory,
we obtained rich phase diagrams in the four-parameter space. The $g$-AT model
is a generalization of the original AT model AT , in which all the links are
regarded as overlapping links; therefore, a single exponent
$\lambda_{\textrm{o}}$ is considered in this context.

We note that the different species of spins represent individuals from two
different communities formed based on friendship and business relations,
respectively. Each pair of individuals may be connected solely via friendship
links, solely via business relations, or via both relationships. The formation
of a majority opinion across both layer is indicated by the magnetizations
$\langle\sigma\rangle>0$, $\langle s\rangle>0$, and $\langle\sigma s\rangle>0$
in the spin model, which can be accomplished through non-overlapping or
overlapping links. The diversity of individual opinions is reflected by
thermal fluctuations.

We investigated PTs arising from the competition between the consensus
formation of each community and that of the entire society, and obtained rich
phase diagrams including diverse types of PTs. These findings are expected to
be beneficial in understanding the underlying mechanisms of local and global
formation of a majority opinion in a society.

Similar to the voter models on multiplex networks Chmiel1 , the $g$-AT model
shows that a majority opinion emerges abruptly thanks to the interactions
across two layers induced by overlapping links. In particular, in the $g$-AT
model, we can control the strength of 4-body interactions among replica nodes
connected by overlapping multilinks with respect of the strenght of 2-body
interactions $x=J_{4}/J_{2}$. This allows us to assess the role of tuning the
strength of the many-body AT-interactions by modulating $x$ and study how the
phase digram change with respect to the original AT model AT . Moreover we can
tune the power-law exponent of the overlapping multidegree distributions of
non-overlapping and overlapping multilinks
$(\lambda_{\textrm{n}},\lambda_{\textrm{o}})$, and investigate the role of
this topological modifications on PTs, in the same spirit of the analysis
conducted for percolation problems in Ref. Cellai1 ; Cellai2 . In the future
our work can be expanded in many directions, investigating further the role
that higher-order interactions battiston ; HO  have in opinion dynamics
defined on multiplex networks and exploring realistic spin models of opinion
dynamics defined on duplex networks.

##  Appendix A Self-consistency equation with an external magnetic field

With the external magnetic field

$\displaystyle\Omega_{a}\equiv
K_{2}m_{a}+H_{a},\quad\textrm{and}\quad\Omega_{M}\equiv K_{4}M+H_{4},$ (60)

the self-consistency equations Eqs. (32) for $m_{a}$ and $M$ of the $g$-AT
model are replaced as follows:

$\displaystyle m_{a}\langle k_{a}\rangle$
$\displaystyle=m_{a\textrm{I}}\equiv\int_{k_{\rm
min}^{\textrm{n}}}^{\infty}\int_{k_{\rm
min}^{\textrm{o}}}^{\infty}dk_{\textrm{n}}dk_{\textrm{o}}P_{d}(k_{\textrm{n}})P_{d}(k_{\textrm{o}})\,\dfrac{\tanh\left(\Omega_{\textrm{o}}k_{\textrm{o}}+\Omega_{\textrm{n}}k_{\textrm{n}}\right)\left[1+\tanh\left(\Omega_{4}k_{\textrm{o}}\right)\right]}{1+\tanh^{2}\left(\Omega_{\textrm{o}}k_{\textrm{o}}+\Omega_{\textrm{n}}k_{\textrm{n}}\right)\tanh\left(\Omega_{4}k_{\textrm{o}}\right)}k_{a}$
(61)

and

$\displaystyle M\langle k_{\textrm{o}}\rangle$
$\displaystyle=M_{\textrm{I}}\equiv\int_{k_{\rm
min}^{\textrm{n}}}^{\infty}\int_{k_{\rm
min}^{\textrm{o}}}^{\infty}dk_{\textrm{n}}dk_{\textrm{o}}P_{d}(k_{\textrm{n}})P_{d}(k_{\textrm{o}})\dfrac{\tanh\left(\Omega_{4}k_{\textrm{o}}\right)+\tanh^{2}\left(\Omega_{\textrm{o}}k_{\textrm{o}}+\Omega_{\textrm{n}}k_{\textrm{n}}\right)}{1+\tanh^{2}\left(\Omega_{\textrm{o}}k_{\textrm{o}}+\Omega_{\textrm{n}}k_{\textrm{n}}\right)\tanh\left(\Omega_{4}k_{\textrm{o}}\right)}k_{\textrm{o}}\,.$
(62)

##  Appendix B Definitions of the $\mathcal{A}$ terms in exact susceptibility
formula

The $\mathcal{A}$ terms are defined as follows:

$\displaystyle\mathcal{A}_{aa}$ $\displaystyle=\int_{k_{\rm
min}^{\textrm{n}}}^{\infty}\int_{k_{\rm
min}^{\textrm{o}}}^{\infty}dk_{\textrm{n}}dk_{\textrm{o}}P_{d}(k_{\textrm{n}})P_{d}(k_{\textrm{o}})\dfrac{(1-\mathcal{T}_{2}^{2}\mathcal{T}_{4})(1+\mathcal{T}_{4})}{(1+\mathcal{T}_{2}^{2}\mathcal{T}_{4})^{2}\cosh^{2}\left(K_{2}(m_{\textrm{o}}k_{\textrm{o}}+m_{\textrm{n}}k_{\textrm{n}})\right)}k_{a}^{2}\,,$
(63) $\displaystyle\mathcal{A}_{a\bar{a}}=\mathcal{A}_{\bar{a}a}$
$\displaystyle=\int_{k_{\rm min}^{\textrm{n}}}^{\infty}\int_{k_{\rm
min}^{\textrm{o}}}^{\infty}dk_{\textrm{n}}dk_{\textrm{o}}P_{d}(k_{\textrm{n}})P_{d}(k_{\textrm{o}})\dfrac{(1-\mathcal{T}_{2}^{2}\mathcal{T}_{4})(1+\mathcal{T}_{4})}{(1+\mathcal{T}_{2}^{2}\mathcal{T}_{4})^{2}\cosh^{2}\left(K_{2}(m_{\textrm{o}}k_{\textrm{o}}+m_{\textrm{n}}k_{\textrm{n}})\right)}k_{\textrm{n}}k_{\textrm{o}}\,,$
(65) $\displaystyle\mathcal{A}_{aM}$ $\displaystyle=\int_{k_{\rm
min}^{\textrm{n}}}^{\infty}\int_{k_{\rm
min}^{\textrm{o}}}^{\infty}dk_{\textrm{n}}dk_{\textrm{o}}P_{d}(k_{\textrm{n}})P_{d}(k_{\textrm{o}})\dfrac{(1-\mathcal{T}_{2}^{2})\mathcal{T}_{4}}{(1+\mathcal{T}_{2}^{2}\mathcal{T}_{4})^{2}\cosh^{2}(K_{4}Mk_{\textrm{o}})}k_{a}k_{\textrm{o}}\,,$
(67) $\displaystyle\mathcal{A}_{Ma}$ $\displaystyle=\int_{k_{\rm
min}^{\textrm{n}}}^{\infty}\int_{k_{\rm
min}^{\textrm{o}}}^{\infty}dk_{\textrm{n}}dk_{\textrm{o}}P_{d}(k_{\textrm{n}})P_{d}(k_{\textrm{o}})\dfrac{2\mathcal{T}_{2}(1-\mathcal{T}_{4}^{2})}{(1+\mathcal{T}_{2}^{2}\mathcal{T}_{4})^{2}\cosh^{2}\left(K_{2}(m_{\textrm{o}}k_{\textrm{o}}+m_{\textrm{n}}k_{\textrm{n}})\right)}k_{\textrm{o}}k_{a}\,,$
(69) $\displaystyle\mathcal{A}_{MM}$ $\displaystyle=\int_{k_{\rm
min}^{\textrm{n}}}^{\infty}\int_{k_{\rm
min}^{\textrm{o}}}^{\infty}dk_{\textrm{n}}dk_{\textrm{o}}P_{d}(k_{\textrm{n}})P_{d}(k_{\textrm{o}})\dfrac{1-\mathcal{T}_{2}^{4}}{(1+\mathcal{T}_{2}^{2}\mathcal{T}_{4})^{2}\cosh^{2}(K_{4}Mk_{\textrm{o}})}k_{\textrm{o}}^{2}\,.$
(71)

##  Appendix C Definitions of coefficients in the free energy density

The coefficients $C_{\textrm{o}}$, $C_{\textrm{n}}$, $D_{m}$, $D_{M}$,
$D_{0}$, and $C_{\textrm{o}}(r_{0},\lambda_{\textrm{o}})$ used in the Landau
free energy formulas in Appendix D are defined as follows:

$\displaystyle C_{\textrm{o}}(\lambda_{\textrm{o}})$
$\displaystyle=C_{M}(\lambda_{\textrm{o}})=-N_{\textrm{o}}\int_{0}^{\infty}\left[\ln(\cosh
y)-\frac{1}{2}y^{2}\right]y^{-\lambda_{\textrm{o}}}dy,$ (72) $\displaystyle
C_{\textrm{n}}(\lambda_{\textrm{n}})$
$\displaystyle=-N_{\textrm{n}}\int_{0}^{\infty}\left[\ln(\cosh
y)-\frac{1}{2}y^{2}\right]y^{-\lambda_{\textrm{n}}}dy,$ (73) $\displaystyle
D_{m}(\lambda_{\textrm{o}})$
$\displaystyle=-N_{\textrm{o}}\int_{0}^{\infty}y\ln\left(1+\tanh^{2}y\right)y^{-\lambda_{\textrm{o}}}dy,$
(74) $\displaystyle D_{M}(\lambda_{\textrm{o}})$
$\displaystyle=-N_{\textrm{o}}\int_{0}^{\infty}y^{2}\ln\left(1+\tanh
y\right)y^{-\lambda_{\textrm{o}}}dy,$ (75) $\displaystyle
D_{0}(r_{0},\lambda_{\textrm{o}})$
$\displaystyle=-N_{\textrm{o}}\int_{0}^{\infty}\ln\left(1+\tanh^{2}y\tanh(r_{0}y)\right)y^{-\lambda_{\textrm{o}}}dy,$
(76) $\displaystyle C_{\textrm{o}}(r_{0},\lambda_{\textrm{o}})$
$\displaystyle=C_{\textrm{o}}(\lambda_{\textrm{o}})+D_{0}(r_{0},\lambda_{\textrm{o}})+\frac{1}{2}C_{M}(\lambda_{\textrm{o}}),$
(77)

where $N_{\textrm{o}}$ and $N_{\textrm{n}}$ are normalization factors written
as ${1}/({\lambda_{\textrm{o}}-1})$ and ${1}/({\lambda_{\textrm{n}}-1})$,
respectively.

##  Appendix D Landau Free Energy Formula

To investigate the critical behavior near the critical temperature, we expand
the free energy density as a function of the order parameters $m_{a}$ and then
analyze the leading terms when $m_{a}$ and $M$ converge to $0$. To proceed, it
is necessary to derive the relation between $m_{a}$ and $M$, which turns out
to depend on the ratio $x$. For values of $x$ smaller and bigger than the
characteristic ratio $x_{M}$, we observe different behaviors. Here we discuss
in details the cases $0<x<x_{M}$, $x>x_{M}$, and $x=x_{M}$, separately.

###  D.1 Case $0<x<x_{M}$

For $x\in(0,x_{M})$ we can expand Eqs. (29) and (30) in terms of $m_{a}$ and
$M$ within the lowest order terms as follows:

$\displaystyle m_{\textrm{n}}\langle
k_{\textrm{n}}\rangle\left(1-K_{2}\dfrac{\langle
k_{\textrm{n}}^{2}\rangle}{\langle k_{\textrm{n}}\rangle}\right)$
$\displaystyle\simeq K_{2}m_{\textrm{o}}\langle k_{\textrm{n}}\rangle\langle
k_{\textrm{o}}\rangle-(\lambda_{\textrm{n}}-1)C_{\textrm{n}}(\lambda_{\textrm{n}})(K_{2}m_{\textrm{n}})^{\lambda_{\textrm{n}}-2}+\hbox{h.o.,}$
(78) $\displaystyle m_{\textrm{o}}\langle
k_{\textrm{o}}\rangle\left(1-K_{2}\dfrac{\langle
k_{\textrm{o}}^{2}\rangle}{\langle k_{\textrm{o}}\rangle}\right)$
$\displaystyle\simeq K_{2}m_{\textrm{n}}\langle k_{\textrm{n}}\rangle\langle
k_{\textrm{o}}\rangle-(\lambda_{\textrm{o}}-1)C_{\textrm{o}}(\lambda_{\textrm{o}})(K_{2}m_{\textrm{o}})^{\lambda_{\textrm{o}}-2}+(\lambda_{\textrm{o}}-2)D_{m}(\lambda_{\textrm{o}})(K_{4}M)(K_{2}m_{\textrm{o}})^{\lambda_{\textrm{o}}-3}+\hbox{h.o.,}$
(80) $\displaystyle M\langle k_{\textrm{o}}\rangle\left(1-xK_{2}\dfrac{\langle
k_{\textrm{o}}^{2}\rangle}{\langle k_{\textrm{o}}\rangle}\right)$
$\displaystyle\simeq
D_{m}(\lambda_{\textrm{o}})(K_{2}m_{\textrm{o}})^{\lambda_{\textrm{o}}-2}-(\lambda_{\textrm{o}}-1)C_{M}(\lambda_{\textrm{o}})(K_{4}M)^{\lambda_{\textrm{o}}-2}+\hbox{h.o.,}$
(82)

where the coefficients of the entropy terms $C_{a}(\lambda_{a})$ and
$C_{M}(\lambda_{\textrm{o}})$ are presented in Appendix C and the coefficient
of the interlayer interaction term $D_{m}(\lambda_{\textrm{o}})$ of the r.h.s.
of Eq. 80 is also presented in Appendix C. This term needs to be considered,
because it is negative and contributes to the first-order transition.

To obtain $T_{s}$, we first consider the lowest-order terms of Eqs. (78) and
(80) and obtain the following:

$\displaystyle\left(1-\dfrac{\langle k_{\textrm{o}}^{2}\rangle}{\langle
k_{\textrm{o}}\rangle}\dfrac{1}{T}\right)\left(1-\dfrac{\langle
k_{\textrm{n}}^{2}\rangle}{\langle
k_{\textrm{n}}\rangle}\dfrac{1}{T}\right)-\dfrac{1}{T^{2}}\langle
k_{\textrm{n}}\rangle\langle k_{\textrm{o}}\rangle=0\,.$ (83)

This equation has two solutions of $T$, denoted as $T_{\ell}$ and $T_{h}$.

It is guaranteed that the l.h.s of Eq. (82) is positive as $T\to T_{h}^{-}$ as
long as $x<x_{M}$ with

$\displaystyle x_{M}\equiv T_{h}{\langle k_{\textrm{o}}\rangle}/{\langle
k_{\textrm{o}}^{2}\rangle}.$ (84)

Then $M$ is written within a leading order as

$\displaystyle M\simeq\frac{D_{m}(\lambda_{\textrm{o}})}{\langle
k_{\textrm{o}}\rangle[1-xT_{h}/(x_{M}T)]}(K_{2}m_{\textrm{o}})^{\lambda_{\textrm{o}}-2}+\hbox{h.o.}\,.$
(85)

This implies that $O(M)\ll O(m_{a})$ near $T_{h}^{-}$ for
$\lambda_{\textrm{o}}>3$. Using this relation (85), we obtain the self-
consistency relations for $m_{a}$ with leading terms as follows:

$\displaystyle m_{a}\langle k_{\textrm{n}}\rangle\langle
k_{\textrm{o}}\rangle\left(1-\frac{T_{\ell}}{T}\right)\left(1-\frac{T_{h}}{T}\right)\simeq-(\lambda_{a}-1)E_{\bar{a}}C_{a}(\lambda_{a})(K_{2}m_{a})^{\lambda_{a}-2}-(\lambda_{\bar{a}}-1)F_{a\bar{a}}C_{\bar{a}}(\lambda_{\bar{a}})(K_{2}m_{\bar{a}})^{\lambda_{\bar{a}}-2}\,,$
(86)

where $E_{\bar{a}}=\langle k_{\bar{a}}\rangle\left(1-\dfrac{\langle
k_{\bar{a}}^{2}\rangle/\langle k_{\bar{a}}\rangle}{T}\right)$ and
$F_{a\bar{a}}=\dfrac{1}{T}\langle k_{\textrm{n}}\rangle\langle
k_{\textrm{o}}\rangle$.

As $T\to T_{h}^{-}$, $m_{a}\to 0$ and thus $M\to 0$ in Eq. (85). This implies
that the PT from Baxter to Para phase is continuous. We confirm that a second-
order transition occurs at $T_{h}$. This temperature is denoted as a critical
temperature $T_{s}\equiv T_{h}$.

Using Eq. (85), we expand the free energy density of Eq. (26) with respect to
$m_{a}$ up to the three lowest order terms:

  * (i)

For $\lambda_{\textrm{n}}>\lambda_{\textrm{o}}$,

$\displaystyle f(m_{\textrm{o}})\simeq
A_{\textrm{o}}K_{2}m_{\textrm{o}}^{2}\left(1-\frac{T_{s}}{T}\right)+C_{\textrm{o}}(\lambda_{\textrm{o}})(K_{2}m_{\textrm{o}})^{\lambda_{\textrm{o}}-1}+C_{\textrm{n}}(\lambda_{\textrm{n}})(B_{\textrm{o}}m_{\textrm{o}})^{\lambda_{\textrm{n}}-1}-\frac{1}{2}\frac{K_{4}[D_{m}(\lambda_{\textrm{o}})]^{2}}{\langle
k_{\textrm{o}}\rangle[1-xT_{s}/(x_{M}T)]}(K_{2}m_{\textrm{o}})^{2(\lambda_{\textrm{o}}-2)}+\hbox{h.o.},$
(87)

where $A_{\textrm{o}}$ and $B_{\textrm{o}}$ are functions of
$\lambda_{\textrm{a}}$ and $K_{2}$, for which the explicit formula is as
follows:

$\displaystyle A_{\textrm{o}}(\lambda_{\textrm{o}})=\dfrac{K_{2}\langle
k_{\textrm{n}}\rangle^{2}K_{2}\langle k_{\textrm{o}}\rangle^{2}}{\langle
k_{\textrm{n}}\rangle-K_{2}\langle k_{\textrm{n}}^{2}\rangle},\qquad
B_{\textrm{o}}(\lambda_{\textrm{o}})=\dfrac{K_{2}\langle k_{\textrm{n}}\rangle
K_{2}\langle k_{\textrm{o}}\rangle}{\langle k_{\textrm{n}}\rangle-K_{2}\langle
k_{\textrm{n}}^{2}\rangle}.$ (88)

There is an additional negative higher order term as follows:

$\displaystyle\dfrac{D_{m}(\lambda_{\textrm{o}})}{\langle
k_{\textrm{o}}\rangle[1-xT_{h}/(x_{M}T)]}(K_{2}m_{\textrm{o}})^{\lambda_{\textrm{o}}}\left[\dfrac{\lambda_{\textrm{o}}-1}{\lambda_{\textrm{o}}-3}(k^{\textrm{o}}_{\textrm{min}})^{3-\lambda_{\textrm{o}}}+2\dfrac{B_{\textrm{o}}}{K_{2}}\langle(k^{\textrm{n}})\rangle\langle(k^{\textrm{o}})^{2}\rangle+(\dfrac{B_{\textrm{o}}}{K_{2}})^{2}\langle(k^{\textrm{n}})^{2}\rangle\langle(k^{\textrm{o}})\rangle\right].$
(89)

Note that as $x\to x_{M}^{-}$, the magnitude of Eq. (89) becomes comparable to
that of the term with $C_{3}$.

  * (ii)

For $\lambda_{\textrm{n}}<\lambda_{\textrm{o}}$,

$\displaystyle f(m_{\textrm{n}})\simeq
A_{\textrm{n}}K_{2}m_{\textrm{n}}^{2}\left(1-\frac{T_{s}}{T}\right)+C_{\textrm{n}}(\lambda_{\textrm{n}})(K_{2}m_{\textrm{n}})^{\lambda_{\textrm{n}}-1}+C_{\textrm{o}}(\lambda_{\textrm{o}})(B_{\textrm{n}}m_{\textrm{n}})^{\lambda_{\textrm{o}}-1}-\frac{1}{2}\frac{K_{4}[D_{m}(\lambda_{\textrm{o}})]^{2}}{\langle
k_{\textrm{o}}\rangle[1-xT_{s}/(x_{M}T)]}(B_{\textrm{n}}m_{\textrm{n}})^{2(\lambda_{\textrm{o}}-2)}+\hbox{h.o.},$
(90)

where $A_{\textrm{n}}$ and $B_{\textrm{n}}$ are functions of
$\lambda_{\textrm{a}}$ and $K_{2}$, for which the explicit formula is as
follows:

$\displaystyle A_{\textrm{n}}(\lambda_{\textrm{n}})=\dfrac{K_{2}\langle
k_{\textrm{n}}\rangle^{2}K_{2}\langle k_{\textrm{o}}\rangle^{2}}{\langle
k_{\textrm{o}}\rangle-K_{2}\langle k_{\textrm{o}}^{2}\rangle}\qquad
B_{\textrm{n}}(\lambda_{\textrm{n}})=\dfrac{K_{2}\langle k_{\textrm{n}}\rangle
K_{2}\langle k_{\textrm{o}}\rangle}{\langle k_{\textrm{o}}\rangle-K_{2}\langle
k_{\textrm{o}}^{2}\rangle}.$ (91)

Note that $C_{a}(\lambda_{a})$ and $D_{m}(\lambda_{\textrm{o}})$ are always
positive.

The phase diagram in the space of $[x,T^{-1}]$ depends on
$\lambda_{\textrm{n}}$ and $\lambda_{\textrm{o}}$ as shown in Fig. 4. These
phase diagrams reveal the nature of the observed PTs and can be obtained by
examining the profiles of the free energy density for different $x$ and $T$
values for given $\lambda_{\textrm{n}}$ and $\lambda_{\textrm{o}}$. To be
concrete, here we consider the case of $\lambda_{\textrm{n}}=3.53$ and
$\lambda_{\textrm{o}}=3.90$, for which we obtain the phase diagram similar to
that of the original AT model with the exponent of degree distribution
$\lambda>\lambda_{c}$.

###  D.2 Case $x>x_{M}$

The self-consistency relations Eqs. (29) and (30) are expanded in terms of
$m_{a}$ and $M$ as follows:

$\displaystyle m_{\textrm{n}}\langle
k_{\textrm{n}}\rangle\left(1-K_{2}\dfrac{\langle
k_{\textrm{n}}^{2}\rangle}{\langle k_{\textrm{n}}\rangle}\right)$
$\displaystyle\simeq K_{2}m_{\textrm{o}}\langle k_{\textrm{n}}\rangle\langle
k_{\textrm{o}}\rangle-(\lambda_{\textrm{n}}-1)C_{\textrm{n}}(\lambda_{\textrm{n}})(K_{2}m_{\textrm{n}})^{\lambda_{\textrm{n}}-2}+(K_{4}M)(K_{2}m_{\textrm{n}})\langle
k_{\textrm{n}}^{2}\rangle\langle
k_{\textrm{o}}\rangle+(K_{4}M)(K_{2}m_{\textrm{o}})\langle
k_{\textrm{n}}\rangle\langle k_{\textrm{o}}^{2}\rangle+\hbox{h.o.}\,$ (92)
$\displaystyle m_{\textrm{o}}\langle
k_{\textrm{o}}\rangle\left(1-K_{2}\dfrac{\langle
k_{\textrm{o}}^{2}\rangle}{\langle k_{\textrm{o}}\rangle}\right)$
$\displaystyle\simeq K_{2}m_{\textrm{n}}\langle k_{\textrm{n}}\rangle\langle
k_{\textrm{o}}\rangle-(\lambda_{\textrm{o}}-1)C_{\textrm{o}}(\lambda_{\textrm{o}})(K_{2}m_{\textrm{o}})^{\lambda_{\textrm{o}}-2}$
$\displaystyle+D_{M}(\lambda_{\textrm{o}})(K_{4}M)^{\lambda_{\textrm{o}}-3}(K_{2}m_{\textrm{o}})-\left[\int_{0}^{1}dk_{\textrm{o}}P_{d}(k_{\textrm{o}})\tanh(K_{4}Mk_{\textrm{o}})k_{\textrm{o}}^{2}\right](K_{2}m_{\textrm{o}})\,+(K_{4}M)(K_{2}m_{\textrm{n}})\langle
k_{\textrm{n}}\rangle\langle k_{\textrm{o}}^{2}\rangle+\hbox{h.o.}\,$ (93)
$\displaystyle M\langle k_{\textrm{o}}\rangle\left(1-xK_{2}\dfrac{\langle
k_{\textrm{o}}^{2}\rangle}{\langle k_{\textrm{o}}\rangle}\right)$
$\displaystyle\simeq-(\lambda_{\textrm{o}}-1)C_{M}(\lambda_{d})(K_{4}M)^{\lambda_{\textrm{o}}-2}+(\lambda_{\textrm{o}}-3)D_{M}(\lambda_{\textrm{o}})(K_{4}M)^{\lambda_{\textrm{o}}-4}(K_{2}m_{\textrm{o}})^{2}+\hbox{h.o.,}\,$
(95)

where $D_{M}(\lambda_{\textrm{o}})>0$ increases monotonically with
$\lambda_{\textrm{o}}$. This coefficient is explicitly derived in Appendix C.
These expansions are valid for
$3<(\lambda_{\textrm{n}},\lambda_{\textrm{o}})<4$ due to the power of the
third term of the r.h.s of Eq. (93).

When $x>x_{M}$, $x\langle k_{\textrm{o}}^{2}\rangle/\langle
k_{\textrm{o}}\rangle>T_{h}$ and the l.h.s. of Eq (95) becomes negative for
$T>T_{h}$. On the other hand, the first term of the r.h.s. of Eq. (95) is also
negative; however, the second term is positive. So, the first term is
comparable to the l.h.s., leading to
$M\sim(T_{s,M}/T-1)^{1/(\lambda_{\textrm{o}}-3)}$, where $T_{s,M}\equiv
x{\langle k_{\textrm{o}}^{2}\rangle}/{\langle k_{\textrm{o}}\rangle}$. Thus,
$M$ exhibits a continuous transition at $T_{s,M}$, corresponding to the
continuous transition curve starting from CE2 in Fig. 3. Note that this
formula is the same as the one of the Ising model on a single SF network
Doro_book . For further discussions, $M_{*}$ is defined as

$M_{*}\equiv\dfrac{1}{K_{4}}\Big{[}\dfrac{\langle
k_{\textrm{o}}\rangle(T_{s,M}/T-1)}{(\lambda_{\textrm{o}}-1)K_{4}C_{M}(\lambda_{\textrm{o}})}\Big{]}^{1/(\lambda_{\textrm{o}}-3)}.$

Next, to determine a critical temperature (denoted as $T_{s,m}$) for $m_{a}$,
we first rewrite Eq. (95) as

$\displaystyle M\simeq M_{*}+\frac{D_{M}(\lambda_{\textrm{o}})}{\langle
k_{\textrm{o}}\rangle(T_{s,M}/T-1)}(K_{4}M_{*})^{\lambda_{\textrm{o}}-4}(K_{2}m_{\textrm{o}})^{2}+\hbox{h.o}.$
(96)

We consider the linear terms of $m_{a}$ in Eqs. 92 and 93, and substitute $M$
with $M_{*}$. Using a similar technique used in Eq. (83), we obtain the
following:

$\displaystyle\langle k_{\textrm{n}}\rangle\langle
k_{\textrm{o}}\rangle\left(1-\dfrac{\langle k_{\textrm{n}}^{2}\rangle/\langle
k_{\textrm{n}}\rangle+g_{\textrm{n}}(M_{*})}{T}\right)\left(1-\dfrac{\langle
k_{\textrm{o}}^{2}\rangle/\langle
k_{\textrm{o}}\rangle+g_{\textrm{o}}(M_{*})}{T}\right)-\left(\dfrac{\langle
k_{\textrm{n}}\rangle\langle
k_{\textrm{o}}\rangle}{T}+\dfrac{K_{4}M_{*}\langle
k_{\textrm{n}}\rangle\langle k_{\textrm{o}}^{2}\rangle}{T}\right)^{2}=0\,,$
(97)

where

$\displaystyle g_{\textrm{n}}(M_{*})\langle
k_{\textrm{n}}\rangle=K_{4}M_{*}\langle k_{\textrm{o}}\rangle\langle
k_{\textrm{n}}^{2}\rangle\,,\quad\textrm{and}\quad
g_{\textrm{o}}(M_{*})\langle
k_{\textrm{o}}\rangle=D_{M}(\lambda_{\textrm{o}})(K_{4}M_{*})^{\lambda_{\textrm{o}}-3}-\int_{0}^{1}dk_{\textrm{o}}P_{d}(k_{\textrm{o}})\tanh(K_{4}M_{*}k_{\textrm{o}})k_{\textrm{o}}^{2}\,.$

Eq. (97) has two solutions for $T$, denoted as $T_{\ell}^{\prime}$ and
$T_{h}^{\prime}$ ($T_{\ell}^{\prime}<T_{h}^{\prime}$). Using the relation
(96), we can obtain a self-consistency relation for $m_{a}$ within the leading
order as follows:

$\displaystyle m_{a}\langle k_{\textrm{o}}\rangle\langle
k_{\textrm{n}}\rangle\left(1-\frac{T_{\ell}^{\prime}}{T}\right)\left(1-\frac{T_{h}^{\prime}}{T}\right)\simeq-(\lambda_{a}-1)E_{\bar{a}}^{\prime}C_{a}(\lambda_{a})\,(K_{2}m_{a})^{\lambda_{a}-2}-(\lambda_{\bar{a}}-1)F_{a\bar{a}}^{\prime}C_{\bar{a}}(\lambda_{\bar{a}})\,(K_{2}m_{\bar{a}})^{\lambda_{\bar{a}}-2}\,,$
(98)

where $E_{\bar{a}}^{\prime}=\langle k_{\bar{a}}\rangle\left(1-\dfrac{\langle
k_{\bar{a}}^{2}\rangle/\langle
k_{\bar{a}}\rangle+g_{\bar{a}}(M_{*})}{T}\right)$ and
$F_{a\bar{a}}^{\prime}=\dfrac{1}{T}\langle k_{\textrm{o}}\rangle\langle
k_{\textrm{n}}\rangle+\dfrac{1}{T}K_{4}M_{*}\langle
k_{\textrm{o}}^{2}\rangle\langle k_{\textrm{n}}\rangle$. We find that near
$T_{h}^{\prime}$, $m_{a}$ converges to zero continuously, whereas $M$ remains
in $\mathcal{O}(1)$. Hence, we regard $T_{h}^{\prime}$ as the critical
temperature $T_{s}^{\prime}$ of $m_{a}$. Note that $M$ has the critical
temperature $T_{s,M}$ separately, given as $x\langle
k_{\textrm{o}}^{2}\rangle/\langle k_{\textrm{o}}\rangle$, which is higher than
$T_{s}^{\prime}$.

Using Eq. (96), we expand the free energy density of Eq. (26) with respect to
$m_{a}$ up to the three lowest order terms:

  * (i)

For $\lambda_{\textrm{n}}>\lambda_{\textrm{o}}$,

$\displaystyle f(m_{\textrm{o}})$ $\displaystyle\simeq
f_{0}(M_{*})+A_{\textrm{o}}^{\prime}K_{2}m_{\textrm{o}}^{2}\left(1-\frac{T_{s}^{\prime}}{T}\right)+C_{\textrm{o}}(\lambda_{\textrm{o}})(K_{2}m_{\textrm{o}})^{\lambda_{\textrm{o}}-1}+C_{\textrm{n}}(\lambda_{\textrm{n}})(B_{\textrm{o}}^{\prime}m_{\textrm{o}})^{\lambda_{\textrm{n}}-1}$
(99)
$\displaystyle-\frac{K_{2}[(\lambda_{\textrm{n}}-1)C_{\textrm{n}}(\lambda_{\textrm{n}})]^{2}}{\langle
k_{\textrm{n}}\rangle\left[1-\left(\langle k_{\textrm{n}}^{2}\rangle/\langle
k_{\textrm{n}}\rangle+g_{\textrm{n}}(M_{*})\right)/T\right]}(B_{\textrm{o}}^{\prime}m_{\textrm{o}})^{2(\lambda_{\textrm{n}}-2)}-\frac{1}{2}(\lambda_{\textrm{o}}-3)\dfrac{K_{4}\left[D_{M}(\lambda_{\textrm{o}})(K_{4}M_{*})^{\lambda_{\textrm{o}}-4}\right]^{2}}{\langle
k_{\textrm{o}}\rangle\left(T_{s,M}/T-1\right)}(K_{2}m_{\textrm{o}})^{4}+\hbox{h.o.}\,,$
(100)

where $A_{\textrm{o}}^{\prime}$ and $B_{\textrm{o}}^{\prime}$ are functions of
$\lambda_{a}$, $K_{2}$ and $K_{4}M_{*}$, They are explicitly derived as
follows:

$\displaystyle
A_{\textrm{o}}^{\prime}(\lambda_{\textrm{o}})=\dfrac{\left(1+K_{4}M_{*}\langle
k_{\textrm{o}}^{2}\rangle/\langle k_{\textrm{o}}\rangle\right)^{2}K_{2}\langle
k_{\textrm{n}}\rangle^{2}K_{2}\langle k_{\textrm{o}}\rangle^{2}}{\langle
k_{\textrm{n}}\rangle\left[1-K_{2}\left(\langle
k_{\textrm{n}}^{2}\rangle/\langle
k_{\textrm{n}}\rangle+g_{\textrm{n}}(M_{*})\right)\right]},\qquad
B_{\textrm{o}}^{\prime}(\lambda_{\textrm{o}})=\dfrac{\left(1+K_{4}M_{*}\langle
k_{\textrm{o}}^{2}\rangle/\langle k_{\textrm{o}}\rangle\right)K_{2}\langle
k_{\textrm{n}}\rangle K_{2}\langle k_{\textrm{o}}\rangle}{\langle
k_{\textrm{n}}\rangle\left[1-K_{2}\left(\langle
k_{\textrm{n}}^{2}\rangle/\langle
k_{\textrm{n}}\rangle+g_{\textrm{n}}(M_{*})\right)\right]}.$ (101)

There is an additional negative higher order term as follows:

$\displaystyle\dfrac{D_{M}(\lambda_{\textrm{o}})}{\langle
k_{\textrm{o}}\rangle(T_{s,M}/T-1)}(K_{4}M_{*})^{\lambda_{\textrm{o}}-4}(K_{2}m_{\textrm{o}})^{4}\left[\dfrac{\lambda_{\textrm{o}}-1}{\lambda_{\textrm{o}}-3}(k^{\textrm{o}}_{\textrm{min}})^{3-\lambda_{\textrm{o}}}+2\dfrac{B_{\textrm{o}}^{\prime}}{K_{2}}\langle(k^{\textrm{n}})\rangle\langle(k^{\textrm{o}})^{2}\rangle+(\dfrac{B_{\textrm{o}}^{\prime}}{K_{2}})^{2}\langle(k^{\textrm{n}})^{2}\rangle\langle(k^{\textrm{o}})\rangle\right]$
(102)

Note that as $x\to x_{M}^{+}$, Eq. (89) becomes close in its magnitude to the
term with $C_{3}$ and thus, play a similar role to the term with $C_{3}$ near
$x_{M}^{+}$.

  * (ii)

For $\lambda_{\textrm{n}}<\lambda_{\textrm{o}}$,

$\displaystyle f(m_{\textrm{n}})$ $\displaystyle\simeq
f_{0}(M_{*})+A_{\textrm{n}}^{\prime}K_{2}m_{\textrm{n}}^{2}\left(1-\frac{T_{s}^{\prime}}{T}\right)+C_{\textrm{n}}(\lambda_{\textrm{n}})(K_{2}m_{\textrm{n}})^{\lambda_{\textrm{n}}-1}+C_{\textrm{o}}(\lambda_{\textrm{o}})(B_{\textrm{n}}^{\prime}m_{\textrm{n}})^{\lambda_{\textrm{o}}-1}$
(103)
$\displaystyle-\frac{K_{2}[(\lambda_{\textrm{o}}-1)C_{\textrm{o}}(\lambda_{\textrm{n}})]^{2}}{\langle
k_{\textrm{o}}\rangle\left[1-\left(\langle k_{\textrm{o}}^{2}\rangle/\langle
k_{\textrm{o}}\rangle+g_{\textrm{o}}(M_{*})\right)/T\right]}(B_{\textrm{n}}^{\prime}m_{\textrm{n}})^{2(\lambda_{\textrm{o}}-2)}-\frac{1}{2}(\lambda_{\textrm{o}}-3)\dfrac{K_{4}\left[D_{M}(\lambda_{\textrm{o}})(K_{4}M_{*})^{\lambda_{\textrm{o}}-4}\right]^{2}}{\langle
k_{\textrm{o}}\rangle\left(T_{s,M}/T-1\right)}(B_{\textrm{n}}^{\prime}m_{\textrm{n}})^{4}+\hbox{h.o.}\,,$
(104)

where $A_{\textrm{n}}^{\prime}$ and $B_{\textrm{n}}^{\prime}$ are functions of
$\lambda_{a}$, $K_{2}$ and $K_{4}M_{*}$. They are explicitly derived as
follows:

$\displaystyle
A_{\textrm{n}}^{\prime}(\lambda_{\textrm{o}})=\dfrac{\left(1+K_{4}M_{*}\langle
k_{\textrm{o}}^{2}\rangle/\langle k_{\textrm{o}}\rangle\right)^{2}K_{2}\langle
k_{\textrm{n}}\rangle^{2}K_{2}\langle k_{\textrm{o}}\rangle^{2}}{\langle
k_{\textrm{o}}\rangle\left[1-K_{2}\left(\langle
k_{\textrm{o}}^{2}\rangle/\langle
k_{\textrm{o}}\rangle+g_{\textrm{o}}(M_{*})\right)\right]},\qquad
B_{\textrm{n}}^{\prime}(\lambda_{\textrm{o}})=\dfrac{\left(1+K_{4}M_{*}\langle
k_{\textrm{o}}^{2}\rangle/\langle k_{\textrm{o}}\rangle\right)K_{2}\langle
k_{\textrm{n}}\rangle K_{2}\langle k_{\textrm{o}}\rangle}{\langle
k_{\textrm{o}}\rangle\left[1-K_{2}\left(\langle
k_{\textrm{o}}^{2}\rangle/\langle
k_{\textrm{o}}\rangle+g_{\textrm{o}}(M_{*})\right)\right]}.$ (105)

Here, first two $C_{a}$ terms are positive, like the case $x<x_{M}$,

The first two $C_{a}$ terms in Eqs.(100) and (104) are positive, as for the
case $x<x_{M}$, whereas the next two terms containing $C_{a}$ and $D_{M}$ are
negative. The $2(\lambda_{a}-2)$-order terms with $C_{a}$ are finite, whereas
the $D_{M}$ term diverges as $T\to T_{s,M}$. Thus, the $D_{M}$ term
contributes to the formation of a global minimum of $f(m_{a})$ as $T$ is
decreased $T_{f}$ and $x\to x_{M}^{+}$.

###  D.3 Case $x\approx x_{M}$

####  D.3.1 Case $\lambda_{\rm n}\geq\lambda_{\rm o}$

In this case, $O(m_{a})\sim O(M)$ near $T_{s}$ and the free energy density of
Eq. (26) is expanded with respect to $m_{\textrm{o}}$ as follows:

$\displaystyle f(m_{\textrm{o}})$ $\displaystyle\simeq
A_{\textrm{o}}K_{2}m_{\textrm{o}}^{2}\left(1-\frac{T_{s}}{T}\right)+\frac{1}{2}K_{4}M^{2}\langle
k_{\textrm{o}}\rangle\left(1-\frac{T_{s}}{T}\right)+C_{\textrm{o}}(\lambda_{\textrm{o}},r_{\textrm{o}})(K_{2}m_{\textrm{o}})^{\lambda_{\textrm{o}}-1}+C_{\textrm{n}}(\lambda_{\textrm{n}})(B_{\textrm{o}}m_{\textrm{o}})^{\lambda_{\textrm{n}}-1}+C_{3}(\lambda_{\textrm{n}},\lambda_{\textrm{o}},r_{0})(K_{2}m_{\textrm{o}})^{3}+\hbox{h.o.}\,,$
(106)

where $A_{\textrm{o}}$ and $B_{\textrm{o}}$ are functions of $\lambda_{a}$ and
$K_{2}$ that are explicitly derived in the $x<x_{M}$ case.
$C_{\textrm{o}}(\lambda_{\textrm{o}},r_{0})$ with $r_{0}\equiv
K_{4}M/K_{2}m_{\textrm{o}}$ is $\mathcal{O}(1)$. Explicit formulas of the
coefficients are given in Appendix C. There is another negative term with
$C_{3}$, which is defined as follows:

$\displaystyle
C_{3}(r_{0},\lambda_{\textrm{o}},\lambda_{\textrm{n}})=r_{0}\left[\dfrac{\lambda_{\textrm{o}}-1}{\lambda_{\textrm{o}}-3}(k^{\textrm{o}}_{\textrm{min}})^{3-\lambda_{\textrm{o}}}+2\dfrac{B_{\textrm{o}}}{K_{2}}\langle(k^{\textrm{n}})\rangle\langle(k^{\textrm{o}})^{2}\rangle+(\dfrac{B_{\textrm{o}}}{K_{2}})^{2}\langle(k^{\textrm{n}})^{2}\rangle\langle(k^{\textrm{o}})\rangle\right].$
(107)

Note that this term does not appear in the original AT model defined on SF
network.

We note that the $C_{\textrm{o}}$ term is a leading order term at $T_{s}$, and
$C_{\textrm{o}}$ decreases monotonically with $\lambda_{\textrm{o}}$. Thus the
sign of $C_{\textrm{o}}$ can change depending on the magnitude of
$\lambda_{\textrm{o}}$. This feature does not appear for both $0<x<x_{M}$ and
$x>x_{M}$ cases. However, it occurs when $x=x_{M}$. Numerically
$C_{\textrm{o}}=0$ at $\lambda_{\textrm{o}}\approx 3.503$, equivalent to
$\lambda_{c}$ introduced earlier in Secs. II and III: $C_{\textrm{o}}$ becomes
positive for $\lambda_{\textrm{o}}<\lambda_{c}$, whereas it is negative for
$\lambda_{\textrm{o}}>\lambda_{c}$. On the other hand, $C_{\textrm{n}}$ and
$C_{3}$ are always positive and negative, respectively.

####  D.3.2 Case $\lambda_{\rm n}=\lambda_{\rm o}$

When $\lambda_{\textrm{n}}=\lambda_{\textrm{o}}$, the $C_{\textrm{o}}$ and
$C_{\textrm{n}}$ terms in Eq. (106) are of the same order. Thus, the two terms
are combined and denoted as
$C_{\textrm{o}}^{\prime}(\lambda_{\textrm{o}},r_{0})(K_{2}m_{\textrm{o}})^{\lambda_{\textrm{o}}-1}$.
The sign of $C_{\textrm{o}}^{\prime}$ depends on $\lambda_{\textrm{o}}$,
equivalently $\lambda_{\textrm{n}}$. Numerically $C_{\textrm{o}}^{\prime}$ can
be zero at a certain $\lambda_{\textrm{o}}$, denoted as $\lambda_{e}$,
estimated to be $\approx 3.605$. $C_{\textrm{o}}^{\prime}$ becomes positive
for $\lambda_{\textrm{o}}<\lambda_{e}$ and negative otherwise.

Similar to the previous case $\lambda_{\textrm{n}}>\lambda_{\textrm{o}}$, a
discontinuous transition always occurs for $C_{\textrm{o}}^{\prime}<0$.
However, depending on relative magnitude between two terms
$C_{\textrm{o}}^{\prime}$ and $C_{3}$, either a discontinuous or continuous
transition occurs for $C_{\textrm{o}}^{\prime}>0$. Note that successive
discontinuous transitions do not occur for
$\lambda_{\textrm{n}}=\lambda_{\textrm{o}}$.

####  D.3.3 Case $\lambda_{\rm n}<\lambda_{\rm o}$

When $\lambda_{\textrm{n}}<\lambda_{\textrm{o}}$,
$\mathcal{O}(m_{a})\ll\mathcal{O}(M)$, the free energy density of Eq. (26) is
expanded with respect to $m_{\textrm{n}}$ as follows:

$\displaystyle f(m_{\textrm{n}})$ $\displaystyle\simeq
f_{0}(M_{*})+A_{\textrm{n}}^{\prime}K_{2}m_{\textrm{n}}^{2}\left(1-\frac{T_{s}}{T}\right)+C_{\textrm{n}}(\lambda_{\textrm{n}})(K_{2}m_{\textrm{n}})^{\lambda_{\textrm{n}}-1}+C_{\textrm{o}}(\lambda_{\textrm{o}})(B_{\textrm{n}}^{\prime}m_{\textrm{n}})^{\lambda_{\textrm{o}}-1}$
(108)
$\displaystyle-\frac{K_{2}[(\lambda_{\textrm{o}}-1)C_{\textrm{o}}(\lambda_{\textrm{n}})]^{2}}{\langle
k_{\textrm{o}}\rangle\left[1-\left(\langle k_{\textrm{o}}^{2}\rangle/\langle
k_{\textrm{o}}\rangle+g_{\textrm{o}}(M_{*})\right)/T\right]}(B_{\textrm{n}}^{\prime}m_{\textrm{n}})^{2(\lambda_{\textrm{o}}-2)}-\frac{1}{2}(\lambda_{\textrm{o}}-3)\frac{K_{4}\left[D_{M}(\lambda_{\textrm{o}})(K_{4}M_{*})^{\lambda_{\textrm{o}}-4}\right]^{2}}{\langle
k_{\textrm{o}}\rangle\left(T_{s,M}/T-1\right)}(B_{\textrm{n}}^{\prime}m_{\textrm{n}})^{4}+\hbox{h.o.}.\,$

This formula is exactly the same as Eq. (104), derived in $x>x_{M}$ case for
continuous transitions.

##  Appendix E The susceptibility near the critical temperature

###  E.1 $m_{a}$ \- magnetization

####  E.1.1 for $x\leq x_{M}$

Now, we consider the susceptibility at the critical temperature $T_{s}$ for
weak interlayer interaction $x<x_{M}$ case. We can omit the higher order terms
in $m_{a}$ and $M$ when $m_{a}$ and $M$ are very small, we expand the self-
consistency relations for $m_{a}$ (61) with respect to $m_{a}$ and $M$ as
follows:

$\displaystyle m_{a}\langle k_{a}\rangle$
$\displaystyle\simeq\Omega_{a}\langle k_{a}^{2}\rangle+\Omega_{\bar{a}}\langle
k_{\textrm{o}}\rangle\langle
k_{\textrm{n}}\rangle-(\lambda_{a}-1)C_{a}(\lambda_{a})(\Omega_{a})^{\lambda_{a}-2}\,.$
(109)

To obtain critical exponent $\gamma$ for each $m_{a}$-magnetization,
repectively, we consider the lowest-order terms of the self-consistency
relations Eqs. (109) and then we obtain the following:

$\displaystyle m_{a}\langle k_{\textrm{n}}\rangle\langle
k_{\textrm{o}}\rangle\left(1-\dfrac{T_{\ell}}{T}\right)\left(1-\dfrac{T_{s}}{T}\right)$
$\displaystyle\simeq E_{\bar{a}}\left(H_{a}\langle
k_{a}^{2}\rangle+H_{\bar{a}}\langle k_{\textrm{o}}\rangle\langle
k_{\textrm{n}}\rangle-(\lambda_{a}-1)C_{a}(\lambda_{a})(\Omega_{a})^{\lambda_{a}-2}\right)$
(110) $\displaystyle+F_{a\bar{a}}\left(H_{\bar{a}}\langle
k_{\bar{a}}^{2}\rangle+H_{a}\langle k_{a}\rangle\langle
k_{\bar{a}}\rangle-(\lambda_{\bar{a}}-1)C_{\bar{a}}(\lambda_{\bar{a}})(\Omega_{\bar{a}})^{\lambda_{\bar{a}}-2}\right)+\hbox{h.o.}\,,$
(111)

where $E_{\bar{a}}=\langle k_{\bar{a}}\rangle\left(1-\dfrac{\langle
k_{\bar{a}}^{2}\rangle/\langle k_{\bar{a}}\rangle}{T}\right)$ and
$F_{a\bar{a}}=\dfrac{1}{T}\langle k_{\textrm{o}}\rangle\langle
k_{\textrm{n}}\rangle$.

In order to derive the susceptilbility from magnetization, we take partial
derivative with respect to $H_{a}$ and then take the $H_{a}$ and
$H_{4}\rightarrow 0$ limit. We have two equations for the susceptibility as
follows:

$\displaystyle\chi_{a}\langle k_{a}\rangle\langle
k_{\bar{a}}\rangle(1-\dfrac{T_{\ell}}{T})(1-\dfrac{T_{s}}{T})$
$\displaystyle\simeq E_{\bar{a}}\left(\langle
k_{a}^{2}\rangle-(\lambda_{a}-1)(\lambda_{a}-2)C_{a}(\lambda_{a})(K_{2}m_{a})^{\lambda_{a}-3}K_{2}\chi_{a}\right)$
(112) $\displaystyle+F\left(\langle k_{a}\rangle\langle
k_{\bar{a}}\rangle-(\lambda_{\bar{a}}-1)(\lambda_{\bar{a}}-2)C_{\bar{a}}(\lambda_{\bar{a}})(K_{2}m_{\bar{a}})^{\lambda_{\bar{a}}-3}K_{2}\dfrac{\partial
m_{\bar{a}}}{\partial H_{2}^{a}}\right)+\hbox{h.o.}\,.$ (113)

When $\lambda_{\textrm{o}}<\lambda_{\textrm{n}}$, the susceptibility near the
critical temperature $T_{s}^{-}$ is written as follows:

$\displaystyle\chi_{\textrm{o}}$ $\displaystyle\langle
k_{\textrm{n}}\rangle\langle
k_{\textrm{o}}\rangle(1-\dfrac{T_{\ell}}{T})(1-\dfrac{T_{s}}{T})\approx
E_{\textrm{n}}\langle k_{\textrm{o}}^{2}\rangle+F\langle
k_{\textrm{o}}\rangle\langle
k_{\textrm{n}}\rangle+(\lambda_{\textrm{o}}-2)\langle
k_{\textrm{o}}\rangle\langle
k_{\textrm{n}}\rangle(1-\dfrac{T_{\ell}}{T})(1-\dfrac{T_{s}}{T})\chi_{\textrm{o}}+\hbox{h.o.}\,,$
(114) $\displaystyle\chi_{\textrm{n}}$ $\displaystyle\langle
k_{\textrm{n}}\rangle\langle
k_{\textrm{o}}\rangle(1-\dfrac{T_{\ell}}{T})(1-\dfrac{T_{s}}{T})\approx
E_{\textrm{o}}\langle k_{\textrm{n}}^{2}\rangle+F\langle
k_{\textrm{o}}\rangle\langle
k_{\textrm{n}}\rangle+(\lambda_{\textrm{o}}-2)\langle
k_{\textrm{o}}\rangle\langle
k_{\textrm{n}}\rangle(1-\dfrac{T_{\ell}}{T})(1-\dfrac{T_{s}}{T})\dfrac{E_{\textrm{o}}}{E_{\textrm{n}}}\dfrac{\partial
m_{\textrm{o}}}{\partial H_{\textrm{n}}}+\hbox{h.o.}\,.$ (115)

To obtain the susceptibility near $T_{s}^{-}$, we use the following relation

$\displaystyle
E_{\textrm{n}}\left(-(\lambda_{\textrm{o}}-1)K_{2}C_{\textrm{o}}(\lambda_{\textrm{o}})(K_{2}m_{\textrm{o}})^{\lambda_{\textrm{o}}-3}\right)\approx{\langle
k_{\textrm{o}}\rangle\langle
k_{\textrm{n}}\rangle}(1-\dfrac{T_{\ell}}{T})(1-\dfrac{T_{s}}{T})$
$\displaystyle~{}~{}\hbox{ for }~{}~{}T\to T_{s}^{-}.$ (116)

To get $\chi_{\textrm{n}}$, we need to compute the $\partial
m_{\textrm{o}}/\partial H_{\textrm{n}}$ term. The partial derivative of
$m_{\textrm{o}}$ in terms of $H_{\textrm{n}}$ is given as follows:

$\displaystyle\dfrac{\partial m_{\textrm{o}}}{\partial H_{\textrm{n}}}\langle
k_{\textrm{o}}\rangle\langle
k_{\textrm{n}}\rangle(1-\dfrac{T_{\ell}}{T})(1-\dfrac{T_{s}}{T})$
$\displaystyle\simeq E_{\textrm{n}}\left(\langle k_{\textrm{o}}\rangle\langle
k_{\textrm{n}}\rangle+(\lambda_{\textrm{o}}-1)(\lambda_{\textrm{o}}-2)C_{\textrm{o}}(\lambda_{\textrm{o}})(K_{2}m_{\textrm{o}})^{\lambda_{\textrm{o}}-3}K_{2}\dfrac{\partial
m_{\textrm{o}}}{\partial H_{\textrm{n}}}\right)+F\langle
k_{\textrm{n}}^{2}\rangle+\hbox{h.o.}$ (117) $\displaystyle\simeq
E_{\textrm{n}}\langle k_{\textrm{o}}\rangle\langle
k_{\textrm{n}}\rangle+F\langle
k_{\textrm{n}}^{2}\rangle+(\lambda_{\textrm{o}}-2)\langle
k_{\textrm{o}}\rangle\langle
k_{\textrm{n}}\rangle(1-\dfrac{T_{\ell}}{T})(1-\dfrac{T_{s}}{T})\dfrac{\partial
m_{\textrm{o}}}{\partial H_{\textrm{n}}}+\hbox{h.o.}\,.$ (118)

From Eq. (118), we obtain that $\partial m_{\textrm{o}}/\partial
H_{\textrm{n}}\approx(T_{s}-T)^{-1}$ for $T\to T_{s}^{-}$. Using this, the
susceptibility of $m_{a}$ is obtained as

$\displaystyle\chi_{\textrm{o}}\approx\left\\{\begin{array}[]{cc}(T-T_{s})^{-1}&~{}~{}\textrm{for}~{}T_{s}^{+},\\\
\\\
(T_{s}-T)^{-1}&~{}~{}\textrm{for}~{}T_{s}^{-},\end{array}\right.\quad\quad\hbox{and}\quad\quad\chi_{\textrm{n}}\approx\left\\{\begin{array}[]{cc}(T-T_{s})^{-1}&~{}~{}\textrm{for}~{}T_{s}^{+},\\\
\\\ (T_{s}-T)^{-1}&~{}~{}\textrm{for}~{}T_{s}^{-}.\end{array}\right.\,$ (125)

Here, we take the limit $m_{\textrm{o}}\to 0$ for near $T_{s}^{+}$ to Eq.
(113). For the $\lambda_{\textrm{n}}<\lambda_{\textrm{o}}$ case, the
susceptibility of $m_{a}$ is obtained by similar computation process same as
for the $\lambda_{\textrm{o}}<\lambda_{\textrm{n}}$ case. Thus, the critical
exponent of $\gamma_{m\pm}$ of magnetization is always $1$ for all cases.
Then, the scaling relation, $a+2\beta_{m}+\gamma_{m-}=2$, is satisfied for
each $m_{a}$-magnetization, respectively.

Now, we compute the susceptibility at CE point as boundary point of continuous
PT regime. Since the location of magnetization jumps to a certain finite at
the CE, the magnitude of the magnetization is much greater than 0, the
perturbative expansions with respect to $m,M$ is not valid any longer at CE
point. Thus, we should keep the integral formula written in self-consistent
relation Eqs. (61) and (62) as follows. In order to obtain the susceptibility
for $m_{a}$ magnetization, we take a partial derivative of self-consistent
relation for $m_{a}$ (61) with respect ot $H_{2}$ and take
$H_{2},H_{4}\rightarrow 0$ limit, then the susceptibility is written as
follows:

$\displaystyle\chi_{\textrm{o}}=\dfrac{\mathcal{A}_{\textrm{o}\textrm{o}}+\mathcal{A}_{\textrm{o}\textrm{n}}K_{2}\partial
m_{\textrm{n}}/\partial H_{\textrm{o}}+\mathcal{A}_{\textrm{o}M}K_{4}\partial
M/\partial H_{\textrm{o}}}{\langle k_{\textrm{o}}\rangle-
K_{2}\mathcal{A}_{\textrm{o}\textrm{o}}}\,.$ (126)

To evaluate Eq. (126), we also should compute the

$\displaystyle\left.\dfrac{\partial m_{\textrm{n}}}{\partial
H_{\textrm{o}}}\right|_{H_{a},H_{4}\rightarrow 0}\textrm{and
}\left.\dfrac{\partial M}{\partial
H_{\textrm{o}}}\right|_{H_{a},H_{4}\rightarrow 0}$ (127)

terms. Thus, we firstly take a derivative of self-consistent relation for
$m_{\textrm{n}}$ and $M$ (61),(62) with respect to $H_{\textrm{o}}$ and then
take the limit $H_{a}$ and $H_{4}\rightarrow 0$, we obtain as follows:

$\displaystyle\dfrac{\partial m_{\textrm{n}}}{\partial
H_{\textrm{o}}}=\dfrac{\mathcal{A}_{\textrm{n}\textrm{o}}+\mathcal{A}_{\textrm{n}\textrm{o}}K_{2}\chi_{\textrm{o}}+\mathcal{A}_{\textrm{n}M}K_{4}\partial
M/\partial H_{\textrm{o}}}{\langle k_{\textrm{n}}\rangle-
K_{2}\mathcal{A}_{\textrm{n}\textrm{n}}},\quad\dfrac{\partial M}{\partial
H_{\textrm{o}}}=\dfrac{\mathcal{A}_{M\textrm{o}}+\mathcal{A}_{M\textrm{o}}K_{2}\chi_{\textrm{o}}+\mathcal{A}_{M\textrm{n}}K_{2}\partial
m_{\textrm{n}}/\partial H_{\textrm{o}}}{\langle k_{\textrm{o}}\rangle-
K_{4}\mathcal{A}_{MM}}.$ (128)

At the CE point, $\chi_{\textrm{o}}$ is computed similarly to Eq. (113) at
$T_{s}^{+}$, where $m_{a}=M=0$. For $T_{s}^{-}$, $\chi_{m}$ can be obtained
numerically from Eqs. (126), (128). We can confirm that the susceptibility has
a certain finite value at $T_{s}^{-}$ by numerical computations.

####  E.1.2 Case $x>x_{M}$

Otherwise, for $x>x_{M}$, we expand the self-consistency relations for $m_{a}$
with respect to $m_{a}$ and $M$ as follows:

$\displaystyle m_{a}\langle k_{\textrm{n}}\rangle\langle
k_{\textrm{o}}\rangle(1-\dfrac{T_{\ell}^{\prime}}{T})(1-\dfrac{T_{s}^{\prime}}{T})$
$\displaystyle\simeq E_{\bar{a}}^{\prime}\left[H_{a}\langle
k_{a}^{2}\rangle+H_{\bar{a}}\langle k_{\textrm{o}}\rangle\langle
k_{\textrm{n}}\rangle-(\lambda_{a}-1)C_{a}(\lambda_{a})(\Omega_{a})^{\lambda_{a}-2}\right]$
(129) $\displaystyle+F_{a\bar{a}}^{\prime}\left[H_{\bar{a}}\langle
k_{\bar{a}}^{2}\rangle+H_{a}\langle k_{a}\rangle\langle
k_{\bar{a}}\rangle-(\lambda_{\bar{a}}-1)C_{\bar{a}}(\lambda_{\bar{a}})(\Omega_{\bar{a}})^{\lambda_{\bar{a}}-2}\right]+\hbox{h.o.}\,,$
(130)

where $E_{\bar{a}}^{\prime}=\langle k_{\bar{a}}\rangle\left(1-\dfrac{\langle
k_{\bar{a}}^{2}\rangle/\langle
k_{\bar{a}}\rangle+g_{\bar{a}}(M_{*})}{T}\right)$ and
$F_{a\bar{a}}^{\prime}=\dfrac{1}{T}\langle k_{\textrm{o}}\rangle\langle
k_{\textrm{n}}\rangle+\dfrac{1}{T}K_{4}M_{*}\langle
k_{\textrm{o}}^{2}\rangle\langle k_{\textrm{n}}\rangle$.

It can be checked easily that Eq. (130) is similar case to Eq. (111) except
the critical temperature $T_{s}^{\prime}$ and coefficients
$E_{\bar{a}}^{\prime}$ and $F^{\prime}$. Thus, we performed similar
calculations as for the $x<x_{M}$ case considering minor differences between
the cases $x<x_{M}$ and $x>x_{M}$. By performing similar calculations as for
the $x<x_{M}$ case, we obtain the susceptibility of $m_{a}$ as follows:

$\displaystyle\chi_{\textrm{o}}\approx\left\\{\begin{array}[]{cc}(T-T_{s}^{\prime})^{-1}&~{}~{}\textrm{for}~{}{T>T_{s}^{\prime}},\\\
\\\
(T_{s}^{\prime}-T)^{-1}&~{}~{}\textrm{for}~{}{T_{s}^{\prime}>T},\end{array}\right.\quad\quad\hbox{and}\quad\quad\chi_{\textrm{n}}\approx\left\\{\begin{array}[]{cc}(T-T_{s}^{\prime})^{-1}&~{}~{}\textrm{for}~{}{T>T_{s}^{\prime}},\\\
\\\
(T_{s}^{\prime}-T)^{-1}&~{}~{}\textrm{for}~{}{T_{s}^{\prime}>T}.\end{array}\right.\,$
(137)

###  E.2 $M$ \- magnetization

Likewise, the self-consistency relation for $M$ (62) can be expanded as

$\displaystyle M\langle k_{\textrm{o}}\rangle\simeq\Omega_{4}\langle
k_{\textrm{o}}^{2}\rangle-(\lambda_{\textrm{o}}-1)C_{M}(\lambda_{\textrm{o}})\Omega_{4}^{\lambda_{\textrm{o}}-2}-D_{m}(\lambda_{\textrm{o}})\Omega_{\lambda_{\textrm{o}}-2}+\hbox{h.o.}.$
(138)

To obtain the susceptibility of $M$, we take partial derivative of the above
self-consistency relation with respect to $H_{4}$ and then taking $H_{2}$ and
$H_{4}\rightarrow 0$.

$\displaystyle\chi_{M}\langle
k_{\textrm{o}}\rangle\simeq(K_{4}\chi_{M}+1)\langle
k_{\textrm{o}}^{2}\rangle-(\lambda_{\textrm{o}}-2)C_{M}(\lambda_{\textrm{o}})(K_{4}\chi_{M})(K_{4}M)^{\lambda_{\textrm{o}}-3}-(\lambda_{\textrm{o}}-2)D_{m}(\lambda_{\textrm{o}})(K_{2}m_{\textrm{o}})^{\lambda_{\textrm{o}}-3}K_{2}\dfrac{\partial
m_{\textrm{o}}}{\partial H_{4}}+\hbox{h.o.}\,.$ (139)

For $x<x_{M}$, because $M$ is
$\mathcal{O}(m_{\textrm{o}}^{\lambda_{\textrm{o}}-2})$, $\partial
m_{\textrm{o}}/\partial H_{4}$ is very small compared with the
$\mathcal{O}(1)$ term. Taking this limit, we obtain the susceptibility of $M$.
We also take the limit $m_{\textrm{o}}=0$ for $T_{s}^{+}$.

$\displaystyle\chi_{M}\approx\left(T-x\dfrac{\langle
k_{\textrm{o}}^{2}\rangle}{\langle k_{\textrm{o}}\rangle}\right)^{-1}.$ (140)

Otherwise, $x\geq x_{M}$ and $T\to T_{s,M}^{-}$, $M$ can be approximated to
$M_{*}$, where $M_{*}$ becomes Ising spin in single SF networks, and $m_{a}$
is negligible to $M$. Taking this limit, we can obtain the susceptibility of
$M$ as follow. We also take the limit $M=0$ for $T_{s,M}^{+}$.

$\displaystyle\chi_{M}\approx\left\\{\begin{array}[]{cc}(T-T_{s,M})^{-1}&~{}~{}\textrm{for}~{}~{}T>T_{s,M},\\\
\\\ (T_{s,M}-T)^{-1}&~{}~{}\textrm{for}~{}~{}T_{s,M}>T.\end{array}\right.$
(144)

###### Acknowledgements.

This research was supported by the NRF Grant (No. NRF-2014R1A3A2069005), the
KENTECH Research Grant (KRG2021-01-007) (BK), and the KIAS individual Grants
(No. PG064901) (JSL) at Korea Institute for Advanced Study.

## References

  * (1) A.-L. Barabási, Network Science (Cambridge, Cambridge University Press, 2016).
  * (2) S. N. Dorogovtsev, Lectures on complex networks (Oxford, Oxford University Press, 2010).
  * (3) M. E. J. Newman, Networks: an introduction (Oxford, Oxford University Press, 2010).
  * (4) D. Lee, B. Kahng, Y.S. Cho, K.-I. Goh, and D.-S. Lee, J. Korean Phys. Soc. 73, 152-164 (2018).
  * (5) S. N. Dorogovtsev, A. V. Goltsev, and J. F. F. Mendes, Rev. Mod. Phys. 80, 1275 (2008).
  * (6) G. Bianconi, Multilayer networks:structure and function (Oxford University Press, 2018).
  * (7) S. V. Buldyrev, R. Parshani, G. Paul, H. E. Stanley, and S.Havlin, Nature (London) 464, 1025 (2010).
  * (8) S. Boccaletti, G. Bianconi, R. Criado, C. I. del Genio, J. Gómez-Gardeñes, M. Romance, I. Sendiña-Nadal, Z. Wang, and M. Zanin, Phys. Rep. 544, 1 (2014).
  * (9) M. Kivelä, A. Arenas, M. Barthelemy, J. P. Gleeson, Y. Moreno, and M. A. Porter, J. Complex Netw. 2, 203 (2014).
  * (10) K.-M. Lee, B. Min, and K.-I. Goh, Eur. Phys. Jour. B 88, 48 (2015).
  * (11) M. Szell, R. Lambiotte, and S. Thurner, Proc. Natl. Acad. Sci. U.S.A. 107, 13636 (2010).
  * (12) G. Menichetti, D. Remondini, P. Panzarasa, R. J. Mondragón, and G. Bianconi, PloS one 9, e97857 (2014).
  * (13) A. Cardillo, J. Gómez-Gardeñes, M. Zanin, M. Romance, D. Papo, F. del Pozo, and S. Boccaletti, Sci. Rep. 3, 1344 (2013).
  * (14) E. Bullmore and O. Sporns, Nat. Rev. Neurosci. 10, 186 (2009).
  * (15) S. D. S. Reis, Y. Hu, A. Babino, J. S. Andrade Jr., S. Canals, M. Sigman, and H. A. Makse, Nat. Phys. 10, 762 (2014).
  * (16) G. Bianconi, Phys. Rev. E 87, 062806 (2013).
  * (17) B. Min, S. D. Yi, K.-M. Lee, and K.-I. Goh, Phys. Rev. E 89, 042811 (2014).
  * (18) V. Nicosia, V. Latora, Phys. Rev. E 92, 032805 (2015).
  * (19) H. Wu, R. G. James, and R. M. D’Souza Jour. Comp. Net. 8, 2 (2020).
  * (20) G. J. Baxter, S. N. Dorogovtsev, A. V. Goltsev, and J. F. F. Mendes, Phys. Rev. Lett. 109, 248701 (2012).
  * (21) D. Cellai, E. Lopez, J. Zhou, J. P. Gleeson, and G. Bianconi, Phys. Rev. E 88, 052811 (2013).
  * (22) S. Hwang, S. Choi, D. Lee, and B. Kahng, Phys. Rev. E 91, 022814 (2015).
  * (23) G. Bianconi, S. N. Dorogovtsev, and J. F. F. Mendes, Phys. Rev. E 91, 012804 (2015).
  * (24) D. Cellai, S. N. Dorogovtsev, and G. Bianconi Phys. Rev. E 94, 032301 (2016).
  * (25) S. Gómez, A. Díaz-Guilera, J. Gómez-Gardeñes, C. J. Pérez-Vicente, Y. Moreno, and A. Arenas, Phys. Rev. Lett. 110, 028701 (2013).
  * (26) M. De Domenico, A. Solé-Ribalta, S. Gómez, and A. Arenas, Proc. Natl. Acad. Sci. U.S.A. 111, 8351 (2014).
  * (27) A. Saumell-Mendiola, M. A. Serrano, and M. Boguna, Phys. Rev. E 86, 026106 (2012).
  * (28) E. Cozzo, R. A. Banos, S. Meloni, and Y. Moreno, Phys. Rev. E 88, 050801(R) (2013).
  * (29) C. Granell, S. Gómez, and A. Arenas, Phys. Rev. Lett. 111, 128701 (2013).
  * (30) Z. Wang, L. Wang, A. Szolnoki, and M. Perc, Eur. Phys. Jour. B 88, 124 (2015).
  * (31) M. Perc, J. Gómez-Gardeñes, A. Szolnoki, L. M. Floría, and Y. Moreno. J. R. Soc. Interface 10, 20120997 (2013).
  * (32) N. Masuda, Phys. Rev. E 90, 012802 (2014).
  * (33) M. Diakonova, M. San Miguel, and V. M. Eguíluz, Phys. Rev. E 89, 062818 (2014).
  * (34) M. Diakonova, V. Nicosia, V. Latora, and M. San Miguel, New J. Phys. 18, 023010 (2016).
  * (35) A. Chmiel and K. Sznajd-Weron, Phys. Rev. E 92, 052812 (2015).
  * (36) A. Chmiel, J. Sienkiewicz, A. Fronczak, and P. Fronczak, Entropy 22 (9), 1018 (2020).
  * (37) A. Halu, K. Zhao, A. Baronchelli, and G. Bianconi, Europhys. Lett. 102, 16002 (2013).
  * (38) F. Battiston, A. Cairoli, V. Nicosia, A. Baule, and V. Latora, Physica D 323-324, 12 (2016).
  * (39) J. Ashkin and E. Teller, Phys. Rev. 64, 178 (1943).
  * (40) S. Jang, J.S. Lee, S. Hwang, and B. Kahng, Phys. Rev. E 92, 022110 (2015).
  * (41) R.V. Ditzian and L. P. Kadanoff, J. Phys. A 12, L229 (1979).
  * (42) F. Iglói and L. Turban, Phys. Rev. E 66, 036140 (2002).
  * (43) S. N. Dorogovtsev, A. V. Goltsev, and J. F. F. Mendes, Eur. Phys. J. B 38, 177 (2004).
  * (44) S. H. Lee, M. Ha, H. Jeong, J. D. Noh and H. Park, Phys. Rev. E 80, 051127 (2009).
  * (45) F. Battiston, G. Cencetti, I. Iacopini, V. Latora, M. Lucas, A. Patania, J. G. Young, and G. Petri, Phys. Rep. 874, 1-92 (2020).
  * (46) G. Bianconi, Higher-order networks: An introduction to simplicial complexes (Cambridge University Press, 2021).