# Multilayer social reinforcement induces bistability on multiplex networks

Longzhao Liu1,2,4,6, Xin Wang1,3,6, Shaoting Tang1,6,7, Hongwei Zheng5,7 and
Zhiming Zheng1,5,6 1 LMIB, NLSDE, BDBC, School of Mathematical Sciences,
Beihang University, Beijing 100191, China  2 ShenYuan Honor School, Beihang
University, Beijing 100191, China  3 Department of Mathematics, Dartmouth
College, Hanover, NH 03755, USA  4 Northwestern Institute on Complex Systems,
Northwestern University, Evanston, IL 60208, USA  5 Institute of Artificial
Intelligence and Blockchain, Guangzhou university, Guangdong province 510006,
China  6 PengCheng Laboratory, Shenzhen, 518055 , China  7 Author to whom any
correspondence should be addressed<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Social reinforcement mechanism, which characterizes the promoting effects when
exposing to multiple sources in social contagion process, is ubiquitous in
information-technology ecosystem and has aroused great attention in recent
years. While the impacts of social reinforcement on single-layer networks are
well-documented, extension to multilayer networks is needed to study how
reinforcement from different social circles influences the spreading dynamics.
To this end, we incorporate multilayer social reinforcement into ignorant-
spreader-ignorant (SIS) model on multiplex networks. Our theoretical analysis
combines pairwise method and mean-field theory and agrees well with large-
scale simulations. Surprisingly, we find this complex social contagion
mechanism triggers the emergence of bistability phenomena, where extinction
and outbreak states coexist. In particular, the hysteresis loop of stationary
prevalence occurs in this bistable region, explaining why the fight towards
the spread of rumors is protracted and difficult in modern society. Further,
we show that the final state of bistable regions depends on the initial
density of adopters, the critical value of which decreases as the contagion
transmissibility or the multilayer reinforcement increases. In particular, we
highlight two possible conditions for the outbreak of social contagion: to
possess large contagion transmissibility, or to possess large initial density
of adopters with strong multilayer reinforcement. Our results unveil the non-
negligible power of social reinforcement on multiplex networks, which sheds
lights on designing efficient strategies in spreading behaviors such as
marketing and promoting innovations.

Keywords: complex social contagion, multiplex networks, multilayer social
reinforcement, bistability, hysteresis loop.

##  1 Introduction

Social contagion describes a variety of behavioral imitations caused by social
influence and is particularly ubiquitous in the digital age [1, 2, 3]. In
order to predict and control the collective contagion phenomena on large-scale
social networks, such as the spread of rumors, social norms and online
behaviors, scientists have made great efforts on understanding the underlying
dynamical mechanisms [4, 5, 6, 7]. Early studies argued that the dynamical
processes of social contagion are similar to disease spreading in the sense
that a simple contact with a single infected individual can trigger diffusion
[8, 9, 10]. Therefore, many epidemic-like models are proposed, which are
called simple contagion models. These models provide profound insights into
many physical phenomena on social networks, such as the cascading process of
information diffusion, the co-contagion dynamics on multiplex networks, and
etc [11, 12].

However, simple contagion model can not deal with the spread of more complex
social behaviors, especially when the behaviors are risky, costly or
polarized, ranging from the spread of public health behaviors (e.g., vaccine,
vaping, diet) to social movement [13, 14, 15]. For instance, a well-known
complex mechanism of social contagion is called the social reinforcement. It
corresponds to the fact that exposure to multiple sources with a same stimuli
would give individuals more confidence to participate than multiple exposures
to the same source, leading to a significant promotion on social transmission
[16, 17]. As McAdam and Paulsen concluded [18], “the fact that we are embedded
in many relationships means that any major decision we are contemplating will
likely be mediated by a significant subset of those relationships”.
Nevertheless, this empirically confirmed phenomenon is fundamentally different
from the view of simple contagion model [19, 20]. Hence, more complex
mechanisms were further incorporated into dynamical models. These complex
contagion models lead to abundant intriguing findings, some of which even
conflict with conclusions of simple contagion models [21, 22]. For example,
threshold models assumed that social contagion only happens when the number
(or fraction) of exposures to multiple sources exceeds a given threshold and
suggested the dependence of diffusion results on initial density of adopters
[23]. In addition, Centola et al. proved that the clustering network structure
suppressed simple contagion while surprisingly facilitated the spread of
behaviors that require social reinforcement [24]. Another framework modified
epidemic-like models by increasing transmission rate when exposing to multiple
sources, and showed the promoting effect of social reinforcement on contagion
processes [25, 26]. Another well-known complex mechanism is caused by the
interplay between spreading dynamics and network topology. For example, Gross
et al. studied contagion process on adaptive networks where the connection
between nodes relies on their states, and showed that this dynamics-topology
interaction could lead to the emergence of bistable region where healthy and
endemic states co-exist [27]. Besides, Iacopini et al. noted that the
existence of high-order structures embedded on single-layer networks, such as
full triangles, could also induce bistable phenomena [28].

Moreover, complex mechanisms depicting interactions between multiple contagion
processes on networked systems, such as the spreading of interacting diseases,
were widely investigated [29]. Chen et al. explored cooperative contagion
processes where individuals infected by one disease are more susceptible to
the other, and showed that large degree of cooperation could lead to the
occurrence of abrupt phase transition and bistability [30]. Soriano-Paños et
al. described the interaction between cooperative or competitive diseases by
increasing or decreasing the susceptibility of individuals infected by one of
the diseases, and also found bistability caused by cooperativity between
diseases [31]. Sanz et al. proposed a general framework analyzing simultaneous
spreading of two interacting diseases and derived complex phase diagrams [32].
Pinotti et al. explored a three-player pathogen system where competition and
cooperation coexist, and showed that the presence of cooperative pathogen
could lead to some intriguing phenomena including non-monotonic boundaries
separating phase diagram and bistability [33]. There were also many works
extending the framework from single-layer networks to multilayer networks
[34]. For example, Wu et al. studied discrete-time Markov-Chain model
depicting the spreading of two diseases in multiplex networks and derived
epidemic thresholds [35].

Recently, owing to the development of various social medias, many social
behaviors spread in multiple social circles rather than in a single social
circle, which are naturally modeled as multilayer networks [36, 37, 38, 39].
Each layer represents a single social circle which could be an online social
platform or a network of offline relationships, such as friendships and
colleagues. Similar to social reinforcement in a single-layer network, there
exists significant difference between exposure to multiple social circles
(multilayer reinforcement) and multiple exposures to the same social circle
(intra-layer reinforcement)[40]. Specifically, individuals would be more
convinced and have higher possibility for diffusion when receiving the same
informative stimulus from different social circles. However, previous studies
mainly concentrated on simple interactions between nodes and its replicas in
different layers [41, 42]. The detailed impacts of multilayer reinforcement
mechanism that engineers complex social contagion remain largely unknown.

To fill this gap, here we propose a theoretical framework that incorporates
multilayer reinforcement into ignorant-spreader-ignorant (SIS) model to study
the spreading dynamics on multiplex networks. We find this complex social
contagion mechanism not only expands dissemination, but also results in the
emergence of bistable region, within which extinction and outbreak states
coexist. Within bistable region, the final state depends on the initial
density of adopters. We also detailedly discuss the hysteresis loop and the
unstable equilibrium manifolds occurred in the bistable region. Through
observing phase diagrams, we highlight two conditions for the outbreak of
social contagion: (i) to possess large transmissibility; (ii) to possess large
initial density of adopters with strong multilayer reinforcement. As the
second condition cannot be obtained by simple contagion model and is less
intuitive, our results stress the unneglectable role of complex social
contagion and are in line with the findings of previous experimental studies
that critical masses are required for establishing collective behaviors such
as social changes [43]. In addition, similar phenomena are also observed on
finite-size heterogeneous multiplex networks. Our findings provide valuable
insights toward dynamical evolutions of complex social contagion on multiplex
networks, which are of vital importance for understanding collective online
behaviors in the era of social media. In particular, the emerging bistable
phenomenon on multiplex networks indicates that the information (truth, ideas,
advertisements, etc.) can stride across the trap of extinction and become
widespread by selecting a large proportion of initial spreaders on multiple
platforms, which implies effective strategies for controlling rumors,
promoting innovations and marketing. [44, 45, 46].

##  2 Model

Consider an undirected multiplex network with two layers, denoted as layer 1
and layer 2. Each layer stands for a social circle, which is an online social
platform or an offline relationship network composed of interacting
individuals such as friends, families and colleges. Both layers have the same
nodes and interlayer edges only connect entities with their replicas, as shown
in figure 1(a).

Here we stress that our main purpose is to explore the influence of complex
social contagion aroused by multiple social circles, i.e., multilayer
reinforcement mechanism. Therefore, we adopt ignorant-spreader-ignorant (SIS)
model rather than complex contagion model to characterize the intralayer
spreading process, which provides better analytical insights while at expense
of being less realistic. In this situation, spreader ($S$) represents
individuals who adopt norms (information, cognition, attitudes, behaviors and
so on) and are willing to spread, while ignorant ($I$) stands for nodes who do
not adopt norms or have no motivation to spread, corresponding to infected and
susceptible state in epidemiology, respectively. In each layer $k$ $(k=1,2)$,
spreader has a probability $\lambda_{k}$ to spread norms to its ignorant
neighbors and becomes ignorant with probability $\mu_{k}$, as shown in figure
1(b).

Beyond intralayer contagion processes, norms also diffuse across layers, i.e.,
the interactions between layers. First, complying with previous studies, we
utilize interlayer contagion processes where ignorant has a probability $p$ to
become a spreader once its counterpart is in $S$ state, which depicts the
interaction between agents and their replicas, as shown in figure 1(c) [47].
Then, we introduce multilayer reinforcement mechanism to mimic the fact that
exposure to multiple social circles is more convincing than multiple exposures
to a single social circle. At each time step, if an individual can receive
norms from both social circles, i.e. he/she has at least one spreader neighbor
in each layer, then with probability $\gamma$ the individual would randomly
choose one layer to spread the norms (figure 1(d)). The parameter $\gamma$
reflects the strength of multilayer reinforcement.

In summary, our model is composed of three processes: intralayer contagion,
interlayer contagion and multilayer reinforcement, which is as follows:

  1. 1.

Intralayer contagion. At each time step, in layer $k$, $S$ makes its neighbors
in $I$ state become $S$ with probability $\lambda_{k}$. Meanwhile, $S$ turns
into $I$ with probability $\mu_{k}$.

  2. 2.

Interlayer contagion. At each time step, $I$ becomes $S$ with probability $p$
if its replica is a spreader.

  3. 3.

Multilayer reinforcement. At each time step, a node has a probability $\gamma$
to randomly choose a layer and become a spreader if there exists at least one
spreader in its neighborhood of each layer.

Figure 1: Schematic of modeling framework. (a) Network structure. The
underlying network is composed of two layers with the same nodes, where
interlayer edges only connect entities and their replicas. Node’s state is
represented by different colors: spreader (red), ignorant (black). (b)-(d)
Dynamical model. Our model includes three dynamical processes, whose
propagation path is highlighted by red arrows in each subfigure. (b)
Intralayer contagion. In each layer, an ignorant has a probability to be
infected by its spreader neighbors. (c) Interlayer contagion. An ignorant may
change its state to a spreader if its replica is a spreader. (d) Multilayer
reinforcement. When node $i$ is connected to at least one spreader in each
layer and is able to receive stimuli from both platforms at the same time, $i$
randomly chooses a layer and becomes a spreader with a certain probability.

##  3 Theoretical framework

In this section, we explore dynamical equations of our model on homogeneous
multiplex systems, composed of two Erdös-Rényi (ER) networks with average
degree $\langle k_{1}\rangle$ and $\langle k_{2}\rangle$. The major challenge
is the dynamical correlation between two layers caused by interlayer contagion
and multilayer reinforcement. Specifically, the counterpart of a spreader is
more likely to be a spreader than the counterpart of an ignorant. Classic
Mean-field approximation does not consider the dynamical correlation, which
leads to the deviation of theoretical predictions from Monte Carlo simulations
(see A). To overcome the difficulty, we conduct theoretical analysis by
combining pairwise method and mean-field theory.

Here, we define $\rho_{X,Y}(t)$ as the probability that the individual is $X$
state in layer 1 and $Y$ state in layer 2. Clearly, each individual has four
possible states and we have
$\rho_{S,S}(t)+\rho_{S,I}(t)+\rho_{I,S}(t)+\rho_{I,I}(t)=1$. To begin, we
explore intralayer contagion process. Let $\theta_{k}(t)$ represent the
probability that a node in layer $k$ is not convinced by its neighbors of the
same layer. Utilizing the mean-field theory, $\theta_{k}(t)$ can be
approximated by

$\eqalign{\theta_{1}(t)&=(1-\lambda_{1}(\rho_{S,S}(t)+\rho_{S,I}(t)))^{\langle
k_{1}\rangle}\\\
\theta_{2}(t)&=(1-\lambda_{2}(\rho_{S,S}(t)+\rho_{I,S}(t)))^{\langle
k_{2}\rangle}.}$ (1)

Then, we examine the influence of multilayer reinforcement. Denote
$\delta_{k}(t)$ as the probability that the multilayer reinforcement does not
make ignorant turn into a spreader in layer $k$. We have
$\delta_{1}(t)=\delta_{2}(t)=\delta(t)$, which can be approximated by

$\eqalign{\delta(t)=1-0.5&*\gamma\\\
&*(1-(1-\rho_{S,S}(t)-\rho_{S,I}(t))^{\langle k_{1}\rangle})\\\
&*(1-(1-\rho_{S,S}(t)-\rho_{I,S}(t))^{\langle k_{2}\rangle}),}$ (2)

where $(1-(1-\rho_{S,S}(t)-\rho_{S,I}(t))^{\langle k_{1}\rangle})$ and
$(1-(1-\rho_{S,S}(t)-\rho_{I,S}(t))^{\langle k_{2}\rangle})$ account for the
probabilities that there exists at least one spreader neighbor in layer 1 and
in layer 2, respectively.

Furthermore, let $g^{I}_{k}(t)$ and $g_{k}^{S}(t)$ represent the probabilities
that ignorant in layer $k$ is not convinced if its counterpart is in ignorant
and spreader state, respectively. The probabilities describe the integrated
influence of intralayer contagion, interlayer contagion and multilayer
reinforcement. Because the three processes are approximately independent in
our model, the probabilities read

$\eqalign{g^{I}_{k}(t)&=\theta_{k}(t)*\delta(t)\\\
g^{S}_{k}(t)&=\theta_{k}(t)*\delta(t)*(1-p).}$ (3)

Finally, the temporal evolutionary equations of $\rho_{S,S}(t)$ can be written
as

$\eqalign{\frac{d\rho_{S,S}}{dt}&=-\rho_{S,S}\\{1-(1-\mu_{1})(1-\mu_{2})\\}+\rho_{S,I}(1-\mu_{1})(1-g^{S}_{2})\\\
&+\rho_{I,S}(1-\mu_{2})(1-g^{S}_{1})+\rho_{I,I}(1-g^{I}_{1})(1-g^{I}_{2}),\\\
}$ (4)

where the first term represents the outflow from the $S-S$ class and the last
three terms stand for transition from the other states to $S-S$ state.
Similarly, we can derive the dynamical evolutions of $\rho_{S,I}(t)$,
$\rho_{I,S}(t)$ and $\rho_{I,I}(t)$. For the sake of theoretical analysis and
readability, $\rho_{I,I}(t)$ is automatically substituted for
$(1-\rho_{S,S}(t)-\rho_{S,I}(t)-\rho_{I,S}(t))$. Thus, the evolutionary
equations of our model read as follows:

$\eqalign{\frac{d\rho_{S,S}}{dt}&=-\rho_{S,S}\\{1-(1-\mu_{1})(1-\mu_{2})\\}+\rho_{S,I}(1-\mu_{1})(1-g^{S}_{2})\\\
&+\rho_{I,S}(1-\mu_{2})(1-g^{S}_{1})+(1-\rho_{S,S}-\rho_{S,I}-\rho_{I,S})(1-g^{I}_{1})(1-g^{I}_{2})\\\
\frac{d\rho_{S,I}}{dt}&=\rho_{S,S}(1-\mu_{1})\mu_{2}-\rho_{S,I}\\{1-(1-\mu_{1})g^{S}_{2}\\}\\\
&+\rho_{I,S}(1-g^{S}_{1})\mu_{2}+(1-\rho_{S,S}-\rho_{S,I}-\rho_{I,S})(1-g^{I}_{1})g^{I}_{2}\\\
\frac{d\rho_{I,S}}{dt}&=\rho_{S,S}(1-\mu_{2})\mu_{1}+\rho_{S,I}(1-g^{S}_{2})\mu_{1}\\\
&-\rho_{I,S}\\{1-(1-\mu_{2})g^{S}_{1}\\}+(1-\rho_{S,S}-\rho_{S,I}-\rho_{I,S})g^{I}_{1}(1-g^{I}_{2}).\\\
}$ (5)

Equation (5) can be simply written as

$\frac{d\boldsymbol{\rho}}{dt}=\boldsymbol{f}(\boldsymbol{\rho}),$ (6)

where $\boldsymbol{\rho}(t)=(\rho_{S,S}(t),\rho_{S,I}(t),\rho_{I,S}(t))$. This
indicates that the dynamical system is 3-dimension autonomous. Thus, the
stability of the fixed points directly determines evolutionary results of the
system. We define $\boldsymbol{\rho}_{f}$ as the fixed point of equation (5),
which satisfies

$\boldsymbol{f}(\boldsymbol{\rho}_{f})=\boldsymbol{0}.\\\ $ (7)

In particular, $\boldsymbol{\rho}_{f}$ represents the final state of the
system if and only if the fixed point is stable, i.e.,

$\Lambda_{max}(\frac{d\boldsymbol{f}}{d\boldsymbol{\rho}}|_{\boldsymbol{\rho}=\boldsymbol{\rho}_{f}})<0,\\\
$ (8)

where $\Lambda_{max}(J)$ is the largest eigenvalue of matrix $J$.

We notice that there is at least one fixed point
($\boldsymbol{\rho}_{f}=\boldsymbol{0}$), which represents that all
individuals are in $I-I$ state. When the initial density of spreaders is
small, norms go extinct if and only if $\boldsymbol{\rho}_{f}=\boldsymbol{0}$
is stable, i.e.,

$\Lambda_{max}(\frac{d\boldsymbol{f}}{d\boldsymbol{\rho}}|_{\boldsymbol{\rho}=\boldsymbol{0}})<0.\\\
$ (9)

##  4 Results

###  4.1 Complex social contagion on homogeneous multiplex networks

Figure 2: Emergence of bistable phenomena. Prevalence curves are shown against
(a)-(b) intralayer transmissibility and (c)-(d) interlayer transmissibility
for different values of multilayer reinforcement and initial density of
adopters. In each figure, $\gamma=0$ (red), $\gamma=0.1$ (yellow),
$\gamma=0.3$ (blue) and $\rho_{0}=0.03$ (square), $\rho_{0}=0.2$ (diamond).
Simulation results are averaged over 30 independent runs, whose standard
errors are shown in subplots. In most cases, errors are within 0.01. However,
errors can exceed 0.01 when the system is close to the critical points, as the
stochastic fluctuation of simulations is relatively large around the
thresholds. Theoretical predictions solved by equation (5) are shown by solid
lines. Note that for $\gamma=0.3$ in (a) and $\gamma=0.1$ or $\gamma=0.3$ in
(c), there emerges bistable region where outbreak and extinction states co-
exist. Parameters: $\mu_{1}=0.6$, $\mu_{2}=0.6$. In addition, (a)
$\lambda_{2}=0.02$, $p=0.1$, (b) $\lambda_{2}=0.1$, $p=0.1$, (c)
$\lambda_{1}=0.06$, $\lambda_{2}=0.04$, and (d) $\lambda_{1}=0.14$,
$\lambda_{2}=0.06$.

In this section, we explore how our model behaves on homogeneous multiplex
networks, especially the influence of multilayer reinforcement on the
prevalence of norms and the critical properties that separate outbreak and
extinction. We start from homogeneous multiplex networks composed of two
Erdös-Rényi (ER) graphs with $N=10000$ nodes, the average degree of which are
$\langle k_{1}\rangle=6$ and $\langle k_{2}\rangle=8$, respectively.
Initially, we randomly set a certain fraction of population as adopters,
denoted by $\rho_{0}$, who are spreaders in both layers[48]. Here we run
300-step simulations to ensure that the system has reached the dynamic
equilibrium.

Figure 3: Hysteresis phenomena. Shown are the hysteresis loops of the
stationary prevalence with respect to changes in (a) intralayer
transmissibility and (b) interlayer transmissibility. Two different
evolutionary routes arise: prevalence increases from zero to positive values
along the ascending path (blue squares) with $\lambda_{1}$ or $p$ grows, while
recovers to zero along the descending path (red triangles) as transmissibility
decreases. Theoretical predictions (solid lines) are in good agreement with
simulations in most situations. The relatively large stochastic fluctuation of
simulations only occurs around the thresholds. Parameters: (a)
$\lambda_{2}=0.02$,$p=0.1$,$\gamma=0.3$,$\mu_{1}=\mu_{2}=0.6$; (b)
$\lambda_{1}=0.06$, $\lambda_{2}=0.04$, $\gamma=0.12$, $\mu_{1}=\mu_{2}=0.6$.
Figure 4: Effect of initial density of adopters. (a)-(d) Time evolutions for
the prevalence are presented under different combinations of transmission
parameters. In each figure, a single curve corresponds to one value of
different initial densities of adopters. (b) We find that when $\gamma$ is
relatively large, the initial condition determines the final state of bistable
region, either outbreak or extinction. Parameters: $\mu_{1}=0.6$,
$\mu_{2}=0.6$, $p=0.1$. In addition, (a) $\lambda_{1}=0.05$,
$\lambda_{2}=0.02$, $\gamma=0.1$, (b) $\lambda_{1}=0.05$, $\lambda_{2}=0.02$,
$\gamma=0.3$, (c) $\lambda_{1}=0.2$, $\lambda_{2}=0.02$, $\gamma=0.1$, and (d)
$\lambda_{1}=0.05$, $\lambda_{2}=0.1$, $\gamma=0.1$.

To begin with, we show prevalence curves as a function of intralayer
transmissibility (figure 2(a) and figure 2(b)) and interlayer transmissibility
(figure 2(c) and figure 2(d)) under different combinations of multilayer
reinforcement and initial density of adopters. All subfigures verify the
intuition that multilayer reinforcement promotes social contagion and show
that our theoretical predictions agree well with simulation results. The case
$\gamma=0$ (the red curves) is equivalent to SIS model on multiplex networks ,
which displays continuous phase transition [41]. Nevertheless, the case
$\gamma=0.3$ (the blue curves) shows large differences in dependence of
prevalence on intralayer transmissibility (figure 2(a)) and interlayer
transmissibility(figure 2(c)). Phase transition appears at lower value of
$\lambda_{1}$ or $p$, and becomes discontinuous. Another interesting
phenomenon is the emergence of bistable region, where outbreak and extinction
states coexist. Specifically, for $\lambda_{1}\in(0.03,0.08)$ in figure 2(a)
and $p\in(0,0.12)$ in figure 2(c), norms outbreak if $\rho_{0}=0.2$, while go
extinct if $\rho_{0}=0.03$. It indicates that the final state in bistable
region might depend on the initial density of adopters. Owing to the
spontaneous recovery mechanism ($\mu_{1}$ and $\mu_{2}$), the system would
finally reach a dynamic equilibrium. We conduct additional simulations under
different $\mu_{1}$ and $\mu_{2}$, and obtain similar bistable phenomena and
discontinuous phase transitions even under weak multilayer reinforcement when
$\mu_{1}$ and $\mu_{2}$ are small (see figure 11 in B).

In addition, figure 3 shows the emergence of hysteresis loops of stationary
prevalence with respect to changes in intralayer transmissibility and
interlayer transmissibility. Here two evolutionary routes can be observed. One
is the ascending path (blue squares) describing the stationary prevalence
under small initial density of adopters, which corresponds to diffusion
processes from few early adopters to global dissemination. For the ascending
path, prevalence first remains zero and then rapidly grows to a high level at
a large threshold as $\lambda_{1}$ or $p$ increases. We call the threshold as
the diffusion threshold. The other one is the descending path (red triangles)
which depicts the final results under large initial density of adopters. It
corresponds to the eradication processes that norms transform from a large
initial prevalence to a small one. For the descending path, as $\lambda_{1}$
or $p$ decreases, prevalence remains a high level and would not recover to
zero until it mitigates less than a small threshold. We designate the
threshold as the eradication threshold. These two paths form hysteresis loops
of stationary prevalence, where the diffusion threshold is much larger than
the eradication threshold. This indicates that more efforts of reducing
transmissibility are required to make the prevalence recover to zero. Our
result explains why it is so difficult to eliminate the rumors and misbeliefs
from social networks, especially in the modern media environment where
different online social platforms interact. Besides, noting the dependence of
the system’s final state on its history, we could alter initial density of
adopters to realize the mutual transformation of the two evolutionary routes.

To give a more intuitive illustration about how initial density of adopters
affects the final state, in figure 4, we further present time evolutions of
prevalence for four different combinations of transmission parameters. In each
subfigure, a single realization indicates the temporal evolution under a
certain value of initial condition ($\rho_{0}$), which ranges from 0 to 1 and
is represented by different colors. The initial conditions show no influence
on the prevalence of norms when $\gamma$ is small ($\gamma=0.1$): either go
extinct (figure 4(a)) or outbreak (figure 4(c) and figure 4(d)) no matter what
$\rho_{0}$ is. However, figure 4(b) shows completely different phenomena when
$\gamma$ is relatively large ($\gamma=0.3$). There appears a threshold of
initial density of adopters ($\rho_{0}=0.09$), below which the prevalence
vanishes, while above which the prevalence converges to about 0.3. This
indicates that initial density of adopters plays a main role in determining
the final state of bistable region. This insight is consistent with empirical
studies that critical masses are necessary for initiating social changes, and
can not be observed in simple contagion models on multiplex networks [43, 49].

Figure 5: Conditions of phase transition for bistability. We provide phase
diagrams (a) for the joint effects of multilayer reinforcement and intralayer
transmissibility and (b) for the joint effects of multilayer reinforcement and
interlayer transmissibility. In each subfigure, the phase plane is composed of
three regions: outbreak (region 1), bistable where outbreak and extinction
states coexist (region 2), extinction (region 3). The separatrixes $L_{1}$ and
$L_{2}$ represent the critical parameters when $\rho_{0}$ is small and when
$\rho_{0}=1$, and are numerically solved by equation (9) and equation (5).
Parameters: (a) $\lambda_{2}=0.02$, $\mu_{1}=0.6$, $\mu_{2}=0.6$, $p=0.1$ and
(b) $\lambda_{1}=0.06$, $\lambda_{2}=0.04$, $\mu_{1}=0.6$, $\mu_{2}=0.6$.

Note that the emergence of bistable region is an intriguing physical
phenomenon arising from multilayer reinforcement. Here, we further explore the
detailed conditions of phase transition for bistability. In figure 5, we
present phase diagram under different combinations of multilayer reinforcement
and intralayer transmissibility (figure 5(a)) or combinations of multilayer
reinforcement and interlayer transmissibility. (figure 5(b)). In all
subfigures, the phase plane is divided into three regions by two separatrixes
($L_{1}$ and $L_{2}$): outbreak, bistable, extinction. The outbreak state
means that the prevalence is positive as long as $\rho_{0}>0$, and the
extinction state represents that the prevalence vanishes no matter what
$\rho_{0}$ is. The bistable state means that norms outbreak for large
$\rho_{0}$ while go extinct for small $\rho_{0}$. Thus, two separatrixes
$L_{1}$ and $L_{2}$ are critical values of the parameters under very small
$\rho_{0}$ and $\rho_{0}=1$, which correspond to the diffusion threshold and
the eradication threshold of transmissibility, respectively. $L_{1}$ and
$L_{2}$ can be numerically solved by equation (9) and equation (5),
respectively. We find that the solution of equation (9) is uncorrelated with
multilayer reinforcement ($\gamma$), which is directly reflected by the
parallel relationship between $L_{1}$ and y-axis. Meanwhile, $L_{2}$ is
determined by the joint effects of intralayer transmissibility, interlayer
transmissibility and multilayer reinforcement. To sum up, region 1 and region
2 highlight two opportunities for outbreak, which are: (i) to own large
transmissibilities and (ii) to own large initial density of adopters with
strong multilayer reinforcement. These conclusions can be directly applied to
many scenarios such as promoting marketing, designing interventions for rumor
spreading and establishing new social norms [50, 43]. We also discuss the
impacts of system size on the critical value of multilayer reinforcement that
leads to bistability, the diffusion threshold of transmissibility as well as
the eradication threshold of transmissibility in C. We find that all the
critical values first decrease as $N$ grows when the system is small, and then
become stable when $N$ is large ($N>8000$), which also proves the reliability
of our large-scale simulations.

Figure 6: Emergence of unstable manifold in bistable region. The prevalence of
multiple equilibriums solved by equation (5) is shown as a function of (a)
intralayer transmissibility and (b) interlayer transmissibility for different
multilayer reinforcement. For $\gamma=0.3$ or $\gamma=0.2$ where bistable
region exists, the unstable equilibrium manifold (dashed lines) emerges
between the stable equilibrium manifolds (solid lines). Vertical line
represents the threshold when the system has no multilayer reinforcement.
Parameters: (a) $\lambda_{2}=0.02$, $p=0.1$, $\mu_{1}=\mu_{2}=0.6$; (b)
$\lambda_{1}=0.06$, $\lambda_{2}=0.04$, $\mu_{1}=\mu_{2}=0.6$.

Furthermore, we show the existence of unstable equilibrium manifold in
bistable region in Figure 6. By solving equation (5), we present the
prevalence of multiple equilibriums with respect to changes in intralayer or
interlayer transmissibility. Again, our results illustrate that strong
multilayer reinforcement can result in the emergence of discontinuous phase
transition and bistability. Moreover, we find that the unstable manifold
appears between two stable manifolds in bistable region and the prevalence of
unstable equilibriums decreases as $\lambda_{1}$ or $p$ grows.

Figure 7: Critical density of initial adopters ($\rho_{c}$) which separates
the bistable region into outbreak and extinction. The boxplots present
$\rho_{c}$ as a function of (a) intralayer transmissibility, (b) interlayer
transmissibility and (c) multilayer reinforcement. Each box contains 30
independent runs. Theoretical predictions are numerically solved by equation
(5) and are represented by yellow triangles. Parameters: $\mu_{1}=0.6$,
$\mu_{2}=0.6$. In addition, (a) $\lambda_{2}=0.02$, $p=0.1$, $\gamma=0.4$, (b)
$\lambda_{1}=0.06$, $\lambda_{2}=0.04$, $\gamma=0.3$, (c) $\lambda_{1}=0.06$,
$\lambda_{2}=0.04$, $p=0.1$.

The bistability phenomenon provides a profound insight that large initial
density of adopters could trigger the outbreak even at low values of
transmissibility. Here we explore the thresholds of initial conditions
($\rho_{c}$) under different circumstances, which describe at least how many
initial adopters are required for the outbreak of norms in bistable region. In
figure 7(a)-(c), we use boxplot to present critical value $\rho_{c}$ as a
function of intralayer transmissibility ($\lambda_{1}$), interlayer
transmissibility ($p$) and multilayer reinforcement ($\gamma$), respectively.
Each box contains 30 independent simulations. We find that the differences
among all simulation results in each box are within 0.02, which implies the
stability of $\rho_{c}$. Besides, the critical density of initial adopters
($\rho_{c}$) decreases as intralayer transmissibility, interlayer
transmissibility or multilayer reinforcement increases. Also, our theoretical
solutions well predict simulation results. In particular, we notice that all
critical values are small (lower than 0.12), which indicates that manipulating
initial condition is a powerful and effective way to control social contagions
under the existence of multilayer reinforcement.

###  4.2 Complex social contagion on heterogeneous multiplex networks

Figure 8: Complex social contagion on heterogeneous multiplex networks. The
prevalence is presented as a function of intralayer transmissibility for
different combinations of multilayer reinforcement and initial density of
adopters on (a) SF-SF multiplex networks, exponent parameters of which are
$\alpha_{1}=2.7$ and $\alpha_{2}=2.9$, and (b) WS-WS multiplex networks with
rewiring probability equal to 0.01. Simulation results are averaged over 30
independent runs, whose standard errors are shown in subplots. Parameters:
$\lambda_{2}=0.02$, $\mu_{1}=0.6$, $\mu_{2}=0.6$, $p=0.1$. In addition,
$\gamma=0$ (red), $\gamma=0.1$ (yellow), $\gamma=0.3$ (blue) and
$\rho_{0}=0.01$ (square), $\rho_{0}=0.3$ (diamond).

Real-world networks often have power-law degree distribution, i.e, $p_{k}\sim
k^{-\alpha}$, or small-world characteristics [51, 52]. Studies have revealed
that the network structures play an important role in dynamical evolutions
[53]. Thus, to mimic real situations, in this section, we examine how our
model behaves on multiplex systems composed of two scale-free (SF) networks or
two Watts-Strogatz small-world (WS) networks, respectively. All networks have
10000 nodes. The rewiring probability of the two WS networks are both 0.01 and
the exponent parameters of the two SF networks are $\alpha_{1}=2.7$ and
$\alpha_{2}=2.9$.

Figure 8 presents the prevalence as a function of intralayer transmissibility
for different values of multilayer reinforcement and initial density of
adopters in SF-SF multiplex networks (figure 8(a)) and WS-WS multiplex
networks (figure 8(b)). Results show similar phenomena to ER-ER multiplex
networks, including promoting effect of multilayer reinforcement and the
emergence of bistability when $\gamma$ is large ($\gamma=0.3$).

Figure 9: Comparison of different strategies for distributing initial seeds.
Critical initial density leading to outbreak is presented as a function of
intralayer transmissibility when random selection (blue) and target selection
(black) are utilized to distribute the initial seeds, respectively. Each box
is derived by 30 independent simulations on SF-SF multiplex networks.
Parameters: $\lambda_{2}=0.01$, $p=0.1$, $\gamma=0.3$, $\mu_{1}=\mu_{2}=0.6$.

Finally, we explore the influence of strategies for selecting initial seeds on
SF-SF multiplex networks. Two strategies are considered: random selection
(randomly selecting initial seeds) and target selection (prefer to select
those with higher sum of degree as initial seeds). Figure 9 presents the
critical density of initial adopters leading to outbreak with respect to
changes in intralayer transmissibility when random selection and target
selection are adopted, respectively. Results show that a small initial density
is enough to result in an outbreak on SF-SF multiplex networks for both
strategies. Specifically, the threshold for target selection is much smaller
than random selection.

##  5 Conclusions and Discussions

Complex social contagion describes complicated behavioral evolutions that
cannot be characterized by simple contact-contagion models, such as public
opinion formation on controversial events, the spread of conspiracies and the
establishment of new social norms [54, 55]. These risky or polarized
collective behaviors are ubiquitous on various social networks and are largely
determined by the high-order interactions regarding to the enhancement effects
of multiple exposures, i.e., the social reinforcement [56]. However, it
remains unclear how multilayer reinforcement which mimics complex contagion
process among multiple social circles influences the final diffusion results.

In this paper, we propose a modified SIS model which incorporates multilayer
reinforcement to describe spreading dynamics on multiplex networks. Firstly,
we examine how our model behaves on homogeneous multiplex networks. In
particular, we stress our efforts on the detailed impact of multilayer
reinforcement. A theoretical framework combining pairwise method and mean-
field theory is proposed and is verified by large-scale simulations.
Interestingly, we find that the multilayer reinforcement induces the emergence
of bistability, where extinction and outbreak states coexist. Furthermore, we
show that the final state of bistable region is determined by the initial
density of adopters whose threshold decreases as intralayer transmissibility,
interlayer transmissibility or multilayer reinforcement increases. We also
illustrate the hysteresis phenomena where two evolutionary routes occur, known
as an ascending path and a descending path, accounting for the phenomenon that
rumors are hard to dispel among online social platforms. The detailed
conditions of phase transition for bistability are further derived
analytically and are shown in phase diagrams. Within the bistable region, we
show the existence of unstable manifolds. Results highlight that the chance
for the outbreak of social contagion is either to possess large contagion
transmissibility or to own large initial density of adopters with strong
multilayer reinforcement. Our results are also valid on finite-size
heterogeneous multiplex networks.

Our basic yet powerful multilayer reinforcement mechanism reveals the dramatic
promoting impacts of complex social contagion on multiplex networks, which
indicates the possibility of facilitating spread outbreaks via controlling
initial spreaders. These insights are in line with the empirical studies that
only large initial density of adopters can successfully initiate social
changes [43]. While current results are performed on uncorrelated multiplex
networks, further study may extend the framework to partially overlapped
multiple networks where edges in different layers are correlated.

This work is supported by Program of National Natural Science Foundation of
China Grant No. 11871004, 11922102, and National Key Research and Development
Program of China Grant No. 2018AAA0101100.

##  Appendix A Comparison of Mean-field approximation and simulations

Figure 10: Comparison of Mean-field approximation and Pairwise approximation
as well as simulations. Prevalence is shown (a) against intralayer
transmissibility when interlayer interactions are weak and (b) against
interlayer transmissibility. Simulation results (blue squares) are averaged
over 30 runs, while mean-field approximation (black line) and pairwise
approximation (red line) are solved by equation (12) and equation (5),
respectively. Parameters: (a) $\lambda_{2}=0.02$, $p=0.1$, $\gamma=0.1$,
$\mu_{1}=\mu_{2}=0.6$; (b) $\lambda_{1}=0.06$, $\lambda_{2}=0.04$,
$\gamma=0.1$, $\mu_{1}=\mu_{2}=0.6$.

Here we consider Mean-field approach to describe the system. Let
$\rho_{X}^{j}(t),(j=1,2)$ represent the probability that state of node in
layer $j$ is $X$. Then we denote $\theta_{MF}^{j}(t)$ as the probability that
ignorant in layer $j$ is not convinced in intralayer contagion process, which
can be written as

$\theta_{MF}^{j}(t)=(1-\lambda_{j}\rho_{S}^{j}(t))^{\langle k_{j}\rangle}$
(10)

Next, we denote $\delta_{MF}(t)$ as the probability that ignorant in layer 1
or layer 2 is not affected by multilayer reinforcement. $\delta_{MF}(t)$ can
be calculated by

$\delta_{MF}(t)=1-0.5*\gamma*(1-(1-\rho_{S}^{1}(t))^{\langle
k_{1}\rangle})*(1-(1-\rho_{S}^{2}(t))^{\langle k_{2}\rangle})$ (11)

where $(1-(1-\rho_{S}^{j}(t))^{\langle k_{j}\rangle})$ accounts for the
probability that node in layer $j$ has at least one spreader neighbor.

Thus, the evolutionary equations of our dynamical model can be written as

$\eqalign{\frac{d\rho_{S}^{1}}{dt}&=-\mu_{1}\rho_{S}^{1}+(1-\rho_{S}^{1})*(1-\theta_{MF}^{1}(1-p\rho_{S}^{2})\delta_{MF})\\\
\frac{d\rho_{S}^{2}}{dt}&=-\mu_{2}\rho_{S}^{2}+(1-\rho_{S}^{2})*(1-\theta_{MF}^{2}(1-p\rho_{S}^{1})\delta_{MF})\\\
}$ (12)

In figure 10, we conduct comparisons of Mean-field approximation and Pairwise
approximation as well as simulations. When interlayer interactions are weak,
both approximations agree well with simulations (figure 10(a)). Nevertheless,
Mean-field approximation would deviate from simulation results when interlayer
interactions are strong, while Pairwise approximation is still valid (figure
10(b)). These indicate that considering dynamical correlations of two layers
is necessary for an accurate description of the system, especially under
strong interlayer interactions.

##  Appendix B complementary studies for figure 2

Figure 11: The emergence of bistability under different values of recovery. We
set (a) $\mu_{1}=\mu_{2}=0.1$ and (b) $\mu_{1}=\mu_{2}=0.2$, respectively. In
each subfigure, the prevalence is shown against intralayer transmissibility
under different combinations of multilayer reinforcement and initial density
of adopters. Simulation results are averaged over 30 independent runs, whose
standard errors are shown in subplots. Parameters: $\lambda_{2}=0.001$,
$p=0.1$. In addition, $\gamma=0$ (red), $\gamma=0.01$ (yellow), $\gamma=0.1$
(blue) and $\rho_{0}=0.002$ (square), $\rho_{0}=0.1$ (diamond).

Figure 11 presents the simulation results for different values of spontaneous
recovery. In all subfigures, we observe the emergence of bistable phenomena
and discontinuous phase transition when multilayer reinforcement is strong
($\gamma=0.1$). In particular, the case, $\gamma=0.01$, triggers a small
bistable region in figure 11(a), which indicates that even weak multilayer
reinforcement could induce bistability on multiplex networks when the value of
spontaneous recovery is small.

##  Appendix C Effect of system size

Figure 12: Effect of system size. (a) Critical value of multilayer
reinforcement ($\gamma_{c}$) inducing bistability is shown with respect to
changes in system size. As system size increases, $\gamma_{c}$ derived by
numerical simulations (blue squares) decreases and finally converges to our
theoretical predictions (red line). (b)-(c) Shown are the impacts of system
size on the diffusion threshold and the eradication threshold of intralayer
transmissibility ($\lambda_{1}$) and interlayer transmissibility ($p$),
respectively. In all subplots, simulation results are averaged over 30
independent runs and light blue/red region represents 95% confidence interval.
Parameters: $\langle k_{1}\rangle=6$, $\langle k_{2}\rangle=8$,
$\mu_{1}=\mu_{2}=0.6$. In addition, (a) $\lambda_{1}=0.06$,
$\lambda_{2}=0.02$, $p=0.1$; (b) $\lambda_{2}=0.02$, $p=0.1$, $\gamma=0.3$;
(c) $\lambda_{1}=0.06$, $\lambda_{2}=0.04$, $\gamma=0.12$.

In this section, we first explore how system size affects critical multilayer
reinforcement that induces bistability ($\gamma_{c}$) in figure 12(a). When
$N$ is small ($N<8000$), $\gamma_{c}$ decreases as system size grows. This
indicates that small system suppresses the emergence of bistability. When $N$
is large ($N>8000$), $\gamma_{c}$ is almost unchanged and is in good agreement
with our theoretical prediction. It proves the stability of $\gamma_{c}$ and
the validity of our theoretical framework when dealing with large-scale
systems.

As for the intralayer or interlayer transmissibility, there arise two
different thresholds: the diffusion threshold and the eradication threshold.
The diffusion threshold is the minimum transmissibility ensuring that few
early adopters of rumors could result in global diffusion. The eradication
threshold is the maximum transmissibility ensuring that the widespread rumors
could be finally eliminated. Here we further study how the diffusion threshold
and the eradication threshold change by the system sizes (see figure
12(b)-(c)). When $N$ is small ($N<8000$), both thresholds decrease as the
system size increases. When $N$ is large ($N>8000$), both thresholds are
almost unchanged. Meanwhile, our theoretical framework well estimates the
diffusion threshold as well as the eradication threshold of intralayer and
interlayer transmissibility when $N$ is large.

## References

  * [1] Axelrod R 1997 Journal of conflict resolution 41 203–226
  * [2] Del Vicario M, Bessi A, Zollo F, Petroni F, Scala A, Caldarelli G, Stanley H E and Quattrociocchi W 2016 Proceedings of the National Academy of Sciences 113 554–559
  * [3] Stewart A J, Mosleh M, Diakonova M, Arechar A A, Rand D G and Plotkin J B 2019 Nature 573 117–121
  * [4] Zhang Z K, Liu C, Zhan X X, Lu X, Zhang C X and Zhang Y C 2016 Physics Reports 651 1–34
  * [5] Castellano C, Fortunato S and Loreto V 2009 Reviews of modern physics 81 591
  * [6] Pei S and Makse H A 2013 Journal of Statistical Mechanics: Theory and Experiment 2013 P12002
  * [7] Liu L, Wang X, Zheng Y, Fang W, Tang S and Zheng Z 2019 New Journal of Physics
  * [8] Daley D J and Kendall D G 1964 Nature 204 1118–1118
  * [9] Li W, Tang S, Pei S, Yan S, Jiang S, Teng X and Zheng Z 2014 Physica A: Statistical Mechanics and its Applications 397 121–128
  * [10] Vega-Oliveros D A, Da F Costa L and Rodrigues F A 2017 Journal of Statistical Mechanics: Theory and Experiment 2017 023401
  * [11] Wang X, Li W, Liu L, Pei S, Tang S and Zheng Z 2017 Physical Review E 96 032304
  * [12] Moreno Y, Nekovee M and Pacheco A F 2004 Physical Review E 69 066130
  * [13] Centola D 2010 science 329 1194–1197
  * [14] Wang X, Tang S, Zheng Z and Fu F 2020 Physical Review X 10 041042
  * [15] Fu F, Christakis N A and Fowler J H 2017 Scientific reports 7 43634
  * [16] Romero D M, Meeder B and Kleinberg J 2011 Differences in the mechanics of information diffusion across topics: idioms, political hashtags, and complex contagion on twitter Proceedings of the 20th international conference on World wide web pp 695–704
  * [17] Centola D, Eguíluz V M and Macy M W 2007 Physica A: Statistical Mechanics and its Applications 374 449–456
  * [18] McAdam D and Paulsen R 1993 American journal of sociology 99 640–667
  * [19] Hodas N O and Lerman K 2014 Scientific reports 4 4343
  * [20] Aral S and Nicolaides C 2017 Nature communications 8 1–8
  * [21] Granovetter M S 1977 The strength of weak ties Social networks (Elsevier) pp 347–367
  * [22] Melnik S, Ward J A, Gleeson J P and Porter M A 2013 Chaos: An Interdisciplinary Journal of Nonlinear Science 23 013124
  * [23] Watts D J 2002 Proceedings of the National Academy of Sciences 99 5766–5771
  * [24] Centola D and Macy M 2007 American journal of Sociology 113 702–734
  * [25] Zheng M, Lü L, Zhao M et al. 2013 Physical Review E 88 012818
  * [26] Neuhäuser L, Mellor A and Lambiotte R 2020 Physical Review E 101 032310
  * [27] Gross T, D’Lima C J D and Blasius B 2006 Physical review letters 96 208701
  * [28] Iacopini I, Petri G, Barrat A and Latora V 2019 Nature communications 10 1–9
  * [29] Wang W, Liu Q H, Liang J, Hu Y and Zhou T 2019 Physics Reports 820 1–51
  * [30] Chen L, Ghanbarnejad F and Brockmann D 2017 New Journal of Physics 19 103041
  * [31] Soriano-Paños D, Ghanbarnejad F, Meloni S and Gómez-Gardeñes J 2019 Physical Review E 100 062308
  * [32] Sanz J, Xia C Y, Meloni S and Moreno Y 2014 Physical Review X 4 041005
  * [33] Pinotti F, Ghanbarnejad F, Hövel P and Poletto C 2020 Royal Society Open Science 7 190305
  * [34] Bródka P, Musial K and Jankowski J 2020 IEEE Access 8 10316–10341
  * [35] Wu Q and Chen S 2020 Chaos: An Interdisciplinary Journal of Nonlinear Science 30 073115
  * [36] Salehi M, Sharma R, Marzolla M, Magnani M, Siyari P and Montesi D 2015 IEEE Transactions on Network Science and Engineering 2 65–83
  * [37] Soriano-Paños D, Guo Q, Latora V and Gómez-Gardeñes J 2019 Physical Review E 99 062311
  * [38] De Domenico M, Granell C, Porter M A and Arenas A 2016 Nature Physics 12 901–906
  * [39] Bianconi G 2017 Journal of Statistical Mechanics: Theory and Experiment 2017 034001
  * [40] Chen L J, Chen X L, Cai M and Wang W 2018 Physica A: Statistical Mechanics and its Applications 503 516–525
  * [41] Cozzo E, Banos R A, Meloni S and Moreno Y 2013 Physical Review E 88 050801
  * [42] Gomez S, Diaz-Guilera A, Gomez-Gardenes J, Perez-Vicente C J, Moreno Y and Arenas A 2013 Physical review letters 110 028701
  * [43] Centola D, Becker J, Brackbill D and Baronchelli A 2018 Science 360 1116–1119
  * [44] Mønsted B, Sapieżyński P, Ferrara E and Lehmann S 2017 PloS one 12
  * [45] Mahajan V 2010 Wiley International Encyclopedia of Marketing
  * [46] Krapivsky P L, Redner S and Volovik D 2011 Journal of Statistical Mechanics: Theory and Experiment 2011 P12003
  * [47] Li W, Tang S, Fang W, Guo Q, Zhang X and Zheng Z 2015 Physical Review E 92 042810
  * [48] Wang W, Tang M, Stanley H E and Braunstein L A 2018 Physical Review E 98 062320
  * [49] Xie J, Sreenivasan S, Korniss G, Zhang W, Lim C and Szymanski B K 2011 Physical Review E 84 011130
  * [50] Karsai M, Iniguez G, Kaski K and Kertész J 2014 Journal of The Royal Society Interface 11 20140694
  * [51] Barabási A L and Albert R 1999 science 286 509–512
  * [52] Watts D J and Strogatz S H 1998 nature 393 440
  * [53] Moreno Y, Pastor-Satorras R and Vespignani A 2002 The European Physical Journal B-Condensed Matter and Complex Systems 26 521–529
  * [54] Dodds P S and Watts D J 2004 Physical review letters 92 218701
  * [55] Lambiotte R, Rosvall M and Scholtes I 2019 Nature physics 15 313–320
  * [56] Lehmann S and Ahn Y Y 2018 Complex spreading phenomena in social systems (Springer)


Also at ]Moscow Institute of Physics and Technology (State University),
Institutskii pereulok 9, Dolgoprudnyi 141701, Moscow region, Russia

#  An epsilon-near-zero-based Dallenbach absorber

Viacheslav V. Medvedev  [<EMAIL_ADDRESS>Institute of Spectroscopy of
the Russian Academy of Science, Fizicheskaya 5, Troitsk, Moscow 108840, Russia

###### Abstract

We report a theoretical analysis of total absorption conditions in a structure
consisting of a lossy coating layer on top of a specular metal substrate,
which is known as a Dallenbach absorber. All possible combinations of complex
refractive indices of the coating material and thicknesses of the coating
layer providing complete absorption are described for the case of normal
incidence of a plane wave. The work of coatings based on epsilon-near-zero
(ENZ) materials is analyzed in detail. Simple analytical relations are
obtained for the conditions of total absorption. It is shown that the range of
angles where strong absorption occurs is limited by the phenomenon of total
external reflection. The characteristics of coatings based on indium-tin
oxide, which is an example of an ENZ material, are analyzed.

electromagnetic wave absorption, thin film devices, antireflection coatings,
perfect absorber, epsilon-near-zero, conducting oxides

††preprint: AIP/123-QED

##  I Introduction

A Dallenbach absorber is one of the most striking and simplest devices
designed to effectively absorb electromagnetic radiation [1, 2, 3, 4].
According to the original design, it consists of a homogeneous lossy coating
layer placed on top of a reflecting metal substrate (see Fig. 1). Dallenbach
absorbers benefit from the destructive interference that minimize reflection
and, as a result, maximize absorption of radiation incident on a structure.
For a normally incident monochromatic plane wave with a wavelength $\lambda$,
one can write approximate equations for unity absorption conditions that
relate the optimal refractive index ($n_{l}$), extinction coefficient
($\kappa_{l}$) and thickness ($d_{l}$) of the coating layer:

$\displaystyle d_{l}=\frac{\lambda(2m+1)}{4n_{l}},$ (1a)
$\kappa_{l}=\frac{2}{\pi(2m+1)},$ (1b)

where $m=0,1,2,...$. These equations are derived for an idealized case when
the metal substrate is approximated by a perfect electric conductor (PEC). The
integer $m$ in Eq. 1b is referred to as the $m$th absorption mode. From the
practical point of view, absorbers with minimal coating layer thicknesses are
of considerable interest. To this end, the mode $m=0$ is most frequently used
in the literature. In this case, the optimal layer thickness and its
extinction coefficient are $d_{l}=\lambda/4n_{l}$ and $\kappa_{l}=2/\pi\approx
0.64$, respectively.

Figure 1: Schematic diagram of the Dallenbach absorber.

Obviously, the PEC model has its limitations. It works quite well at long
wavelengths, starting from the mid-IR range. In this case, the reflection
coefficient of real metals at normal incidence of a wave is close to unity,
and the phase of the reflected wave is close to $\pi$, as in PEC. In the near-
IR range and especially in the visible range, the difference between these
parameters for PEC and for real metals becomes significant if the interference
phenomena are taken into account. This also affects the conditions of complete
absorption for the Dallenbach layer. For instance, Kats et al. [5] showed that
for germanium coatings on a gold substrate, the optimal thickness in the
visible range turns out to be several times less than $\lambda/4n_{l}$. Park
et al. [6, 7] presented a detailed theoretical analysis for the case of a
substrate with optical constants that modulate real metals and for coating
materials with $n_{l}>1$. The authors of Ref. [6] came to the conclusion that
in the case of real metals, the optimal extinction coefficient of the coating
material, $\kappa_{l}$, exceeds the value of 0.64. These values of the complex
refractive index are well matched by various semiconductor materials. Cleary
et al. [8] have demonstrated strong IR absorption in coatings based on
gallium-doped zinc oxide (GZO) on a silver substrate. GZO is an example of
materials with a refractive index and extinction coefficient less than unity
and therefore can be classified as an ENZ materias. Badsha et al. [9] analyzed
in detail the conditions of total absorption for such coatings in the case of
grazing incidence. However, the case of normal incidence, which is of greatest
practical interest, has never been analyzed theoretically.

In this paper, we consider absorption of electromagnetic radiation by a
single-layer lossy coating on top of a metal substrate. Using numerical
calculations, we find all possible combinations of the complex refractive
index ($n_{l}+i\kappa_{l}$) of the coating material and its thickness
($d_{l}$), providing complete absorption of a plane monochromatic
electromagnetic wave for the idealized case of the PEC substrate. Based on
this analysis, we show that the condition of unity absorption described by Eq.
(1) is a special case corresponding to large values of $n_{l}$, i.e. $n_{l}\gg
1$. We also demonstrate that the unity absorption can be achieved for any
values of $n_{l}$, including $n_{l}$ near zero. In the latter case, analytical
expressions are also obtained for the conditions of total absorption:
$d_{l}=\lambda/2n_{l}$ and $\kappa_{l}=n_{l}^{2}/\pi$. An example of the
design of light-absorbing structures based on real index-near-zero material is
provided.

##  II Conditions for unity absorption

Consider the interaction of a plane electromagnetic wave having a wavelength
$\lambda$ with the model structure shown in Fig. 1. The wave is incident on
the structure at an angle $\theta_{a}$ measured from the normal. In the case
of opaque metal substrates, absorptance can be calculated by the expression

$A=1-R,$ (2)

where $R$ is the reflectance. In turn, reflectace can be calculated via the
amplitude reflection coefficient $r$ as $R=|r|^{2}$. The expression for $r$
has the form

$r=\frac{r_{al}+r_{ls}e^{2i\phi_{l}}}{1+r_{al}r_{ls}e^{2i\phi_{l}}},$ (3)

where $r_{al}$ and $r_{ls}$ are the Fresnel reflection coefficients for the
interfaces between the air and the coating layer and between the coating layer
and the substrate, respectively;
$\phi_{l}=2\pi\widetilde{n}_{l}d_{l}\cos\theta_{l}/\lambda$; and $\theta_{l}$
is the complex angle which is determined through the angle of incidence
$\theta_{a}$ via Snell’s law
$\widetilde{n}_{a}\sin\theta_{a}=\widetilde{n}_{l}\sin\theta_{l}$. We also
define $\delta_{l}=n_{l}d_{l}/\lambda$ as a reduced coating thickness which
will be used below in the analysis.

Let us analyze an idealized structure with a PEC substrate. In this case,
$r_{ls}=-1$ in Eq. 3. First, we consider the normal incidence of
electromagnetic radiation, i.e. $\theta_{a}=0$. The condition for unity
absorption is obtained from the roots of the numerator of Eq. 3:

$r_{al}-e^{2i\phi_{l}}=0,$ (4)

which can be rewritten in the form of two equations that must be
simultaneously satisfied:

$\displaystyle|r_{al}|=\exp(-4\pi\delta_{l}\kappa_{l}/n_{l}),$ (5a)
$\arg(r_{al})+2\pi m=4\pi\delta_{l},m=0,1,2,...$ (5b)

Equations 5b can be solved numerically for $\delta_{l}$ and $\kappa_{l}$ at a
fixed value of $n_{l}$. Thus, by varying $n_{l}$, one can find all possible
combinations of ($n_{l}$, $\kappa_{l}$, $\delta_{l}$) corresponding to unity
absorption. Figures 2(a)-(b) show a solution space for $m=0$. The green
horizontal dot-dashed lines in Figs. 2(a)-(b) denote the previously known
solution defined by Eqs. 1b, i.e. $\kappa_{l}=2/\pi$ and $\delta_{l}=0.25$. It
is seen that the numerical solution obtained here asymptotically tends to the
indicated values only for large values of $n_{l}$. For $n_{l}<2$, the
numerical solution differs significantly from that defined by Eqs. 1b. As
$n_{l}$ tends to zero, $\delta_{l}$ tends to a constant value of 0.5, and
$\kappa_{l}$ tends to zero approximately as $0.32n_{l}^{2}$ (see black dashed
curve in Fig. 2(a)).

|
---|---
(a) | (b)
Figure 2: (a)-(b): Solid red curves demonstrate the calculated combinations of parameters ($n_{l},\kappa_{l},\lambda_{l}$) of the Dallenbach layer on a PEC substrate, providing complete absorption of a normally incident plane wave. Horizontal green dot-dashed lines are classic solutions for the Dallenbach layer represented by Eqs. 1 at $m=0$. (a) Dependence of the optimal extinction coefficient $\kappa_{l}$ of the coating material on its refractive index $n_{l}$. Black dashed curve shows the fitting function $\kappa_{l}=0.32n_{l}^{2}$. (b) Dependence of the optimal normalized coating thickness $\delta_{l}$ on the refractive index of the coating material $n_{l}$.  |  |
---|---|---
(a) | (b) | (c)
Figure 3: (a)-(c): Calculated profiles of the electric field for Dallenbach layers with different values of the refractive index of the coating material: (a) $n_{l}=10$, (b) $n_{l}=1$, and (c) $n_{l}=0.1$. |  |
---|---|---
(a) | (b) | (c)
Figure 4: (a)-(c): Calculated angular dependencies of the absorbance for
Dallenbach layers with different values of the refractive index of the
coating: (a) $n_{l}=10$, (b) $n_{l}=1$, and (c) $n_{l}=0.1$. Dashed orange
curves correspond to TM polarization, and solid purple curves correspond to TE
polarization.

Let us analyze Eq. 5b(a) under the conditions $n_{l}\rightarrow 0$,
$\kappa_{l}<n_{l}$ and $\delta_{l}\rightarrow 0.5$. In this case, we obtain
the following approximation for the Fresnel reflection coefficient for the
upper interface of the structure:
$r_{al}=(1-\widetilde{n}_{l})/(1+\widetilde{n}_{l})\approx(1-n_{l})/(1+n_{l})\approx
1-2n_{l}$. Equation 5a then transforms into
$1-2n_{l}=\exp(-2\pi\kappa_{l}/n_{l})$. Expanding the exponent in this
equation into a Taylor series in powers of $\kappa_{l}/n_{l}$, we obtain the
solution $\kappa_{l}\approx n_{l}^{2}/\pi\approx 0.32n_{l}^{2}$.

It is also interesting to analyze the distribution of the electric field
inside the considered absorber. The field within a layer can be calculated
using the expressions:

$\displaystyle E_{a}(z)=E_{0}e^{-ik_{a}z}+rE_{0}e^{ik_{a}z},z>d_{l}$ (6a)
$\displaystyle
E_{l}(z)=\frac{2it_{al}e^{i(k_{l}-k_{a})d_{l}}E_{0}}{1+r_{al}r_{lm}e^{2ik_{l}d_{l}}}\sin(k_{l}z),0<z\leq
d_{l}$ (6b) $\displaystyle E_{s}(z)=0,z\leq 0,$ (6c)

where $E_{0}$ is the amplitude of the incident wave;
$k_{j}=2\pi\widetilde{n}_{j}/\lambda$; and $j=a,l,s$. Figures 3(a)-(c) show
the calculated profiles of the modulus of the electric field for three
structures with different values of the refractive index: $n_{l}=10$,
$n_{l}=1$, and $n_{l}=0.1$. The values of $\delta_{l}$ and $\kappa_{l}$ were
chosen in such a way that all three structures provide complete absorption of
the incident wave. It follows from Fig. 3(c) that in the case of $n_{l}<1$,
the field is enhanced inside the absorbing layer. Analysis of Eq. 6 at
$n_{l}\rightarrow 0$, $\kappa_{l}\approx n_{l}^{2}/\pi$ and
$\delta_{l}\rightarrow 0.5$ shows that the field amplitude inside the layer
inversely depends on $n_{l}$.

Let us now analyze the angular dependence of absorptance by the structures in
question at different values of the refractive index. Figures 4(a)-(c) show
the calculated angular dependences of the absorption coefficient for three
structures with $n_{l}=10$, $n_{l}=1$, and $n_{l}=0.1$. The values of
$\delta_{l}$ and $\kappa_{l}$ were chosen in such a way that all three
structures provide complete absorption of the wave at normal incidence. The
purple curves in Figs. 4(a)-(c) correspond to the TE polarization of the
incident plane wave, and orange ones, to the TM polarization. It follows from
Figs. 4(a)-(c) that for $n_{l}=10$ and $n_{l}=1$, absorption decreases
relatively slow with increasing angle $\theta_{a}$. When $n_{l}=0.1$, strong
absorption occurs only at incidence angles of several degrees, which is due to
the fact that at $n_{l}<1$, total external reflection at oblique incidence
occurs at angles $\theta_{a}>\arcsin(n_{l})$. In particular, for $n_{l}=0.1$,
total external reflection occurs for angles of incidence exceeding
$\arcsin(0.1)\approx 5.74^{\circ}$. In the angular dependence of the
absorptance for the TM wave, a narrow peak is observed in the region of the
angle $\theta_{a}\approx 5.73^{\circ}$, which is close in value to the
Brewster angle for the upper structure interface, i.e. $\arctan(0.1)\approx
5.71^{\circ}$.

##  III Real materials

|  |
---|---|---
(a) | (b) | (c)
Figure 5: (a) Refractive index (solid green curve) and extinction coefficient
(dot-dashed orange curve) of ITO at $N_{e}=10^{21}$ cm-3. (b) Reflectance of
the ITO coating on the Al substrate as a function of the wavelength and
thickness of the ITO ($N_{e}=10^{21}$ cm-3) coating. (c) Absorption spectrum
of the ITO ($N_{e}=10^{21}$ cm-3) coating on the Al substrate for the coating
having a thickness of 0.56 $\mu$m.

The results presented in Fig. 1(b) show that light-absorbing coatings can be
made of a variety of materials. Coatings based on materials with refractive
indices greater than unity have been extensively studied in previous works [5,
7, 6, 10, 11, 12, 13, 14, 15, 16]. Here we consider in detail coatings based
on ENZ materials. Examples of such materials are transparent conductive oxides
such as indium tin oxide (ITO) and aluminum zinc oxide (AZO). Such oxides are
usually characterized by high electron concentrations, $N_{e}\sim
10^{19}-10^{21}$ cm -3, resulting from the deposition conditions. Such high
concentrations of charge carriers make it possible to obtain the ENZ mode in
the infrared range. The optical properties of such materials are well
described by the Drude model:

$\varepsilon=\varepsilon_{\infty}-\frac{\omega_{p}^{2}}{\omega^{2}+i\omega\Gamma},$
(7)

where $\varepsilon$ denotes the complex permittivity of the material at a
frequency $\omega=2\pi c/\lambda$, $\varepsilon_{\infty}$ is the permittivity
at an infinite frequency, $\omega_{p}$ is the electron plasma frequency, and
$\Gamma$ is the electron collision frequency. The last two quantities can be
calculated by the formulas:
$\omega_{p}=\sqrt{N_{e}e^{2}/m^{*}\varepsilon_{0}}$ and $\Gamma=e/m^{*}\mu$.
Here $e$ is the elementary charge, $m^{*}$ is the electron effective mass,
$\mu$ is the electron mobility, and $\varepsilon_{0}$ is the permittivity of
free space.

Let us analyze the characteristics of ITO-based coatings. The parameters
required for calculating the ITO optical constants by the Drude model were
borrowed from Ref. [17]. We choose aluminum (Al) as the substrate material.
The optical constants of Al for the calculations were borrowed from Ref. [18].
To begin with, let the concentration of electrons be $N_{e}=10^{21}$ cm-3.
Figure 5(a) shows the complex refractive index of ITO, calculated according to
the Drude model at a given electron concentration for wavelengths in the range
of 0.5 - 1.5 $\mu$m. One can see from Fig. 5(a) that the ENZ mode is realized
near $\lambda=1$ $\mu$m. Figure 5(b) demonstrates the calculated reflectance
of the ITO/Al structure as a function of the wavelength and thickness of the
ITO coating. According to the calculations, at point with a wavelength of
$\lambda_{0}\approx 0.99$ $\mu$m and a thickness of $d_{0}\approx 0.56$
$\mu$m, the reflection vanishes, which is equivalent to total absorption. The
complex refractive index of ITO at a given electron concentration at a
wavelength of 0.99 $\mu$m is
$\widetilde{n}_{\text{ITO}}=n_{\text{ITO}}+i\kappa_{\text{ITO}}=0.74+0.22i$.
According to the above predictions, the optimal extinction coefficient for a
structure with a PEC substrate is $n_{\text{ITO}}^{2}/\pi\approx 0.18$, which
is close to the extinction coefficient of the considered material
$\kappa_{\text{ITO}}=0.22$. The normalized optical thickness of the coating in
this case is $\delta=d_{0}n_{\text{ITO}}/\lambda_{0}\approx 0.42$, which is
also close to the analytical predictions for the case with the PEC substrate,
i.e. $\delta=0.5$. Note that there are also other points with a larger coating
thickness and lower wavelength values, at which zero reflection is
implemented. However, we restrict ourselves to considering structures with a
minimum resonant coating thickness. The calculated absorption spectrum for a
structure with an ITO coating having a thickness of 0.56 $\mu$m is shown in
Fig. 5(c).

What happens when the concentration of electrons in the coating material
changes? For each value of Ne, we can also find a point at which the
reflection vanishes. Figure 6 shows the dependences of $\lambda_{0}$ and
$d_{0}n_{ITO}/\lambda_{0}$ on $N_{e}$ in the range $10^{19}-10^{21}$ cm-3. It
can be seen from Fig. 6 that with a decrease in the electron concentration
from $10^{21}$ cm-3 to $10^{19}$ cm-3, the resonance wavelength increases from
approximately 1 $\mu$m to 10 $\mu$m, and the normalized optical thickness of
the coating will decrease from approximately 0.42 to 0.29. Note also that
changes in the substrate material do not lead to significant quantitative
changes in the values of $\lambda_{0}$ and $d_{0}$.

|
---|---
(a) | (b)
Figure 6: (a) Dependence of the resonant wavelength $\lambda_{0}$ for the ITO
coating atop the Al substrate on the electron concentration in ITO. (b)
Dependence of the normalized optical thickness
$d_{0}n_{\text{ITO}}/\lambda_{0}$ of the ITO coating atop the Al substrate on
the electron concentration in ITO.

##  IV Conclusions

We have investigated absorption of electromagnetic radiation in the so-called
Dallenbach layer, resulting from the phenomenon of destructive interference.
We have shown that the Dallenbach layer can provide complete absorption of a
normally incident plane wave in a wide range of complex refractive indices of
the coating layer material, including close-to-zero indices. This fact makes
it possible to design efficient ENZ-based absorbers. Simple analytical design
rules have been derived, demonstrating that such structures provide high
angular absorption selectivity, which is determined by total external
reflection. An example of the design of a tunable wavelength absorber based on
an ITO layer on an Al substrate have been presented.

## References

  * Dallenbach and Kleinsteuber [1938] W. Dallenbach and W. Kleinsteuber, “Reflection and absorption of decimeter-waves by plane dielectric layers,” Hochfrequenztechnik und Elektroakustik 51, 152–156 (1938).
  * Knott, Shaeffer, and Tuley [2004] E. F. Knott, J. F. Shaeffer, and M. T. Tuley, “Radar cross section,” (SciTech Publishing, Inc, 2004) Chap. Radar Absorbing Materials, pp. 297–360, 2nd ed.
  * Saville [2005] P. Saville, “Review of radar absorbing materials,” Tech. Rep. (Defence Research and Development Atlantic Dartmouth (Canada), 2005).
  * Tretyakov [2016] S. Tretyakov, “Thin absorbers: operational principles and various realizations,” IEEE Electromagnetic Compatibility Magazine 5, 61–66 (2016).
  * Kats  _et al._ [2013] M. A. Kats, R. Blanchard, P. Genevet, and F. Capasso, “Nanometre optical coatings based on strong interference effects in highly absorbing media,” Nature materials 12, 20 (2013).
  * Park, Kim, and Brongersma [2015] J. Park, S. J. Kim, and M. L. Brongersma, “Condition for unity absorption in an ultrathin and highly lossy film in a gires–tournois interferometer configuration,” Optics letters 40, 1960–1963 (2015).
  * Park  _et al._ [2014] J. Park, J.-H. Kang, A. P. Vasudev, D. T. Schoen, H. Kim, E. Hasman, and M. L. Brongersma, “Omnidirectional near-unity absorption in an ultrathin planar semiconductor layer on a metal substrate,” Acs Photonics 1, 812–821 (2014).
  * Cleary  _et al._ [2015] J. W. Cleary, N. Nader, K. D. Leedy, and R. Soref, “Tunable short-to mid-infrared perfectly absorbing thin films utilizing conductive zinc oxide on metal,” Optical Materials Express 5, 1898–1909 (2015).
  * Badsha, Jun, and Hwangbo [2014] M. A. Badsha, Y. C. Jun, and C. K. Hwangbo, “Admittance matching analysis of perfect absorption in unpatterned thin films,” Optics Communications 332, 206–213 (2014).
  * Wang  _et al._ [2018] H.-C. Wang, C. H. Chu, P. C. Wu, H.-H. Hsiao, H. J. Wu, J.-W. Chen, W. H. Lee, Y.-C. Lai, Y.-W. Huang, M. L. Tseng, _et al._ , “Ultrathin planar cavity metasurfaces,” Small 14, 1703920 (2018).
  * Liu  _et al._ [2016] D. Liu, H. Yu, Z. Yang, and Y. Duan, “Ultrathin planar broadband absorber through effective medium design,” Nano Research 9, 2354–2363 (2016).
  * Wen and Xiong [2016] X. Wen and Q. Xiong, “A large scale perfect absorber and optical switch based on phase change material (ge 2 sb 2 te 5) thin film,” Science China Materials 59, 165–172 (2016).
  * Mirshafieyan and Guo [2014] S. S. Mirshafieyan and J. Guo, “Silicon colors: spectral selective perfect light absorption in single layer silicon films on aluminum surface and its thermal tunability,” Optics express 22, 31545–31554 (2014).
  * Mirshafieyan, Guo, and Guo [2016] S. S. Mirshafieyan, H. Guo, and J. Guo, “Zeroth order fabry–perot resonance enabled strong light absorption in ultrathin silicon films on different metals and its application for color filters,” IEEE Photonics Journal 8, 1–12 (2016).
  * Krayer  _et al._ [2018] L. J. Krayer, E. M. Tennyson, M. S. Leite, and J. N. Munday, “Near-ir imaging based on hot carrier generation in nanometer-scale optical coatings,” ACS Photonics 5, 306–311 (2018).
  * Dias  _et al._ [2018] M. R. S. Dias, C. Gong, Z. A. Benson, and M. S. Leite, “Lithography-free, omnidirectional, cmos-compatible alcu alloys for thin-film superabsorbers,” Advanced Optical Materials 6, 1700830 (2018).
  * Anopchenko  _et al._ [2018] A. Anopchenko, L. Tao, C. Arndt, and H. W. H. Lee, “Field-effect tunable and broadband epsilon-near-zero perfect absorbers with deep subwavelength thickness,” ACS Photonics 5, 2631–2637 (2018).
  * Rakić [1995] A. D. Rakić, “Algorithm for the determination of intrinsic optical constants of metal films: application to aluminum,” Applied optics 34, 4755–4767 (1995).


# Machine learning for potion development at Hogwarts

Christoph F. Kurz1,2, Adriana N. König1,2

(1Institute of Health Economics and Health Care Management, Helmholtz Zentrum
München, Ingolstädter Landstraße 1, 85764 Neuherberg, Germany
2Munich School of Management and Munich Center of Health Sciences, Ludwig-
Maximilians-Universität München, Geschwister-Scholl-Platz 1, 80539 Munich,
Germany
)

###### Abstract

Objective: To determine whether machine learning methods can generate useful
potion recipes for research and teaching at Hogwarts School of Witchcraft and
Wizardry.
Design: Using deep neural networks to classify generated recipes into a
standard drug classification system.
Setting: Hogwarts School of Witchcraft and Wizardry.
Data sources: 72 potion recipes from the Hogwarts curriculum, extracted from
the Harry Potter Wiki.
Results: Most generated recipes fall into the categories of psychoanaleptics
and dermatologicals. The number of recipes predicted for each category
reflected the number of training recipes. Predicted probabilities were often
above 90% but some recipes were classified into 2 or more categories with
similar probabilities which complicates anticipating the predicted effects.
Conclusions: Machine learning powered methods are able to generate potentially
useful potion recipes for teaching and research at Hogwarts. This corresponds
to similar efforts in the non-magical world where such methods have been
applied to identify potentially effective drug combinations.

## Introduction

Potions are a required subject for students at Hogwarts School of Witchcraft
and Wizardry from the first to the fifth year [1]. They are optional to
students in their sixth and seventh-years if they achieved a high score on
their Ordinary Wizarding Level exam [2]. Potions classes are considered to be
among the most difficult lessons at Hogwarts because the nuances of timing,
ageing, bottling, and stirring techniques are difficult to acquire even with
the guidance of experienced teachers such as Professor Severus Snape. Brewing
potions requires glass vials, weighting scales and a cauldron. Ingredients
range from plants such as belladonna and shrivelfig to magical components such
as unicorn hair or fairy wings. The brewing process often requires some degree
of wand work to complete [1].

Potions can be used as medicines, antidotes, poisons, or to provide the
drinker with any magical effect ranging from increased strength to flame
immunity. They are not always consumed by drinking them; some, like the
Regeneration potion, might be applied by physical touch or have an effect
merely by being created [3]. Certain magical effects can only be achieved
through the use of potions. Some potions mimic the effects of spells and
charms, but a few (such as the Polyjuice Potion, a potion that allows the
drinker to take the form of someone else, and Felix Felicis, a luck potion)
have effects that can not be achieved in any other way [4, 2].

Because brewing is so difficult and the smallest deviations from the recipe
can have serious consequences, there are countless reports of accidents and
undesirable side effects happening during class at Hogwarts [1, 4, 5]. For
example, in 1992, there was a well-documented case of the student Neville
Longbottom who, while improperly brewing the Cure for Boils potion, infected
himself with red boils all over his body [1]. Nevertheless, some deviations
from instructions have proven successful for Harry Potter in his fifth year at
Hogwarts [2]. Accurately following the brewing instructions is already
difficult, but the discovery and development of new potions is an even more
complex and dangerous process.

Recent advances in the field of artificial intelligence (AI) have led to
increased interest in the use of machine learning approaches within the
pharmaceutical industry. Advances in new algorithms, such as deep neural
networks, demonstrated its utility in addressing diverse problems in drug
discovery such as bioactivity prediction or novel molecular designs [6, 7].

In this work, we explore the usefulness of machine learning for generating
recipes for magic potions. For this, we randomly generated new magic potion
recipes with various ingredients and predicted their most likely effect using
an artificial neural network.

## Methods

We collected the recipes for all known potions from the Harry Potter Wiki [8]
and classified them according to the Anatomical Therapeutic Chemical (ATC)
classification system [9] in one of the following categories: anesthetics;
antiinfectives for systemic use; antiparasitic products, insecticides and
repellents; dermatologicals; musculo-skeletal system; psychoanaleptics;
psycholeptics; respiratory system; sensory organs; and various. These
categories represent the first and second level of the ATC classification
which describes pharmacological or therapeutic subgroups [9]. Recipes in the
musculo-skeletal system category include, for example, the pompion potion that
temporarily turns the drinkers head into a pumpkin, or the skelegro potion
that regrows bones. Dermatologicals are, among others, potions that make your
skin immune to fire, grow hair or cure boils. Recipes in the psychoanaleptics
category include, for example, the forgetfulness potion which causes memory
loss in the drinker. Others are the wit-sharpening potion which improves clear
thinking or the befuddlement draught that provokes belligerence and
recklessness in the drinker. The various category contains several antidotes
as well as potions that boost spell-casting, such as the exstimulo potion.

We additionally added a category for poisons because many recipes fall into
it. However, poisons are not associated with an ATC code for obvious reasons.
Usually, administered drugs aim at improving and not deteriorating an
individual’s health. See Table 1 for an overview of the number of recipes in
each of the 11 categories. In total, the training set contained 72 recipes.
Each recipe includes instructions for adding ingredients and brewing. We then
generated 10,000 new potion recipes by randomly picking between 3 to 8 single
ingredients (e.g., “Add 4 horned slugs to your cauldron”) and mixing
instructions (e.g., “Stir 5 times, clockwise.”). We used a custom BioBERT
neural network [10] for predicting the class of potion for each newly
generated recipe. This language model has been pre-trained on large-scale
biomedical corpora comprising over 18 billion words from PubMed abstracts and
full-text articles. We fine-tuned the model to all known Hogwarts potion
recipes so it would input a recipe and output the probabilities of belonging
to each of the 11 classes. The top probability is the most likely effect. This
method is often referred to as feature extraction transfer learning [11]. All
computations were done in Mathematica 13 [12]. Code and data are available on
our GitHub page [13].

Category | ATC Code | Count
---|---|---
Anesthetics | N01 | 2
Antiinfectives for systemic use | J | 4
Antiparasitic products, insecticides and repellants | P | 2
Dermatologicals | D | 8
Musco-sceletal system | M | 12
Poison | – | 11
Psychoanaleptics | N06 | 21
Psycholeptics | N05 | 2
Respiratory system | R | 1
Sensory organs | S | 1
Various | V | 8
Table 1: Number of recipes and ATC codes for each of the 11 categories in the
training data.

## Results

Figure 1 shows the number of predicted potion recipes in each category. Most
of the 10,000 generated recipes fall into the psychoanaleptics category
($n=5549$), followed by the dermatologicals category ($n=1539$) and the
various category ($n=1487$). 225 recipes fall into the newly added poison
category. In contrast, only 3 psycholeptics, 3 respiratory systems, and 1
sensory organs recipes were generated. This corresponds to the number of
available training recipes. All generated recipes differed from the training
set of recipes. Our BioBERT model was generally confident in its predictions.
Predicted probabilities of belonging to a certain ATC category were often
above 90%, see Figure 2.

Figure 1: Predicted potion categories of the 10,000 randomly generated
recipes. Figure 2: Count histogram of the predicted probabilities of belonging
to one of the 11 ATC categories (including poisons). Only the top probability
for each generated recipe is shown.

For example, Table 2 shows a generated recipe where its predicted effect is in
the psychoanaleptics category with a probability of 99.9%. In contrast, the
effects of some recipes are difficult to predict for our model. Table 3 shows
a generated recipe that could be dermatological with 58.4% probability, be a
psychoanaleptic with 10% probability, or an antiinfective for systemic use
with 24.1% probability.

Add Baneberry. Add 2 bundles of knotgrass to the cauldron. Add Dogbane. Add
syrup of hellebore until the potion turns turquoise. Add a sprig of Peppermint
to counteract side-effects. Add Honey water until it turns back to a turquoise
colour. Stir four times anti-clockwise. Add the Infusion of Wormwood.
---
Table 2: Generated recipe that is predicted to work as a psychoanaleptic with
99.9% probability. Add Stewed Mandrake. Add Wormwood. Add a dash of
Flobberworm Mucus and stir vigorously. Leave to brew and return in 8 hours
(Copper), 14 hours (Brass), or 23 hours (Pewter). Shake and add wormwood until
the potion turns green. Slice bursting mushrooms with knife, add to cauldron
and stir clockwise until potion turns blue.
---
Table 3: Generated recipe that is predicted to work as a dermatological with
58.4% probability, as a psychoanaleptic with 10% probability, and as an
antiinfective for systemic use with 24.1% probability.

## Discussion

Our findings suggest that AI powered methods are able to generate potentially
useful potion recipes for teaching and research at Hogwarts School of
Witchcraft and Wizardry. We were able to produce many previously unknown
combinations of ingredients and stirring instructions that were predicted to
belong to a specific bioactivity class with high probability. Previously, AI
methods have also been used to identify potentially effective drug
combinations [14, 15]. In the magical world, our research could be extended to
not only detect new combinations of ingredients but also new combinations of
potions. Apart from new effective combinations, AI methods could also be
applied to identify potentially harmful drug combinations [14]. This
complements our predictions in the category for poisons and could be extended
to harmful potion combinations.

Still, our results are not without limitations. In general, AI models need
very large training sets. We only had a set of 72 recipes available for
training. For this reason, we used a model that has been pre-trained on large
medical corpora. Still, potions belonging to the same ATC category often have
very different effects. For example, both the Babbling Potion, a potion that
causes the drinker to babble nonsense, and Baruffio’s Brain Elixir, a potion
that increases the drinker’s brain power, are part of the nervous system
category. This makes it extremely difficult to predict specific potion
effects, other than the organ or system on which they act.

Furthermore, our AI approach for drug discovery and potion generation could
potentially be misused. For example, Urbina et al. [16] trained an AI model
that generated new molecules that were predicted to be more toxic than
publicly known chemical warfare agents. In this sense, machine learning could
be used to support the Dark Arts. The Dark Arts refer to spells and actions
that could harm others, such as powerful curses, as well as brewing dark
potions and breeding dark creatures. Our approach could lead to the discovery
of new spells and potions that would enable Dark Wizards or Witches becoming
even more powerful than Lord Voldemort, considered to have been the most
capable and dangerous practitioner of the Dark Arts of all time [2].

At last, two muggles with (presumably) no magical abilities performed the
study. Thus, it is difficult to assess the validity and classification quality
of the generated recipes.

### Competing interests statement

All authors have completed the Unified Competing Interest form (available on
request from the corresponding author) and declare: no support from any
organisation for the submitted work; no financial relationships with any
organisations that might have an interest in the submitted work in the
previous three years, no other relationships or activities that could appear
to have influenced the submitted work.

### Details of contributors

CFK and ANK collected the data and wrote the manuscript. CFK analysed the
data. ANK is the guarantor.

### Transparency declaration

The lead author (the manuscript’s guarantor) affirms that this manuscript is
an honest, accurate, and transparent account of the study being reported; that
no important aspects of the study have been omitted; and that any
discrepancies from the study as planned have been explained.

### Ethical approval

The data were obtained from publicly available data sources.

### Details of funding

No funding was received.

### Patient and public involvement statement

Not applicable.

## References

  * Rowling [1997] Rowling JK.  Harry Potter and the philosopher’s stone.  Bloomsbury Publishing; 1997.
  * Rowling [2005] Rowling JK.  Harry Potter and the half-blood prince.  Bloomsbury publishing; 2005.
  * Rowling [2000] Rowling JK.  Harry Potter and the goblet of fire.  Bloomsbury publishing; 2000.
  * Rowling [1998] Rowling JK.  Harry Potter and the chamber of secrets.  Bloomsbury publishing; 1998.
  * Rowling [1999] Rowling JK.  Harry Potter and the prisoner of Azkaban.  Bloomsbury publishing; 1999.
  * Chen et al. [2018] Chen H, Engkvist O, Wang Y, Olivecrona M, Blaschke T.  The rise of deep learning in drug discovery.  Drug discovery today 2018;23(6):1241–1250.
  * Vamathevan et al. [2019] Vamathevan J, Clark D, Czodrowski P, Dunham I, Ferran E, Lee G, et al.  Applications of machine learning in drug discovery and development.  Nature reviews Drug discovery 2019;18(6):463–477.
  * Harry Potter Wiki [2022] Harry Potter Wiki, List of Potions; 2022.  https://harrypotter.fandom.com/wiki/List_of_potions, accessed on 2022-03-22.
  * WHO Collaborating Centre for Drug Statistics Methodology [2022] WHO Collaborating Centre for Drug Statistics Methodology.  ATC classification index with DDDs 2022; 2022.
  * Lee et al. [2019] Lee J, Yoon W, Kim S, Kim D, Kim S, So CH, et al.  BioBERT: a pre-trained biomedical language representation model for biomedical text mining.  Bioinformatics 2019 09;36(4):1234–1240.
  * Goodfellow et al. [2016] Goodfellow I, Bengio Y, Courville A.  Deep Learning.  Adaptive Computation and Machine Learning series, MIT Press; 2016.
  * Wolfram Research [2021] Wolfram Research, Mathematica, Version 13.0.0; 2021.  https://www.wolfram.com/mathematica, Champaign, IL.
  * Kurz [2022] Kurz C, Harry Potter potion discovery.  GitHub; 2022.  https://github.com/krz/hp-potion-discovery.
  * Güvenç Paltun et al. [2021] Güvenç Paltun B, Kaski S, Mamitsuka H.  Machine learning approaches for drug combination therapies.  Briefings in Bioinformatics 2021;22(6):bbab293.
  * Wu et al. [2022] Wu L, Wen Y, Leng D, Zhang Q, Dai C, Wang Z, et al.  Machine learning methods, databases and tools for drug combination prediction.  Briefings in bioinformatics 2022;23(1):bbab355.
  * Urbina et al. [2022] Urbina F, Lentzos F, Invernizzi C, Ekins S.  Dual use of artificial-intelligence-powered drug discovery.  Nature Machine Intelligence 2022;p. 1–3.