# Amplitude death in coupled replicator map lattice: averting migration
dilemma

Shubhadeep Sadhukhan<EMAIL_ADDRESS>Department of Physics, Indian Institute
of Technology Kanpur, Uttar Pradesh 208016, India  Rohitashwa Chattopadhyay
<EMAIL_ADDRESS>Department of Physics, Indian Institute of Technology
Kanpur, Uttar Pradesh 208016, India  Sagar Chakraborty<EMAIL_ADDRESS>Department of Physics, Indian Institute of Technology Kanpur, Uttar Pradesh
208016, India

###### Abstract

Populations composed of a collection of subpopulations (demes) with random
migration between them are quite common occurrences. The emergence and
sustenance of cooperation in such a population is a highly researched topic in
the evolutionary game theory. If the individuals in every deme are considered
to be either cooperators or defectors, the migration dilemma can be envisaged:
The cooperators would not want to migrate to a defector-rich deme as they fear
of facing exploitation; but without migration, cooperation can not be
established throughout the network of demes. With a view to studying the
aforementioned scenario, in this paper, we set up a theoretical model
consisting of a coupled map lattice of replicator maps based on two-
player–two-strategy games. The replicator map considered is capable of showing
a variety of evolutionary outcomes, like convergent (fixed point) outcomes and
nonconvergent (periodic and chaotic) outcomes. Furthermore, this coupled
network of the replicator maps undergoes the phenomenon of amplitude death
leading to non-oscillatory stable synchronized states. We specifically explore
the effect of (i) the nature of coupling that models migration between the
maps, (ii) the heterogenous demes (in the sense that not all the demes have
same game being played by the individuals), (iii) the degree of the network,
and (iv) the cost associated with the migration. In the course of
investigation, we are intrigued by the effectiveness of the random migration
in sustaining a uniform cooperator fraction across a population irrespective
of the details of the replicator dynamics and the interaction among the demes.

##  I Introduction

Large spatially distributed populations more often than not form clusters of
several subpopulations connected through migration which is one of the
important mechanisms in shaping the evolution and bringing forth emergence of
cooperation [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12]. It is easy to envisage
that availability of better opportunities [11] elsewhere lead individuals to
abandon their home and migrate. Moreover, individuals may want to migrate to
satisfy their aspirations [5, 7]. The migration can furthermore depend on the
expectations [9] of the individuals. Risk-driven migration [10] and success-
driven migration [3, 6] can also promote cooperation effectively. While a
random migration arguably weakens the emergence of cooperation by favouring
the invasion by defection [13, 14], very mildly incentivizing cooperating
behaviour [15] can overcome this drawback of the random migration.

Theoretically, emergence and sustenance of cooperation in a collection of
subpopulations or demes with active migration between them can be conveniently
studied within the paradigm of the evolutionary game theory by using coupled
map lattice (CML) models [16, 17]. In fact, in such a setting, the phenomena
of cooperation, chaos, and synchronization come together to overcome the
migration dilemma [15]: In a CML of subpopulations of replicators with two
actions—cooperate or defect—if all but a few subpopulations have defectors
exclusively, the cooperators would not want to migrate lest they should be
exploited by the defectors; however, in the absence of any migration,
cooperation would not be established across the network of subpopulations and
therefore the collective utility gain for the population is denied in the
light of cooperators’ not risking their relatively higher payoffs.

The nonlinear dynamics and network dynamics of the evolutionary systems in the
context of the interplay between synchronization and cooperation have
motivated quite a few recent studies, e.g., the ones on the evolutionary
Kuramoto dilemma [18, 19, 20] and the one on chaotic agent dynamics [21]. In
the setting of the CML with chaotic replicator maps, the amplitude variations
of the chaotic oscillations of the fraction of the cooperators in the
subpopulations are suppressed due to synchronization onto a fixed point of the
CML [15]. It is natural to draw an analogy with the amplitude death [22] in
coupled oscillators whenever the CML synchronizes onto a fixed point [23, 24,
25, 26, 27, 28, 29]. Technically, the term amplitude death refers to the
situation when the oscillations—either periodic or quasi-periodic or
chaotic—of an entire system of coupled oscillators ceases, leading to the
stationarity [30, 31, 32]. The amplitude death occurs in a wide variety of
systems, whether interacting systems are identical [33, 34, 35, 36],
mismatched [32, 37], dynamically coupled [38], or nonlinear [39]. Apart from
diffusive coupling, the nonlinear coupling is also used in achieving the
amplitude death [39, 40]. In the presence of the nonlinear coupling, the
amplitude death occurs in the absence of parameter mismatch and also in the
absence of time delay [39, 41].

In this paper, we investigate the migration dilemma in the setup of the CML of
replicators’ subpopulations and investigate the nonlinear dynamics of the
amplitude death induced synchronization. Specifically, we ask the following
relevant questions: Firstly, what happens if the interaction between the demes
is not simple diffusive coupling but some more general nonlinear coupling?
Secondly, does the amplitude death synchronizes the replicator dynamics at the
demes that have non-chaotic but periodic dynamics? Thirdly, how is cooperation
across the population supported as the degree of the network varies?
Furthermore, costly interactions in networks while studying the coevolution of
synchronization and cooperation has gathered a lot of recent interest [18, 21,
20]. In such a scenario, the deme—based on the payoff it receives—can decide
whether to participate in the process of migration or not. Migration—usually
costly for different reasons—is seen in the population of different kind of
insects [42], birds [43], fishes [44], and mammals [45]. One can find
interesting investigations about migration cost in the populations—such as
that of insects [42], white storks [46], and spoonbills [47, 48]. In the light
of these studies, in our setup, another crucial theoretical question can be
addressed: How a costly interdemic interaction affects the cooperation levels
in the network? It is of interest for us to understand the effect (if present
at all) of such costly migration on averting the migration dilemma via the
amplitude death.

Before we present the results that we have found while investigating the
above-mentioned problems, let us first succinctly set up the CML—on which the
results are based—in the following section.

##  II The coupled map lattice

Mathematically, our model comprises of a CML network whose nodal dynamics is
governed by a replicator map, which for the two-player–two-strategy games are
the most convenient yet nontrivial proving ground for our idea. The one-
dimensional replicator map [49, 50, 51, 52, 53, 54, 55, 56, 57] is given by,

$x_{n+1}=f(x_{n}):=x_{n}+x_{n}[({\sf
A}\mathbf{x_{n}})_{1}-(\mathbf{x_{n}})^{T}{\sf A}\mathbf{x_{n}}],$ (1)

where subscript ‘1’ denotes the first component of vector ${\sf
A}\mathbf{x_{n}}$, $n$ denotes the time step, and

exhibits the strategies involved and the real-valued payoff matrix ${\sf A}$
for Player 1 in the two-player–two-strategy symmetric game. ${\bf x}=(x,1-x)$
is the state of the population such that $x$ is the fraction of the
cooperators and $1-x$ is the fraction of the defectors. For consistency, it is
required that the one-dimensional replicator map be such that $0\leq x_{n}\leq
1$ for all $n$. This strictly depends on the values of the elements of the
payoff matrix. It should be borne in mind that the discrete replicator
equation is concerned with the replication-selection across the generations of
a vast well-mixed collection of the cooperators and the defectors.

We judiciously choose the above form of the replicator equation because (a) it
phenomenologically models the replication-selection dynamics that is in line
with the Darwinian tenet of the natural selection, (b) its fixed points
correspond to the Nash equilibria [58] and the evolutionarily stable
strategies/states [59] through the folk and related theorems of the
evolutionary game theory [60, 55], and (c) most importantly, it is endowed
with chaotic attractors [54, 55, 57]. Another more popular form of the
discrete replicator equation possible [61] that, however, is not conducive to
our study as it does not show [55] non-convergent dynamics for the simple two-
player–two-strategy symmetric game. While relatively less in vogue in the
biological systems, Eq. (1)—as itself or in related forms—also appears in
modelling intergenerational cultural transmission [51, 53], boundedly rational
players’ imitational behaviour in bimatrix cyclic games [50], and
reinforcement learning [49]. It is interesting to recall that a good behaviour
rule does not require aggregate population behaviour implicit in the natural
selection to induce the replicator map; instead, the map is arrived at based
on the rational behaviour of the players [52].

The basic CML considered in this work is a linear lattice with $N$ lattice
sites/nodes and periodic boundary condition such that each lattice site (or in
our case deme) is connected to its two nearest neighbours. The individuals
interact with each other within the deme via a strategic interaction modelled
by the replicator map, and there is migration between connected demes. The
dynamics on the CML is given by,

$x^{i}_{n+1}=(1-\epsilon)f(x^{i}_{n})+\frac{\epsilon}{2}[g(x^{i-1}_{n})+g(x^{i+1}_{n})],$
(3)

where, the superscript $i$ denotes the $i$th lattice site and $\epsilon$ is
the coupling strength measuring the rate of migration to the $i$th node from
the $(i-1)$th and $(i+1)$th nodes. We must restrict $\epsilon$ between $0$ to
$1$ so that $x^{i}_{n}$ does not become either negative or greater than one.
In the model under consideration, the migration can be modelled as a linear
term or, more generally, as a nonlinear term depending on the function $g(x)$.

We must remark that the coupled continuous replicator dynamics—as opposed to
the discrete one used herein—with migration included was studied [62] for a
system of two demes in which the individuals play symmetric two-player games
such that there exist two ESSs: one payoff dominant (maximal payoff state
among the equilibria) and the other risk-dominant (less risk of loss for the
players). It was concluded that due to the migration, most of the individuals
play the payoff-dominant strategy even if both the demes are initiated in the
basin of attraction of the risk-dominant equilibrium. Our system is far richer
dynamically since the discrete replicator equation can lead to nonconvergent
attractors that are clearly not connected with the game-theoretic concepts
like the evolutionary stable strategy (and hence, the Nash equilibrium).
Because of the prevalence of non-fixed point dynamics in the dynamical
systems, we must address them in spite of rather limited understanding of
their connection with the game theory. Needless to say, in the scenario where
two or more discrete replicators coupled via migration, evolutionary stable
strategies—whether payoff- or risk-dominant—are rendered unachievable in the
evolutionary dynamics since these equilibrium strategies, by construction, can
only be connected with some convergent fixed point attractors.

##  III Nonlinear coupling

Figure 1: Emergence of cooperation in the CML with only the PD, nonlinear
coupling, and no rewiring. The average cooperation $\langle\bar{x}\rangle$ at
each deme is plotted as a function of coupling strength, $\epsilon$, for the
replicator map. The black, green, blue, and red solid lines correspond to the
average cooperation as obtained from the numerical simulations done with
$\alpha=0.2$, $0.8$, $0.95,$ and $1.0$ respectively. Circular markers are
denoting the homogeneous fixed points for the respective combination of the
parameters $\alpha$ and $\epsilon$. The black, green, blue, and red dashed
lines are the theoretically predicted critical coupling strength beyond which
one gets full cooperation for $\alpha=0.2$, $0.8$, $0.95,$ and $1.0$
respectively. Here we have set $R=1.1,~{}S=0.0,~{}T=1.2,~{}P=0.1$ in the
payoff matrix of PD. Figure 2: Schematic diagram of the coupled map lattice
with dynamic random rewiring. We see in left panel (a) the base CML with eight
demes each having six representative individuals for illustrative purpose.
Every deme has a game—say, the PD or the LG—played by the individuals in it.
The arrowheads point towards the destinations of respective migration. In
right panel (b), as the dynamic random rewiring is employed, some of the
directed edges (shown by blue arrows) of the base CML are randomly broken
(shown by brown arrows with scissors) and new incoming edges (shown by green
arrows) are created.

Now, let us consider a class of nonlinear coupling—a power-law coupling with
$\alpha\geq 0$—such that the dynamics of the CML is given by,

$\displaystyle
x^{i}_{n+1}=f(x^{i}_{n})(1-\epsilon)+\frac{\epsilon}{2}\left[(x^{i+1}_{n})^{\alpha}+(x^{i-1}_{n})^{\alpha}\right],$
(4)

where, $i=1,2,\cdots,N$. This CML has two homogeneous fixed points: $x^{*}=0$
or $x^{*}=1$ at all the demes—they, respectively, correspond to all defector
or all cooperator states. The choice of power law coupling is mostly to
illustrate the possible effects of nonlinear coupling in a mathematical
tractable setting. However, a physical motivation for using the power-law
coupling may be sought in the seminal work by Zipf [63] where such a power law
in the growth of city population—due to immigration of people—was predicted.
Furthermore, in a recent work [64] based on US census data, a power law in the
migration of people from city to city is observed.

We perform the linear stability analysis [65] about a homogeneous fixed point
$x^{*}$ by substituting $x^{i}_{n}=x^{*}+h^{i}_{n}$ and expanding the
resultant equation up to the first order to get:

$h_{n+1}^{i}=\left(1-\epsilon\right)f^{\prime}\left(x^{*}\right)h_{n}^{i}+\frac{\alpha\epsilon}{2}\left(h_{n}^{i+1}+h_{n}^{i-1}\right).$
(5)

Now expressing the small perturbations as a sum of its Fourier components,
i.e., $h_{n}^{i}=\sum_{q}\tilde{h}^{q}_{n}\exp\left(\sqrt{-1}qi\right)$, where
$q$ is the wave number, and substituting in Eq. (5), we arrive at the
following expression:

$\frac{\tilde{h}^{q}_{n+1}}{\tilde{h}^{q}_{n}}=f^{\prime}\left(x^{*}\right)\left(1-\epsilon\right)+\alpha\epsilon\cos
q.$ (6)

For perturbation amplitude to decrease with time, i.e., in order for the fixed
point to be stable, we need

$\bigg{|}\frac{\tilde{h}^{q}_{n+1}}{\tilde{h}^{q}_{n}}\bigg{|}=\left|f^{\prime}\left(x^{*}\right)\left(1-\epsilon\right)+\alpha\epsilon\cos
q\right|<1.$ (7)

Therefore, in order to find the condition for establishing full cooperation in
the CML when same game is played at all the demes, we necessarily need to find
the condition for $x^{*}=1$ to be stable. In the case of $0\leq\alpha<1$, the
critical coupling strength, $\epsilon_{\rm{crit}}$ (the minimum coupling
strength required to impart full cooperation), in the CML of replicator maps
is given by:

$\epsilon_{\rm{crit}}=\frac{1}{1+\frac{1-\alpha}{T-R}},$ (8)

where we have explicitly used that fact that for the replicator map,
$f^{\prime}(1)=T-R+1$ . We have verified this result for different values of
$\alpha$ using numerical simulations as presented in Fig. 1.

From this expression (Eq. (8)), we can clearly see that if we have the CML
where at every deme the prisoner’s dilemma (PD) game ($T>R$ and $P>S$) is
played, the emergence of cooperation is possible because
$\epsilon_{\rm{crit}}$ can be less than unity. However, if $\alpha\geq 1$,
then the CML with only the PD at each deme can not lead to cooperation because
$f^{\prime}(1)=T-R+1>1$ for the PD and hence we cannot have $x^{*}=1$ stable
(see Eq. (7)).

Figure 3: Numerical validation of critical coupling strength for the CML with
dynamic rewiring. We show the order parameter, $\langle r_{G}\rangle$, as a
function of coupling strength $\epsilon$ and the rewiring probability $p$. In
subplot (a), every deme has the LG, whereas in subplot (b), only half of the
demes have the LG (payoff matrix $\tiny{\begin{pmatrix}1&7\\\ 8&0\\\
\end{pmatrix}}$) and the rest have the PD (payoff matrix
$\tiny{\begin{pmatrix}1.1&0.0\\\ 1.2&0.1\\\ \end{pmatrix}}$). The white dashed
line corresponds to Eq. (10) specifying how the critical coupling strength
varies with the rewiring probability in the CML with exclusively leader games
at all demes. The average in- and out-degree of the network is two. Here,
$N=100$ demes and the system has been evolved for $2000$ time steps. Figure 4:
Classification of twelve symmetric games. In subplot (a), we show how the
region of the $S$-$T$ parameter space is separated in twelve well-known games
whose names have been specified in the figure. The interior of the white
dashed leaf-like boundary corresponds to the region where the replicator map
has physical solutions. On taking the payoff matrices from the line $T=S+1$
and evolving the replicator map, we find period doubling route to chaos as
depicted in subplot (b) where the vertical dashed lines, from left to right,
respectively correspond to $\textsf{A}_{\rm P1}$, $\textsf{A}_{\rm P2}$,
$\textsf{A}_{\rm P4}$, and $\textsf{A}_{\rm C}$. Figure 5: Amplitude death:
The demes with the LG induce cooperation in the demes with the PD. When all
the demes of the CML are playing the LG with (a) $\textsf{A}_{\rm P1}$, (b)
$\textsf{A}_{\rm P2}$, (c) $\textsf{A}_{\rm P4}$, or (d) $\textsf{A}_{\rm C}$
after the respective critical values of the coupling strength (vertical dashed
line) all the demes’ trajectories (black, brown, light brown, and red dots
respectively for the four aforementioned LGs) synchronize onto the fixed point
$x^{*}=0.5$ of the CML. As we introduce the PD in some of the demes with no
(or few) cooperators, then the corresponding trajectories (blue dots) are
pulled onto the synchronized state $x^{*}\approx 0.5$ for all demes beyond
$\epsilon_{\rm crit}$ as exhibited in rest of the subplots for the LG game
fraction, $\phi=0.7$, $0.3$, and $0.1$. Subplots (a), (e), (i), and (m)
correspond to $\textsf{A}_{\rm P1}$; subplots (b), (f), (j), and (n)
correspond to $\textsf{A}_{\rm P2}$; subplots (c), (g), (k), and (o)
correspond to $\textsf{A}_{\rm P4}$; and subplots (d), (h), (l), and (p)
correspond to $\textsf{A}_{\rm C}$. Here, $N=100$ demes and the system has
been evolved for $2000$ time steps. Figure 6: Synchronization order parameter
and emergence of cooperation. For the CML with mixed types of demes—some
playing the PD and some the LG (with $\textsf{A}_{\rm P1}$, $\textsf{A}_{\rm
P2}$, $\textsf{A}_{\rm P4}$, and $\textsf{A}_{\rm C}$ from left to right
columns respectively)—we plot order parameter $\langle r_{G}\rangle$ and
average cooperation $\langle\bar{x}_{n}\rangle$ (averaged over 64
realizations) as functions of coupling strength $\epsilon$ and the leader game
fraction $\phi$ for the cases where the LG has period one (or fixed point)
attractor ((a) and (e)), period two attractor ((b) and (f)), period four
attractor ((c) and (g)), and chaotic attractor ((d) and (i)) when plugged into
the replicator map. The vertical black dashed line corresponds to
$\epsilon=\epsilon_{\rm crit}$ when the CML has only the LG at all the lattice
points. Here, $N=100$ demes and the system has been evolved for $2000$ time
steps.

##  IV Random Rewiring

Evidently, in the CML consisting of only one kind of demes—each having
individuals playing the prisoner’s dilemma game and linear coupling
($\alpha=1$) among them, the emergence of cooperation ($x^{*}\neq 0$) requires
dynamic rewiring as a possible mechanism [65]. In order to implement the
random migration in the system under investigation, we modify the couplings in
the CML. At every time step, any node can either allow migration from its two
nearest neighbours or two other demes picked randomly from an uniform
distribution. The probability of remaining coupled to the nearest neighbours
is $1-p$, where $p$ is called the dynamic rewiring probability; ‘dynamic’
emphasizes that the rewiring is happening at every time step.

Mathematically, the mean-field equation for the CML with the random coupling
is given by,

$\displaystyle
x^{i}_{n+1}=(1-\epsilon)f(x^{i}_{n})+\frac{(1-p)\epsilon}{2}(g(x^{i-1}_{n})+g(x^{i+1}_{n}))$
$\displaystyle+\frac{p\epsilon}{2}(g(x^{\xi}_{n})+g(x^{\eta}_{n})),$ (9)

where, $\xi$ and $\eta$ are the indices of the two randomly chosen demes and
are not equal to $i-1$, $i$, or $i+1$. Note that while the in-degree of every
node remains fixed at two, the out-degree of every node may not remain two.
See Fig. 2(b) for a schematic representation. We would like to spell out that,
in effect, there are two edges—incoming and outgoing—between two nearest
neighbours (see, e.g., Fig. 2(a)). From a deme’s perspective, an incoming edge
indicates immigration and an outgoing one denotes emigration. The parameter
$p$ in Eq. (9) can also be interpreted as a measure of how strong the long-
range migration (or non-nearest neighbour migration) is compared to the short-
range migration (nearest neighbour migration) since $\epsilon
p/\epsilon(1-p)$—the ratio of the coefficients in the third and the second
terms in the right hand side of Eq. (9)—is a monotonically increasing function
of $p$ in the range $0$ to $1$.

In the model, there is a possibility that the dynamics at all the demes may be
completely synchronized to an interior fixed point, i.e., $x^{i}=x^{j}=x^{*}$
for all $i$ and $j$. As done in the immediately preceding section (Sec. III),
one could do a linear stability analysis to find if this synchronized state is
at all stable and hence attainable. It can be shown that such a stable state
in fact exists when $\epsilon\geq\epsilon_{\rm crit}$ which can be explicitly
given as,

$\epsilon_{\rm{crit}}=\frac{|f^{\prime}(x^{*})|-1}{|f^{\prime}(x^{*})|-1+p}.$
(10)

In passing, we note that with $\alpha=1$ in Eq. (8) and with rewiring
probability $p=0$ in Eq. (10), both of them reduce to the same value of
critical coupling strength, viz, unity, as they should.

We quantify the extent to which the system is synchronized by defining a
global order parameter [66],
$r_{G}:=|\sum_{i=1}^{N}e^{2\pi\sqrt{-1}x^{i}}|/N$, that asymptotically reaches
unity as the system attains complete synchrony. For large $N$ and uniformly
distributed $x^{i}$ in the interval $[0,1]$, we can easily observe that
$r_{G}=0$. Thus, any partially synchronized state have a non-zero value of
$r_{G}$ that is less than unity. The system is completely synchronized when
$r_{G}=1$ that corresponds to the state of the population where $x^{i}=x^{j}$
for all $i,j\in\\{1,2,\cdots,N\\}$. Fig. 3(a) presents the verification that
the critical coupling strength depends on the rewiring probability $p$ exactly
as predicted by Eq. (10) for the Leader game (LG) [67, 68]—a game with chaotic
solutions [55]. What is even more satisfying is that the prediction is quite
good even when the CML has mixed kind of demes—half playing the PD and the
rest the LG (see Fig. 3(b)). In the figure, and henceforth, $\phi\in[0,1]$
denotes the fraction of the demes at which the LG is played.

###  IV.1 Twelve ordinal class of games: periodic orbits and chaos

There are twelve ordinally equivalent classes of two-player–two-strategy
symmetric normal form games [68, 55] that can be represented by the following
general payoff matrix:

${\sf A}=\begin{pmatrix}1&S\\\ T&0\end{pmatrix}.$ (11)

One must work with only those values of $S$ and $T$ for which $0\leq
x^{i}_{n}\leq 1$ for all values of $i$ and $n$. In the $S$-$T$ plane, the
twelve classes of games are clearly demarcated by the straight lines: $S=0$,
$T=0$, $S=1$, $T=1$, and $T=S$. The physical region for the replicator map
given by Eq. (1) is shown in Fig. 4(a).

Without any loss of generality, for the sake of concreteness, we choose
$S=-0.1$ and $T=1.1$ in the PD where obviously a player is better off
defecting irrespective of what the opponent chooses from the two strategies
available to her. We note that this form of the PD corresponds to the additive
model studied by W. D. Hamilton [69] and R. Trivers [70]. We have already seen
that the replicator map corresponding to the PD game displays solutions that
asymptotically reach a fixed point attractor. Also, when some of the demes
playing the PD, transform into the LG—which has chaotic solutions—the entire
CML tries to synchronize to impart cooperation in all the demes. Historically
[67], the LG’s name comes from the fact that if a player shifts its strategy
from cooperation to defection, it rewards both the players but herself more
so, and hence, is a leader. The LG has a symmetric mixed Nash equilibrium [58]
that is an unstable fixed point of the corresponding replicator map.

In order to explicitly see that the PD with altruism rewarded is the LG, we
look at the following two specific payoff matrices where subscripts indicate
that the matrices are respective that for the PD and the LG:

${\sf A}_{\rm PD}=\begin{pmatrix}1&-0.1\\\ 1.1&0\end{pmatrix}\,{\rm
and}\,\,{\sf A}_{\rm LG}=\begin{pmatrix}1&-0.1+1.15\\\ 1.1&0\end{pmatrix}.$

It is explicitly evident that on giving a cooperator playing against a
defector some extra reward, say, $+1.15$ in the above, can turn the PD game to
the LG (as now $T>S>1>0$). The LG can have different entries in the matrix as
long as the required defining condition is satisfied. It is interesting to
note that as one varies the parameters along the line $T=1+S$, starting from
$(S,T)=(-1,0)$, the replicator map undergoes period-doubling route to chaos as
shown in Fig. 4(b). Hence, it is natural to be curious about how other the
LGs, whose replicator dynamics lead to other non-chaotic asymptotic solutions,
affect the cooperators in the CML with most demes playing the PD.

To this end, we choose three LGs (marked in Fig. 4(b)) with following payoff
matrices:

${\sf A}_{\rm P1}=\begin{pmatrix}1&2\\\ 3&0\end{pmatrix},\,\,{\sf A}_{\rm
P2}=\begin{pmatrix}1&5\\\ 6&0\end{pmatrix},\,{\rm and}\,\,{\sf A}_{\rm
P4}=\begin{pmatrix}1&6.3\\\ 7.3&0\end{pmatrix},$

that respectively correspond to replicator dynamics with asymptotically stable
period one (fixed point) orbit, period-two orbit, and period-four orbit. We
set parameters $S=7$ and $T=8$ to choose another LG (also marked in Fig. 4(b))
endowed with a chaotic attractor. We denote its payoff matrix by
$\textsf{A}_{\rm C}$

###  IV.2 Mixed games in demes

We set the rewiring probability to $p=0.5$. If the CML only has an LG with
fixed payoff matrix at all the demes then synchronization ($r_{G}\approx 1$)
onto $\langle\bar{x}\rangle\approx 0.5$ occurs at threshold critical coupling
parameters that respectively are $0$, $0.5$, $0.7$, and $0.75$ for
$\textsf{A}_{\rm P1}$, $\textsf{A}_{\rm P2}$, $\textsf{A}_{\rm P4}$, and
$\textsf{A}_{\rm C}$. It should be noted that this effectively is amplitude
death because all the periodic or chaotic oscillations die when the
synchronized state is reached.

Just like the case of mixed demes with the PD and the chaotic LG (with ${\sf
A}_{\rm C}$) [15], if the chaotic LG is replaced by any of the other
aforementioned LGs, one can see that, starting from infinitesimal cooperator
fraction, enough cooperation can be established in the demes playing the PD if
the random migration is strong enough. Fig. 5 exhibits this fact in detail and
transparently. The symmetry in the subplots of the figure about the line
$x=0.5$ is easily explained. Along the line $T=1+S$ in the parameter space if
one recasts Eq. (9) using $x_{n}^{i}\rightarrow y_{n}^{i}=x_{n}^{i}-0.5$, the
dynamics in terms of $y$-variable—whose range is $[-0.5,0.5]$ for all
$i$—appears to be symmetric about zero:
$y_{n+1}=y_{n}+{S}y_{n}\left(4{y_{n}}^{2}-1\right)/2$.

The emergence and sustenance of cooperation is further depicted in Fig. 6.
Note that remarkably high degree of synchronization—viz., $\langle
r_{G}\rangle\approx 1$, angular brackets denoting average over many
realizations of random migration—is seen for the very low values of $\phi$ for
all coupling strengths and for coupling strength more than a critical value
for almost all $\phi$. The reason for the former is that the demes with the PD
dominate and defectors cannot be replaced, and the reason for the latter is
that the synchronization of the chaotic LG pulls the cooperator fraction of
all the demes with the PD up. The high values of sustained cooperation—viz,
$\langle\bar{x}\rangle\approx 0.5$, overbar denoting the average over
demes—owe to the high values of coupling strength (mostly around and beyond
$\epsilon_{\rm crit}$) for almost all $\phi$ (except for very low values).
Note that $\langle\bar{x}\rangle\approx 0.5$ at intermediate and lower values
of $\epsilon$ is accompanied with a low degree of synchronization, meaning
that at any given instant cooperator-fraction at a given deme can either be
lower or higher than $\langle\bar{x}\rangle$ and this averages out to give the
high values of cooperation.

##  V varying degree

In the CML that we have considered so far, the in-degree of every node is
always two even when the dynamic rewiring is in action. However, one could
think of a situation where the in-degree is more than two: A straightforward
generalized scenario would be when the in-degree of every node is same
$k\in\\{2,4,6\cdots\\}$ such that each node can be connected with another node
at most once and every stage of rewiring is done with probability $p$. We
conveniently choose $k$ to be an even number so that in absence of any
rewiring each node has $k$ incoming edges—one each from its nearest $k/2$
neighbours on each side. The mean field equation for node $i$ is,

$x_{n+1}^{i}=f(x_{n}^{i})(1-\epsilon)+\frac{\epsilon(1-p)}{k}\sum_{\begin{subarray}{c}l=-k/2\\\
l\neq 0\end{subarray}}^{k/2}x_{n}^{i+l}+\frac{\epsilon
p}{k}\sum_{l=1}^{k}x_{n}^{\xi_{l}}.$ (12)

$\xi_{l}$ are the randomly chosen neighbours other than the node $i$ and its
nearest $k/2$ neighbours on the either sides. For notational convenience, we
denote the $k/2$ neighbours on the clockwise direction by $i+l$ where
$l=1,2,...,k/2$ and on the anticlockwise direction by $i-l$.

Assuming that an interior fixed point exists such that at each deme $x=x^{*}$,
we find its stability by putting $x^{i}_{n+1}=x^{*}+h^{i}_{n+1}$ in Eq. (12),
and subsequently keeping only up to the linear order terms for the
perturbations. We arrive at

$h^{i}_{n+1}=(1-\epsilon)f^{\prime}(x^{*})h^{i}_{n}+\frac{\epsilon(1-p)}{k}\sum_{\begin{subarray}{c}l=-k/2\\\
l\neq 0\end{subarray}}^{k/2}h_{n}^{i+l}.$ (13)

Note that we have neglected the randomly changing neighbours’ contribution to
this equation because that should average out to zero. Again, we make the
ansatz: $h^{i}_{n}=\sum_{q}{\tilde{h}^{q}_{n}e^{\sqrt{-1}qi}}$, and Eq. (13)
yields,

$\frac{\tilde{h}^{q}_{n+1}}{\tilde{h}^{q}_{n}}=(1-\epsilon)f^{\prime}(x^{*})+\frac{2\epsilon(1-p)}{k}\sum_{l=1}^{k/2}{\rm{cos}}~{}lq.$
(14)

A stable fixed point, or in other words, a stationary synchronous state (akin
to ampplitude death in coupled oscillators), exists when the perturbations die
out, i.e., modulus of the right-hand side in Eq. (14) is less than unity.

The critical coupling strength beyond which it so happens can easily be found
numerically, as shown in Fig. 7. We remark that although the above
calculations are for the case of single kind of demes in the CML, it gives
close enough estimates even when we have mixed kinds of demes—some playing the
PD and some playing the LG as evident from Figs. 7(b) and 7(d). The most
important outcome of increasing degree is that it can overcome the necessity
of random wiring for effecting cooperation through amplitude death (contrast
the upper row with the lower row in Figs. 7). However, the effect of degree
saturates quite fast; we note that the figures are almost identical beyond
$k\approx 10$. Since we have already seen in the preceding section that all
types of LGs have effectively similar qualitative outcomes, here we have
illustrated our results only for the chaotic LG; the degree-dependece of the
other LGs are qualitatively similar as expected.

Figure 7: Amplitude death induced cooperation saturates with degree. We plot
the order parameter $\langle r_{G}\rangle$ as a function of the coupling
strength $\epsilon$ and the in-degree of the nodes of the network $k$. In
subplots (a) and (c), every deme in the network has the LG (with ${\sf A}_{\rm
C}$), while in subplots (b) and (d) only half of the nodes have the LG (the
rest have the PD). Viewed from the perspective of the rewiring probability,
subplots (a) and (b), and subplots (c) and (d) respectively correspond to
$p=0$ and $p=0.5$. The white dashed lines exhibit the analytically estimated
critical coupling strength (beyond which amplitude-death induced
synchronization is expected) as a function of the in-degree for the case where
the network has the LG exclusively. Here, $N=100$ demes and the system has
been evolved for $2000$ time steps. Figure 8: Effect of $\phi$ on
synchronization and intrademic cooperation. The order parameter $\langle
r_{G}\rangle$ ((a)-(d)) and the average intra-cooperation
$\langle\bar{x}\rangle$ ((e)-(h)) is plotted for the entire range of the
coupling strength $\epsilon$ for a fixed initial cooperating deme fraction
($C_{0}=0.1$). The four columns present the simulation results for the
rewiring probability, $p=0.25,~{}p=0.5,~{}0.75,$ and $1.0$ respectively. Blue,
red, and green colours correspond to $\phi=0.1,~{}0.5,~{}{\rm and}~{}0.9$
respectively. The dashed lines are used for the case of $C_{0}=1.0$ where
strategy update and hence, costly interdemic interactions is effectively
absent. Here, we use $N=100$ demes and evolve the system for $2000$ time
steps. All the results are averaged over $512$ realizations. The critical
coupling strength for case $\phi=1$ is shown by black dotted lines. Figure 9:
Optimal region for the coevolution of intrademic cooperation and
synchronization. The degree of coevolution of intrademic cooperation and its
synchronization $\langle\bar{x}\rangle\langle r_{G}\rangle$ is plotted with
the fraction, $\phi$, of demes where cooperation is rewarded. We fix the
coupling strength, $\epsilon=0.9$ and study for different rewiring
probability, $p=0.25,~{}0.5,~{}0.75,~{}{\rm and}~{}1.0$ respectively, in four
different columns. In each subplot, we have three different curves
corresponding to three different initial inter-cooperation levels, viz.,
$C_{0}=1.0,~{}0.1,~{}{\rm and}~{}0.01$. All the results are averaged over 512
realizations.

##  VI Strategic interdemic interaction

As mentioned in the introduction of this paper, migration is known to be
costly [42, 46, 47, 48]. The subpopulations may decide whether to cooperate by
participating in the migration or to avoid migration because the migration is
costly. Here we study the effect of costly migration on intrademic cooperation
and synchronization in our model, where cooperation is established through
random migration. As done in the case of the evolutionary Kuramoto dilemma
[18], we introduce a strategic interdemic interaction that models whether the
demes participate in migration based on their acquired payoffs. Note that this
effectively introduces an idea of cooperation at the level of demes on top of
the idea of cooperation between the replicators of any particular deme.

Mathematically, we modify Eq. (9) to incorporate the costly interdemic
interaction:

$\displaystyle
x^{i}_{n+1}=f(x_{n}^{i})(1-s^{i}_{n}\epsilon)+\frac{s^{i}_{n}\epsilon(1-p)}{2}(x^{i+1}_{n}+x^{i-1}_{n})\qquad$
$\displaystyle+\frac{s^{i}_{n}\epsilon
p}{2}(x^{\xi}_{n}+x^{\eta}_{n}).~{}~{}~{}~{}~{}~{}$ (15)

Here $s^{i}_{n}$ is the action of the $i$th deme at the $n$th time step; it
can take binary values—zero and one. From Eq. (15) we can see if $s^{i}_{n}$
is zero then the $i$th deme does not participate in migration. The demes with
strategy $s^{i}=0$ are not cooperators; they do not participate in migration
to set up a fixed synchronized level of intra-cooperation in the entire
network through amplitude death. On the other hand, demes with $s^{i}=1$ are
coupled to other demes through migration; they may be said to be cooperating
in trying to establish cooperation throughout the entire network of the
subpopulations. Therefore, while the intrademic cooperation level is
quantified by $\langle{\bar{x}}\rangle$, the interdemic cooperation level at
$n$th time step can conveniently be defined as,

$\displaystyle C_{n}\equiv\frac{1}{N}\sum_{i=1}^{N}s^{i}_{n}.$ (16)

We are interested in this section to see how $\langle{\bar{x}}\rangle$ and
$r_{G}$ depend on $C_{0}$—the initial interdemic cooperation fraction.

Of course, in order for $s_{n}^{i}$ to be a non-trivial game theoretic action,
one has to associate some payoffs for the corresponding player. In the present
context, it is straightforward because on choosing to cooperate (or
equivalently, participate in migration), the deme has to pay a cost that we
take as the rate of the deviation [18, 21] from its state:

$\displaystyle
c_{n}^{i}\equiv|[x^{i}_{n}-f(x_{n-1}^{i})]-[x^{i}_{n-1}-f(x_{n-2}^{i})]|.$
(17)

Obviously, the demes with $s^{i}=0$ incur no cost. We quantify the benefit of
each deme as a measure of how much in synchrony it is with its neighbours, and
hence, it is aptly measured by the local order parameter [18, 21] given by,

$\displaystyle b^{i}\equiv\frac{\sum_{j}r_{ij}a_{ij}}{\sum_{j}a_{ij}}.$ (18)

Here, $r_{ij}\equiv\frac{1}{2}|(e^{2\pi\iota x_{i}}+e^{2\pi\iota x_{j}})|$
(with $\iota=\sqrt{-1}$) is the pairwise order parameter. The quantity
$a_{ij}$ is adjacency matrix that takes value $1$ if $i$th and $j$th demes
interact through migration; otherwise, it is zero. Therefore the benefit is
nothing but the average pairwise order parameter over the neighbours of the
deme. In conclusion, the total payoff acquired by a deme should be given by,
$U^{i}\equiv b^{i}-\alpha c^{i}$, where $\alpha$ is the relative cost
modulating the effect of the cost.

Like in all evolutionary games, the actions $s^{i}_{n}$ must also evolve in
time in accordance with an update rule. In some of the simplest possible
update rules, any focal player (deme, in our context) would compare (in some
way) its payoff with that of its neighbours to decide whether to change its
action in the subsequent time step so as to reap more payoff. This strategy
update can be done stochastically using the myopic rationality factor $\beta$,
such that the probability of changing the strategy $s^{i}$ of $i$th deme to
the strategy of its neighbour (say, $s^{j}$ of $j$th deme) is given by
$P_{s^{i}\rightarrow s^{j}}=(1+e^{\beta(U^{i}-U^{j})})^{-1}$ [71, 72, 73]. In
the analogous deterministic rule ($\beta\rightarrow\infty$, effectively), the
focal player simply imitates the most successful strategy of its neighbours;
since the deterministic rule alters no result qualitatively, we exclusively
use the deterministic rule for the strategy update in this paper.

We simulate Eq. (15) with $N=100$ demes and employ the aforementioned
deterministic update rule. In general we work with inhomogeneous demes meaning
that $\phi$ is a variable in our simulations. For PD, we set $T=1.1$ and
$S=-0.1$, whereas we set $T=8,~{}S=7$ for the demes with the (chaotic) LG. We
start our simulations with random initial intrademic cooperation. We keep
$\alpha=0.01$ in our simulations. The initial actions $s^{i}_{0}$ for the
interdemic interaction is assigned randomly such that the initial interdemic
cooperation fraction $C_{0}$ is fixed to some pre-decided small value. We have
seen in preceding sections that the random migration helps to synchronize the
intrademic cooperation beyond a critical coupling strength $\epsilon_{\rm
crit}$ in the absence of interdemic interaction. It is clear from the strategy
update rule that the case $C_{0}=1$ effectively corresponds to the absence of
any effect of costly interdemic interaction at any point of time because the
only strategy available is to cooperate.

As presented in Fig. 8, we numerically find contrasting effects of the costly
interdemic interaction on the intrademic cooperation and the synchronization
as $\phi$ changes. We fix initial interdemic cooperation fraction as
$C_{0}=0.1$ for concreteness; and we contrast the results found in this case
with the case when the costly interdemic interaction (or in other words, the
strategy update rule) is absent (realized in the numerical simulation by
putting $C_{0}=1$).

We note that some general gross features are independent of $p$: The degree of
synchronization (see Figs. 8(a)-8(d)) is lower for higher fraction ($\phi$) of
demes with the chaotic LG, whereas the intrademic cooperation level (see Figs.
8(e)-8(h)) is higher for higher $\phi$. This can be attributed to the fact
that when interdemic cooperation level is not unity, it effectively means that
some demes are isolated, i.e., they do not participate in the interdemic
interaction (or migration); therefore, higher values of $\phi$ means
possibility of higher fraction of isolated demes with the LG that ensures some
intrademic cooperation but—being chaotic and uncoupled—asynchronous dynamics.
Moreover, the intrademic cooperation level is always lower (at any $p$,
$\phi$, or $\epsilon$) when strategic interdemic interactions are in action
than when they are absent (see Figs. 8(e)-8(h)). This is again because the
isolated demes with the LG cannot interact with the demes with the PD to
effect higher cooperator-fraction in them. For similar reasons, we also
observe that above $\epsilon_{\rm crit}(p)$ (the critical coupling strength at
$\phi=1$ and $C_{0}=1$ as a function of $p$)—at which we expect synchronized
dynamics (hence, $\langle r_{G}\rangle\rightarrow 1$) in the absence of any
interdemic strategy update rule—the degree of synchronization is lower when
strategic interdemic interactions are in action than when they are absent (see
Figs. 8(a)-8(d)).

Obviously, high values of $\phi$ are good for intrademic cooperation but not
as good for synchronization, whereas low values of $\phi$ are good for
synchronization but not good for intrademic cooperation. Thus, as hinted by
Fig. 8, we expect an optimal regime in the parameter space of $\phi$, where
both the intrademic cooperation and the synchronization are supported. We
quantify the degree of this co-evolution by the product $\langle
r_{G}\rangle\langle\bar{x}\rangle$ which should be low if either $\langle
r_{G}\rangle$ or $\langle\bar{x}\rangle$ is very low. We plot $\langle
r_{G}\rangle\langle\bar{x}\rangle$ versus $\phi$ in Fig. 9 for various initial
interdemic cooperation level, $C_{0}$, including for $C_{0}=1$ which
corresponds to the permanent absence of any effect of costly interdemic
interaction. We do note that—irrespective of the rewiring probability $p$—the
higher is the initial costly interdemic cooperation, the higher is the degree
of co-evolution.

To understand the nature of the curves in Fig. 9, we first consider the
trivial case of $C_{0}=0$ which corresponds to the absence of any interdemic
interaction at all, costly or otherwise; in fact, even at any later time, the
interdemic cooperation level remains zero because the only available strategy
is defection. Thus, for the case of mixed (isolated) demes (with the PD and
the LG) in the CML, $\langle\bar{x}\rangle$ and $\langle r_{G}\rangle$ vanish
for $\phi=0$ and $\phi=1$ respectively. With increase in $\phi$, one expects
$\langle\bar{x}\rangle$ to increase but $\langle r_{G}\rangle$ to decrease as
the number of demes with the LG increases; thus, clearly, one expects a
maximum for $\langle r_{G}\rangle\langle\bar{x}\rangle$ because $\langle
r_{G}\rangle\langle\bar{x}\rangle$ is zero at the two end points, $\phi=0$ and
$\phi=1$. With increase in the value of $C_{0}$, some interdemic connections
are established. At $\phi=0$, one still has zero intrademic cooperation as
only the PD is present. However, at $\phi=1$, because of the interdemic
couplings introduced, some synchronization ($\langle r_{G}\rangle\neq 0$) is
established beyond $\epsilon_{\rm crit}$ (this is the reason we choose high
$\epsilon$ value in Fig. 9 for illustrating our point). Consequently, $\langle
r_{G}\rangle\langle\bar{x}\rangle\neq 0$ at $\phi=1$ although $\langle
r_{G}\rangle\langle\bar{x}\rangle=0$ at $\phi=0$ since the demes with the PD
do not allow for any intrademic cooperation ($\langle\bar{x}\rangle=0$) even
in the presence of interdemic cooperation. The less is the set of isolated
demes (higher $C_{0}$), the more is the degree of synchronization; and hence,
the more is the degree of co-evolution of cooperation and synchronization. One
sees this fact validated in Fig. 9 when one compares the two curves
corresponding to $C_{0}=0.01$ and $C_{0}=0.1$. Of course, $C_{0}=1$ leaves no
demes isolated and hence establishes full synchrony for $\phi=1$, leading to
high value of intrademic cooperation as well.

In summary, the take-home message from the studies presented in this section
would be the following: The co-emergence and sustenance of the intrademic
cooperation and synchronization—effected by incentivizing the cooperators in a
few demes with the PD—is hampered in the presence of costly interdemic
interactions.

##  VII Summary

The interplay between emergence of cooperation and synchronization is an
intriguing phenomenon and a recent fertile direction of research in the
evolutionary game theory. In this paper, we have put forward our insights on
this interplay by using the coupled map lattice of replicator maps where
nonconvergent outcomes, like the periodic orbits and the chaotic orbits, are
present even in the evolutionary games with only two strategies. In the model,
the emergence of cooperation occurs as the replicator maps synchronize onto a
fixed point—a phenomenon very similar to the amplitude death in coupled
nonlinear oscillators. The migration dilemma that arises in the system under
consideration presents a well known situation—just like in other social
dilemmas such as the tragedy of commons [74], the prisoner’s dilemma [75], and
agglomeration dilemma [7]—where an individual’s interest is at odds with the
entire population’s interest as a whole. We have discussed many factors that
can lead to the aversion of the migration dilemma in the model.

The general form of the model, taking into consideration all the factors
investigated in this paper, may be summarized by the following equation:

$\displaystyle
x^{i}_{n+1}=f(x_{n}^{i})(1-s^{i}_{n}\epsilon)+\frac{s^{i}_{n}\epsilon(1-p)}{k}\sum_{\begin{subarray}{c}l=-k/2\\\
l\neq 0\end{subarray}}^{k/2}(x^{i+l}_{n})^{\alpha}\qquad$
$\displaystyle+\frac{s^{i}_{n}\epsilon
p}{k}\sum_{l=1}^{k}(x_{n}^{\xi_{l}})^{\alpha};\,\,(i=1,2,\cdots,N).\qquad$
(19)

There are five main parameters in our model: (i) the coupling strength
$\epsilon$, (ii) the rewiring probability $p$, (iii) the strategy $s^{i}_{n}$
of the $i^{th}$ node at the $n^{th}$ time step, (iv) the degree of the network
$k$, and (v) the exponent $\alpha$ that makes the migration nonlinear. The
model used in Sec. III for nonlinear coupling is obtained by fixing $p=0$,
$k=2$, and $s^{i}_{n}=1\textrm{ }\forall i,n$ in Eq. (VII). The model used in
Sec. IV is obtained by putting $\alpha=1$, $k=2$, and $s^{i}_{n}=1\textrm{
}\forall i,n$ in Eq. (VII). Similarly, in Sec. V, the degree $k$ of the
network is greater than 2; and $\alpha=1$ and $s^{i}_{n}=1\textrm{ }\forall
i,n$. Finally, we have introduced costly migration between the demes in Sec.
VI where we have calculated $s^{i}_{n}$ for each node $i$ at each time step
$n$ with $k=2$ and $\alpha=1$.

We have shown that the migration modelled with a nonlinear power-law function
can impart full cooperation even if the whole population is playing only the
PD. Moreover, when the migration is modelled with a linear function, a little
initiative in rewarding altruism in a few of the demes—thus, transforming the
PD to the LG—goes a long way in raising the cooperator-fraction of the
population. The network’s degree also plays an essential role as we do not
require random migration to establish cooperation across the network when the
network’s degree is greater than two. If the degree is high enough, then the
critical coupling strength for synchronization becomes effectively independent
of the network’s degree. We have also presented in detail how the interdemic
interaction due to the costly migration affects the order of synchrony quite
nontrivially. The extent of the simultaneous emergence of cooperation and
synchrony depends on the fraction of demes where altruism is rewarded.
Interestingly, the magnitude of this simultaneous emergence of cooperation and
synchrony is maximum at a finite value (not unity) of the fraction $\phi$. It
is also important to note that the interdemic cooperation is required for
cooperation to be sustained in the individual subpopulations. Such a
multilevel interaction can be thought of as a potential model for studying
group selection [76, 77, 78, 79].

##  VIII Discussion

The random migration, which has been crucial in arriving at the results in
this paper, may be seen in a different light: For the individuals in a
particular deme when migrating to another deme, there is a chance to get
involved in a different game than what they have been playing before
migration. Especially, for the payoff matrices chosen in this paper, a player
previously playing the PD when migrates to a deme playing the LG, it gets a
chance to reap more payoff. Thus, effectively, the CML with random migration
is somewhat akin to a stochastic game [80, 81]; however, it is more intriguing
as the fraction of players jumping into a neighbouring deme with different
payoff matrix is determined by a nonconvergent deterministic dynamics. In
passing, we mention that the fact that the random migration helps in
establishing a higher cooperation state is interestingly at odds with what
happens in the case of demes with a finite population where an increase in the
migration rate, for the same benefit-to-cost ratio, decreases the probability
of fixation by the cooperators [79].

We would also like to contrast our migration-based model with other diffusion-
based models in the literature. In a model [4] with asymmetric diffusion rate
for the cooperators and the defectors in a population with the ecological
public goods games played in subpopulations leads to pattern formations in the
densities of the cooperators and the defectors. When a third strategy,
loner—agent who does not participate in public goods game—is incorporated [82]
in the optional public goods game along with the two usual strategies (viz.
cooperate and defect), the frequencies of the strategies undergo global
oscillations. In contrast, we use a single (in effect, symmetric) diffusion
coefficient—with the agents in the subpopulations playing two-player games (LG
or PD)—to bring about the coexistence of the cooperators and the defectors
over the entire population homogeneously. Regardless of the differences, it
would be interesting to incorporate multi-player games like the public goods
games and also multi-strategy games [83] within subpopulations in our model.

Our model is technically a multilayer network [84, 85, 86, 87] where inside
each node the individuals have an undirected all-to-all connection being
modelled by the one-shot games, and between two nodes there is directed
connection fashioned by the migration. Moreover, since the edges between any
two nodes are stochastic functions of time, one could consider our directed
multilayer network as temporal. Consequently, the emergence of cooperation in
it could be seen as an extension of similar enhancement of cooperation in
temporal networks [88]. A possible extension of our work could be to
investigate how different multilayered network structures affect the emergence
of cooperation in multiplayer games with many strategies.

The random rewiring scheme used in the CML could remind the readers about
similar works done with the small-world networks. Previous works on small-
world networks, with agents placed on the nodes and playing the PD game, focus
on the role of topology [89, 90], introduction of a new strategy [91], and
aspiration level of agents [5] in promoting cooperation in the network. In
contrast, we focus on building the simplest possible model to showcase the
role of non-convergent dynamics and synchronization in bringing about
cooperation in the entire network. Another key difference is that we have
placed subpopulations (and not individual agents) on the nodes of the network.
Placing individual agents on nodes of a small world network yields a lower
fraction of cooperators for a higher rewiring probability [89, 90]. This is at
odds with our results, where the system attains full cooperation beyond the
critical coupling strength regardless of the rewiring probability when the
power law exponent $\alpha$ is less than 1; and for $\alpha=1$, one can have
stable coexistence of cooperation and defectors, and the stability depends on
rewiring probability $p$.

We emphasize once again that the replicator map used in this paper does not
consider the network structure within the subpopulation, i.e., it is assumed
that an agent interacts with all other agents within a subpopulation.
Interestingly, a formalism of replicator equation on graphs has been proposed
[92] where the underlying structure of the population is taken into
consideration. It would be interesting to extend our model by considering
structured subpopulations in line with this formalism. Also, we realize that
the replicator map being a deterministic equation, should be seen as the mean-
field limit of a microscopic stochastic process [93]; in a finite population,
stochastic effects must have impact on the evolution of cooperation and its
evolutionary stability [94, 95]. A stochastic differential equation approach
[96] was adopted to deal with the finiteness of the deme size where the
direction and the strength of selection varied from deme to deme. Naturally,
in future, one may envisage investigating the stochastic phenomena, like
fixation probability, in the CML with finite subpopulations.

Besides the pairwise interaction, multiplayer interactions where more than two
individuals are involved is ubiquitous and very relevant in the study of the
evolution of cooperation [97, 98, 99, 100, 101]. The multiplayer interactions
can not be reduced to the sum of pairwise interactions that makes it quite
fascinating [95, 102]. Naturally, non-linear fitness gets introduced through
multiplayer interactions. Several studies on the effect of interacting group
size have been performed and they highlight the detrimental impact on the
evolution of cooperation [97, 98]. Realistic interactions are spatially
structured and are conveniently modelled by means of networks without all-to-
all interactions; as could be guessed, their topologies have inevitable
effects on the corresponding results. On a square lattice, enhanced spatial
reciprocity can additionally promote cooperation in large groups [103]. Our
study can be extended for the multiplayer interactions by introducing the
multiplayer replicator maps using public goods game or multiplayer PD [104,
105, 4, 106]. The $d$-player $n$-strategy interaction gives rise to
$(d-1)^{n-1}$ number of internal equilibria at maximum for continuous
replicator dynamics [99, 101]. The distribution and the number of stable
equilibria was investigated for $d$-player two-strategy game with random
payoff matrix; it was found that the number of stable equilibria is
$\sqrt{d-1}$ asymptotically [107]. It is easy to see that the discrete
replicator map can have a maximum of $(d-1)^{n-1}$ internal equilibria if it
is extended for $d$-player $n$-strategy game. Thus, there is a possibility
that with only two strategies, multiplayer PD can give rise to internal
equilibria which imply the coexistence of cooperation and defection, and this
could trigger even higher degree of cooperation in the CML with LG in some
demes.

Before we conclude, we would also like to point out that some variants of the
coupled replicator equations [108, 109, 110, 111, 112] model the dynamics of
collective learning in a group of agents who may be faced with the social
dilemmas [113]. They show quasiperiodicity, limit cycles, intermittency, and
chaos, among other rich dynamical features. In this setting, one should be
able to explore collective phenomena like synchronization onto an attractor.
Furthermore, one could also ponder upon how bounded rationality [114, 115,
116], mutation [117] and delay [118], and hypergames [119] in a deme would
modify the results of this paper.

## References

  * Kimura and Weiss [1964] M. Kimura and G. H. Weiss, Genetics 49, 561 (1964).
  * Latter [1973] B. D. H. Latter, Genetics 73, 147 (1973).
  * Helbing and Yu [2009] D. Helbing and W. Yu, Proc. Natl. Acad. Sci. U.S.A. 106, 3680 (2009).
  * Wakano  _et al._ [2009] J. Y. Wakano, M. A. Nowak, and C. Hauert, Proc. Natl. Acad. Sci. U.S.A. 106, 7910 (2009).
  * Yang  _et al._ [2010] H.-X. Yang, Z.-X. Wu, and B.-H. Wang, Phys. Rev. E 81, 065101 (2010).
  * Yu [2011] W. Yu, Phys. Rev. E 83, 026105 (2011).
  * Roca and Helbing [2011] C. P. Roca and D. Helbing, Proc. Natl. Acad. Sci. U.S.A. 108, 11370 (2011).
  * Fu and Nowak [2012] F. Fu and M. A. Nowak, J. Stat. Phys. 151, 637 (2012).
  * Wu  _et al._ [2012] T. Wu, F. Fu, Y. Zhang, and L. Wang, Phys. Rev. E 85, 066104 (2012).
  * Chen  _et al._ [2012] X. Chen, A. Szolnoki, and M. Perc, Phys. Rev. E 86, 036101 (2012).
  * Buesser  _et al._ [2013] P. Buesser, M. Tomassini, and A. Antonioni, Phys. Rev. E 88, 042806 (2013).
  * Limdi  _et al._ [2018] A. Limdi, A. Pérez-Escudero, A. Li, and J. Gore, Nat. Commun. 9, 2969 (2018).
  * Dugatkin and Wilson [1991] L. A. Dugatkin and D. S. Wilson, The American Naturalist 138, 687 (1991).
  * Enquist and Leimar [1993] M. Enquist and O. Leimar, Animal Behaviour 45, 747 (1993).
  * Sadhukhan  _et al._ [2021] S. Sadhukhan, R. Chattopadhyay, and S. Chakraborty, Phys. Rev. Research 3, 013009 (2021).
  * Kaneko [1992] K. Kaneko, Chaos 2, 279 (1992).
  * Kaneko and Yanagita [2014] K. Kaneko and T. Yanagita, Scholarpedia 9, 4085 (2014), revision #149460.
  * Antonioni and Cardillo [2017] A. Antonioni and A. Cardillo, Phys. Rev. Lett. 118, 238301 (2017).
  * Liu  _et al._ [2018] X.-S. Liu, Z.-X. Wu, and J.-Y. Guan, EPL 122, 20001 (2018).
  * Yang  _et al._ [2018] H.-X. Yang, T. Zhou, and Z.-X. Wu, Phys. Rev. E 98, 022201 (2018).
  * Chattopadhyay  _et al._ [2020] R. Chattopadhyay, S. Sadhukhan, and S. Chakraborty, Chaos 30, 113111 (2020).
  * Saxena  _et al._ [2012] G. Saxena, A. Prasad, and R. Ramaswamy, Phys. Rep. 521, 205 (2012), amplitude Death: The Emergence of Stationarity in Coupled Nonlinear Systems.
  * Jost and Joy [2001] J. Jost and M. P. Joy, Phys. Rev. E 65, 016201 (2001).
  * Martí and Masoller [2003] A. C. Martí and C. Masoller, Phys. Rev. E 67, 056219 (2003).
  * Lind  _et al._ [2004] P. G. Lind, J. A. C. Gallas, and H. J. Herrmann, Phys. Rev. E 70, 056207 (2004).
  * Jalan  _et al._ [2005] S. Jalan, R. E. Amritkar, and C.-K. Hu, Phys. Rev. E 72, 016211 (2005).
  * Amritkar  _et al._ [2005] R. E. Amritkar, S. Jalan, and C.-K. Hu, Phys. Rev. E 72, 016212 (2005).
  * C.  _et al._ [2008] M. P. C., C. Masoller, and A. C. Martí, The European Physical Journal B 67, 83 (2008).
  * Li  _et al._ [2009] C. Li, C. Xu, W. Sun, J. Xu, and J. Kurths, Chaos: An Interdisciplinary Journal of Nonlinear Science 19, 013106 (2009).
  * Turing [1952] A. Turing, Philosophical Transactions of the Royal Society of London. Series B, Biological Sciences 237, 37 (1952).
  * Prigogine and Lefever [1968] I. Prigogine and R. Lefever, The Journal of Chemical Physics 48, 1695 (1968).
  * Bar-Eli [1984] K. Bar-Eli, The Journal of Physical Chemistry 88, 3616 (1984).
  * Ramana Reddy  _et al._ [1998] D. V. Ramana Reddy, A. Sen, and G. L. Johnston, Phys. Rev. Lett. 80, 5109 (1998).
  * Karnatak  _et al._ [2010] R. Karnatak, N. Punetha, A. Prasad, and R. Ramaswamy, Phys. Rev. E 82, 046219 (2010).
  * Karnatak  _et al._ [2007] R. Karnatak, R. Ramaswamy, and A. Prasad, Phys. Rev. E 76, 035201 (2007).
  * Atay [2003] F. M. Atay, Phys. Rev. Lett. 91, 094101 (2003).
  * Aronson  _et al._ [1990] D. Aronson, G. Ermentrout, and N. Kopell, Physica D: Nonlinear Phenomena 41, 403 (1990).
  * Konishi [2003] K. Konishi, Phys. Rev. E 68, 067202 (2003).
  * Prasad  _et al._ [2010a] A. Prasad, M. Dhamala, B. M. Adhikari, and R. Ramaswamy, Phys. Rev. E 81, 027201 (2010a).
  * Prasad  _et al._ [2003] A. Prasad, Y.-C. Lai, A. Gavrielides, and V. Kovanis, Physics Letters A 318, 71 (2003).
  * Prasad  _et al._ [2010b] A. Prasad, M. Dhamala, B. M. Adhikari, and R. Ramaswamy, Phys. Rev. E 82, 027201 (2010b).
  * Rankin and Burchsted [1992] M. A. Rankin and J. C. A. Burchsted, Annual Review of Entomology 37, 533 (1992).
  * Ian [2007] N. Ian,  _The Migration Ecology of Birds_ (Elsevier, 2007).
  * Harden Jones [1968] F. R. Harden Jones,  _Fish migration_ (London : Edward Arnold, 1968) bibliography: p. 275-303.
  * Harris  _et al._ [2009] G. Harris, S. Thirgood, J. Hopcraft, J. Cromsight, and J. Berger, Endangered Species Research 7, 55 (2009).
  * Flack  _et al._ [2016] A. Flack, W. Fiedler, J. Blas, I. Pokrovsky, M. Kaatz, M. Mitropolsky, K. Aghababyan, I. Fakriadis, E. Makrigianni, L. Jerzak, H. Azafzaf, C. Feltrup-Azafzaf, S. Rotics, T. M. Mokotjomela, R. Nathan, and M. Wikelski, Science Advances 2, e1500931 (2016).
  * Lok  _et al._ [2015] T. Lok, O. Overdijk, and T. Piersma, Biology Letters 11, 20140944 (2015).
  * Lok  _et al._ [2017] T. Lok, L. Veldhoen, O. Overdijk, J. M. Tinbergen, and T. Piersma, Journal of Animal Ecology 86, 998 (2017).
  * Börgers and Sarin [1997] T. Börgers and R. Sarin, J. Econ. Theory 77, 1 (1997).
  * Hofbauer and Schlag [2000] J. Hofbauer and K. H. Schlag, J. Evol. Econ. 10, 523 (2000).
  * Bisin and Verdier [2000] A. Bisin and T. Verdier, Q. J. Econ 115, 955 (2000).
  * Cressman [2003] R. Cressman,  _Evolutionary dynamics and extensive form games_ (MIT Press, Cambridge, MA., 2003).
  * Montgomery [2010] J. D. Montgomery, Am. Econ. J. Microecon. 2, 115 (2010).
  * Vilone  _et al._ [2011] D. Vilone, A. Robledo, and A. Sánchez, Phys. Rev. Lett. 107, 038101 (2011).
  * Pandit  _et al._ [2018] V. Pandit, A. Mukhopadhyay, and S. Chakraborty, Chaos 28, 033104 (2018).
  * Mukhopadhyay and Chakraborty [2020a] A. Mukhopadhyay and S. Chakraborty, J. Theor. Biol. 497, 110288 (2020a).
  * Mukhopadhyay and Chakraborty [2020b] A. Mukhopadhyay and S. Chakraborty, Chaos 30, 121104 (2020b).
  * Nash [1950] J. F. Nash, Proc. Natl. Acad. Sci. U.S.A. 36, 48 (1950).
  * Smith and Price [1973] J. M. Smith and G. R. Price, Nature 246, 15 (1973).
  * Cressman and Tao [2014] R. Cressman and Y. Tao, Proc. Natl. Acad. Sci. U.S.A. 111, 10810 (2014).
  * Taylor and Jonker [1978] P. D. Taylor and L. B. Jonker, Math. Biosci. 40, 145 (1978).
  * Miȩkisz and Płatkowski [2005] J. Miȩkisz and T. Płatkowski, J. Theor. Biol. 237, 363 (2005).
  * Zipf [1949] G. Zipf,  _Human Behavior and the Principle of Least Effort: An Introduction to Human Ecology_ (Addison-Wesley Press, 1949).
  * Curiel  _et al._ [2018] R. P. Curiel, L. Pappalardo, L. Gabrielli, and S. R. Bishop, PLoS One 13, e0199892 (2018).
  * Sinha [2002] S. Sinha, Phys. Rev. E 66, 016209 (2002).
  * Wang  _et al._ [2016] S.-J. Wang, R.-H. Du, T. Jin, X.-S. Wu, and S.-X. Qu, Sci. Rep. 6, 23448 (2016).
  * Rapoport [1967] A. Rapoport, Behav. Sci. 12, 81 (1967).
  * Hummert  _et al._ [2014] S. Hummert, K. Bohl, D. Basanta, A. Deutsch, S. Werner, G. Theißen, A. Schroeter, and S. Schuster, Mol. BioSyst. 10, 3044 (2014).
  * Hamilton [1964] W. Hamilton, J. Theor. Biol. 7, 1 (1964).
  * Trivers [1971] R. L. Trivers, Q. Rev. Biol. 46, 35 (1971).
  * Blume [1993] L. E. Blume, Games and Economic Behavior 5, 387 (1993).
  * Szabó and Tőke [1998] G. Szabó and C. Tőke, Phys. Rev. E 58, 69 (1998).
  * Roca  _et al._ [2009] C. P. Roca, J. A. Cuesta, and A. Sánchez, Physics of Life Reviews 6, 208 (2009).
  * Hardin [1968] G. Hardin, Science 162, 1243 (1968).
  * Rapoport and Chammah [1965] A. Rapoport and A. Chammah,  _Prisoner’s Dilemma_ (University of Michigan Press, Ann Arbor, MI., 1965).
  * Aoki [1982] K. Aoki, Evolution 36, 832 (1982).
  * Crow and Aoki [1982] J. F. Crow and K. Aoki, Proc. Natl. Acad. Sci. U.S.A. 79, 2628 (1982).
  * Rogers [1990] A. R. Rogers, Am. Nat. 135, 398 (1990).
  * Traulsen and Nowak [2006] A. Traulsen and M. A. Nowak, Proc. Natl. Acad. Sci. U.S.A. 103, 10952 (2006).
  * Hilbe  _et al._ [2018] C. Hilbe, Š. Šimsa, K. Chatterjee, and M. A. Nowak, Nature 559, 246 (2018).
  * Su  _et al._ [2019] Q. Su, A. McAvoy, L. Wang, and M. A. Nowak, Proc. Natl. Acad. Sci. U.S.A. 116, 25398 (2019).
  * Valverde  _et al._ [2017] P. A. Valverde, R. da Silva, and E. V. Stock, Physica A 474, 61 (2017).
  * Chowdhury  _et al._ [2021] S. N. Chowdhury, S. Kundu, M. Perc, and D. Ghosh, Proc. Math. Phys. Eng. Sci. 477, 20210397 (2021).
  * Boccaletti  _et al._ [2014] S. Boccaletti, G. Bianconi, R. Criado, C. del Genio, J. Gómez-Gardeñes, M. Romance, I. Sendiña-Nadal, Z. Wang, and M. Zanin, Phys. Rep. 544, 1 (2014).
  * Wang  _et al._ [2015] Z. Wang, L. Wang, A. Szolnoki, and M. Perc, Eur. Phys. J. B 88, 124 (2015).
  * Iranzo  _et al._ [2016] J. Iranzo, J. M. Buldú, and J. Aguirre, Nat. Commun. 7, 13273 (2016).
  * Requejo and Díaz-Guilera [2016] R. J. Requejo and A. Díaz-Guilera, Phys. Rev. E 94, 022301 (2016).
  * Li  _et al._ [2020] A. Li, L. Zhou, Q. Su, S. P. Cornelius, Y.-Y. Liu, L. Wang, and S. A. Levin, Nat. Commun. 11, 2259 (2020).
  * Abramson and Kuperman [2001] G. Abramson and M. Kuperman, Phys. Rev. E 63, 030901 (2001).
  * Masuda and Aihara [2003] N. Masuda and K. Aihara, Phys. lett. A 313, 55 (2003).
  * Wu  _et al._ [2005] Z.-X. Wu, X.-J. Xu, Y. Chen, and Y.-H. Wang, Phys. Rev. E 71, 037103 (2005).
  * Ohtsuki and Nowak [2006] H. Ohtsuki and M. A. Nowak, J. Theor. Biol. 243, 86 (2006).
  * Mukhopadhyay and Chakraborty [2021] A. Mukhopadhyay and S. Chakraborty, Chaos 31, 023123 (2021).
  * Nowak  _et al._ [2004] M. A. Nowak, A. Sasaki, C. Taylor, and D. Fudenberg, Nature 428, 646 (2004).
  * Wu  _et al._ [2013] B. Wu, A. Traulsen, and C. Gokhale, Games 4, 182 (2013).
  * Constable and McKane [2014] G. W. Constable and A. J. McKane, J. Theor. Biol. 358, 149 (2014).
  * Joshi [1987] N. V. Joshi, Journal of Genetics 66, 69 (1987).
  * Boyd and Richerson [1988] R. Boyd and P. J. Richerson, J. Theor. Biol. 132, 337 (1988).
  * Gokhale and Traulsen [2010] C. S. Gokhale and A. Traulsen, Proc. Natl. Acad. Sci. U.S.A. 107, 5500 (2010).
  * Han  _et al._ [2012] T. A. Han, A. Traulsen, and C. S. Gokhale, Theor Popul Biol 81, 264 (2012).
  * Gokhale and Traulsen [2014] C. S. Gokhale and A. Traulsen, Dyn Games Appl 4, 468 (2014).
  * Perc  _et al._ [2013] M. Perc, J. Gómez-Gardeñes, A. Szolnoki, L. M. Floría, and Y. Moreno, J R Soc Interface 10, 20120997 (2013).
  * Szolnoki and Perc [2011] A. Szolnoki and M. Perc, Phys. Rev. E 84 (2011).
  * Kurokawa  _et al._ [2010] S. Kurokawa, J. Y. Wakano, and Y. Ihara, Theor Popul Biol 77, 257 (2010).
  * Kurokawa [2018] S. Kurokawa, Games 9, 1 (2018).
  * Wang  _et al._ [2012] B. Wang, X. Chen, and L. Wang, J. Stat. Mech. Theory Exp. 2012, P11017 (2012).
  * Duong and Han [2016] M. H. Duong and T. A. Han, J Math Biol 73, 1727 (2016).
  * Sato and Crutchfield [2003] Y. Sato and J. P. Crutchfield, Phys. Rev. E 67, 015206 (2003).
  * Traulsen  _et al._ [2004] A. Traulsen, T. Röhl, and H. G. Schuster, Phys. Rev. Lett. 93, 028701 (2004).
  * Galla [2009] T. Galla, Phys. Rev. Lett. 103, 198702 (2009).
  * Kianercy and Galstyan [2013] A. Kianercy and A. Galstyan, Phys. Rev. E 88, 012815 (2013).
  * Jost and Li [2014] J. Jost and W. Li, Phys. Rev. E 89, 022113 (2014).
  * Macy and Flache [2002] M. W. Macy and A. Flache, Proc. Natl. Acad. Sci. U.S.A. 99, 7229 (2002).
  * Simon [1955] H. A. Simon, Q. J. Econ. 69, 99 (1955).
  * Rubinstein [1998] A. Rubinstein,  _Modeling bounded rationality_ (MIT Press, Cambridge, MA., 1998).
  * Sethi [2000] R. Sethi, Games Econ. Behav. 32, 85 (2000).
  * Toupo and Strogatz [2015] D. F. P. Toupo and S. H. Strogatz, Phys. Rev. E 91, 052907 (2015).
  * Mittal  _et al._ [2020] S. Mittal, A. Mukhopadhyay, and S. Chakraborty, Phys. Rev. E 101, 042410 (2020).
  * Jiang  _et al._ [2018] J. Jiang, Y.-Z. Chen, Z.-G. Huang, and Y.-C. Lai, Phys. Rev. E 98, 042305 (2018).