$\left.\frac{\mathrm{d}\dot{x}_{2}}{\mathrm{d}x_{2}}\right|_{\mathbf{x}=\mathbf{x}^{(CE)}}=-\left(\frac{rc}{k+1}-c\right)-kx_{3}^{(CE)}\beta.$
(S66)

Therefore, $\mathbf{x}^{(CE)}$ is stable if $x_{1}^{(CE)}<x_{1,\star}^{(CE)}$
(or $x_{3}^{(CE)}>x_{3,\star}^{(CE)}$), where

$x_{1,\star}^{(CE)}=1+\left(\frac{r}{k+1}-1\right)\frac{c}{k\beta}\equiv\left(1+\frac{\alpha}{\beta}\right)x_{2}^{(DE)},$
(S67)

or $x_{3,\star}^{(CE)}=\left(-\frac{r}{k+1}+1\right)\frac{c}{k\beta}$. That
is, while every point on the $CE$-edge is equilibrium, only the points on the
side of $x_{1}^{(CE)}<x_{1,\star}^{(CE)}$ are stable. Moreover, we note that
$x_{1,\star}^{(CE)}>1$ when $r>k+1$. That is,
$x_{1}^{(CE)}<x_{1,\star}^{(CE)}$ always holds and the points on the $CE$-edge
are stable everywhere for $r>k+1$.

####  S3.2.2 The structured population

We notice that the payoff structure given by Eq. (S54) is linear, which means
that we can utilize the simplified method for special linear systems given by
section S2.4.3 for convenience.

Comparing the payoff structure $a_{i|\mathbf{k}}$ in Eq. (S54) with Eqs.
(S29)–(S31), we extract matrices $\mathbf{b}$ and $\mathbf{c}$,

$\mathbf{b}=\begin{pmatrix}\dfrac{rc}{k+1}&0&\dfrac{rc}{k+1}\\\\[10.00002pt]
\dfrac{rc}{k+1}&0&\dfrac{rc}{k+1}-\beta\\\\[10.00002pt]
\dfrac{rc}{k+1}&-\alpha&\dfrac{rc}{k+1}\end{pmatrix},\leavevmode\nobreak\
\mathbf{c}=\begin{pmatrix}\dfrac{rc}{k+1}-c\\\\[10.00002pt] 0\\\
\dfrac{rc}{k+1}-c\end{pmatrix}.$ (S68)

According to Eq. (S35), let us calculate

$\displaystyle 3\sum_{j=1}^{3}x_{j}(b_{11}-b_{1j}-b_{j1}-b_{jj})$
$\displaystyle=-6(x_{1}+x_{3})\frac{rc}{k+1},$ (S69a) $\displaystyle
3\sum_{j=1}^{3}x_{j}(b_{22}-b_{2j}-b_{j2}-b_{jj})$
$\displaystyle=-6(x_{1}+x_{3})\frac{rc}{k+1}+3x_{3}(\alpha+\beta),$ (S69b)
$\displaystyle 3\sum_{j=1}^{3}x_{j}(b_{33}-b_{3j}-b_{j3}-b_{jj})$
$\displaystyle=-6(x_{1}+x_{3})\frac{rc}{k+1}+3x_{2}(\alpha+\beta),$ (S69c)
$\displaystyle 6\sum_{j=1}^{3}\sum_{l=1}^{3}x_{j}x_{l}b_{jl}$
$\displaystyle=6(x_{1}+x_{3})\frac{rc}{k+1}-6x_{2}x_{3}(\alpha+\beta).$ (S69d)

Then, we can substitute them into Eq. (S35). Meanwhile, the expression of
$\bar{\pi}_{i}$ and $\bar{\pi}$ is the same as we give in section S3.2.1 for
well-mixed populations. In this way, we obtain the replicator equations in
structured populations as follows.

$\displaystyle\dot{x}_{1}$
$\displaystyle=\frac{\delta(k-2)}{2(k-1)}x_{1}\left\\{(k+1)\left[(1-x_{1}-x_{3})\left(\frac{rc}{k+1}-c\right)+kx_{2}x_{3}(\alpha+\beta)\right]-6x_{2}x_{3}(\alpha+\beta)\right\\},$
(S70a) $\displaystyle\dot{x}_{2}$
$\displaystyle=\frac{\delta(k-2)}{2(k-1)}x_{2}\left\\{(k+1)\left[-(x_{1}+x_{3})\left(\frac{rc}{k+1}-c\right)+k\left(x_{2}x_{3}(\alpha+\beta)-x_{3}\beta\right)\right]-6x_{2}x_{3}(\alpha+\beta)+3x_{3}(\alpha+\beta)\right\\},$
(S70b) $\displaystyle\dot{x}_{3}$
$\displaystyle=\frac{\delta(k-2)}{2(k-1)}x_{3}\left\\{(k+1)\left[(1-x_{1}-x_{3})\left(\frac{rc}{k+1}-c\right)+k\left(x_{2}x_{3}(\alpha+\beta)-x_{2}\alpha\right)\right]-6x_{2}x_{3}(\alpha+\beta)+3x_{2}(\alpha+\beta)\right\\}.$
(S70c)

Similarly, solving $\dot{\mathbf{x}}=\mathbf{0}$, we obtain equilibrium
points, which can be divided into three categories. The first and second
categories are single equilibrium points: a point on the $D$-vertex,
$\mathbf{x}^{(D)}=(0,1,0)$, and a point on the $DE$-edge,
$\mathbf{x}^{(DE)}=(0,x_{2}^{(DE)},x_{3}^{(DE)})$, where

$\displaystyle x_{2}^{(DE)}$
$\displaystyle=\frac{1}{(k-2)(k+3)(\alpha+\beta)}\left[(k+1)\left(\frac{rc}{k+1}-c+k\beta\right)-3(\alpha+\beta)\right],$
(S71a) $\displaystyle x_{3}^{(DE)}$
$\displaystyle=1-x_{2}^{(DE)}=\frac{1}{(k-2)(k+3)(\alpha+\beta)}\left[(k+1)\left(-\frac{rc}{k+1}+c+k\alpha\right)-3(\alpha+\beta)\right].$
(S71b)

The third category contains infinite equilibrium points on the $CE$-edge,
denoted by $\mathbf{x}^{(CE)}=(x_{1}^{(CE)},0,x_{3}^{(CE)})$, where $0\leq
x_{1}^{(CE)}\leq 1$, $0\leq x_{3}^{(CE)}\leq 1$,
$x_{1}^{(CE)}+x_{3}^{(CE)}=1$. This category covers other two vertex points
$(1,0,0)$ and $(0,0,1)$.

Again, the stability of $\mathbf{x}^{(D)}$ and $\mathbf{x}^{(DE)}$ can be
studied by the regular method. We cancel $x_{2}=1-x_{1}-x_{3}$ and study the
dynamics depicted by $\dot{x}_{1}$ and $\dot{x}_{3}$,

$\displaystyle\dot{x}_{1}$
$\displaystyle=\frac{\delta(k-2)}{2(k-1)}x_{1}(1-x_{1}-x_{3})\left\\{(k+1)\left(\frac{rc}{k+1}-c\right)+x_{3}(k-2)(k+3)(\alpha+\beta)\right\\},$
(S72a) $\displaystyle\dot{x}_{3}$
$\displaystyle=\frac{\delta(k-2)}{2(k-1)}x_{3}(1-x_{1}-x_{3})\left\\{(k+1)\left(\frac{rc}{k+1}-c-k\alpha\right)+x_{3}(k-2)(k+3)(\alpha+\beta)+3(\alpha+\beta)\right\\}.$
(S72b)

The Jacobian matrix of system (S72) is

$J=\begin{pmatrix}\displaystyle{\frac{\partial\dot{x}_{1}}{\partial
x_{1}}}&\displaystyle{\frac{\partial\dot{x}_{1}}{\partial
x_{3}}}\\\\[10.00002pt] \displaystyle{\frac{\partial\dot{x}_{3}}{\partial
x_{1}}}&\displaystyle{\frac{\partial\dot{x}_{3}}{\partial
x_{3}}}\end{pmatrix},$ (S73)

where

$\displaystyle\frac{\partial\dot{x}_{1}}{\partial x_{1}}$
$\displaystyle=\frac{\delta(k-2)}{2(k-1)}(1-2x_{1}-x_{3})\left\\{(k+1)\left(\frac{rc}{k+1}-c\right)+x_{3}(k-2)(k+3)(\alpha+\beta)\right\\},$
(S74) $\displaystyle\frac{\partial\dot{x}_{1}}{\partial x_{3}}$
$\displaystyle=-\frac{\delta(k-2)}{2(k-1)}x_{1}\left\\{(k+1)\left(\frac{rc}{k+1}-c\right)-(1-x_{1}-2x_{3})(k-2)(k+3)(\alpha+\beta)\right\\},$
(S75) $\displaystyle\frac{\partial\dot{x}_{3}}{\partial x_{1}}$
$\displaystyle=-\frac{\delta(k-2)}{2(k-1)}x_{3}\left\\{(k+1)\left(\frac{rc}{k+1}-c-k\alpha\right)+x_{3}(k-2)(k+3)(\alpha+\beta)+3(\alpha+\beta)\right\\},$
(S76) $\displaystyle\frac{\partial\dot{x}_{3}}{\partial x_{3}}$
$\displaystyle=\frac{\delta(k-2)}{2(k-1)}\left\\{(1-x_{1}-2x_{3})\left[(k+1)\left(\frac{rc}{k+1}-c-k\alpha\right)+3(\alpha+\beta)\right]+x_{3}(2-2x_{1}-3x_{3})(k-2)(k+3)(\alpha+\beta)\right\\}.$
(S77)

Substituting the value of $\mathbf{x}^{(D)}$ into Eq. (S73), we have

$\left.J\right|_{\mathbf{x}=\mathbf{x}^{(D)}}=\frac{\delta(k-2)}{2(k-1)}\begin{pmatrix}\displaystyle{(k+1)\left(\frac{rc}{k+1}-c\right)}&\displaystyle{0}\\\\[10.00002pt]
\displaystyle{0}&\displaystyle{(k+1)\left(\frac{rc}{k+1}-c-k\alpha\right)+3(\alpha+\beta)}\end{pmatrix}.$
(S78)

The condition ensuring $J$ negative-definite is $r<k+1$ and
$r<(1+k\alpha/c)(k+1)-3(\alpha+\beta)/c$. It is hard to compare them in size,
but we can express the second condition as $x_{3}^{(DE)}>0$, which holds if
$\mathbf{x}^{(DE)}$ exists. Therefore, if $r<k+1$ and $\mathbf{x}^{(DE)}$
exists, the equilibrium $\mathbf{x}^{(D)}$ is stable.

Substitute the value of $\mathbf{x}^{(DE)}$ into Eq. (S73), we have
$\left.\partial\dot{x}_{1}/\partial
x_{3}\right|_{\mathbf{x}=\mathbf{x}^{(DE)}}=0$. Therefore, the first- and
second-order sequential principal subformulas’ negativity is equivalent to
$\left.\partial\dot{x}_{1}/\partial
x_{1}\right|_{\mathbf{x}=\mathbf{x}^{(DE)}}<0$ and
$\left.\partial\dot{x}_{3}/\partial
x_{3}\right|_{\mathbf{x}=\mathbf{x}^{(DE)}}<0$. We have

$\left.\frac{\partial\dot{x}_{1}}{\partial
x_{1}}\right|_{\mathbf{x}=\mathbf{x}^{(DE)}}=\frac{\delta(k-2)}{2(k-1)}(1-x_{3}^{(DE)})\left[(k+1)k\alpha-3(\alpha+\beta)\right],$
(S79)

whose sign is hard to judge. Instead, let us calculate

$\left.\frac{\partial\dot{x}_{3}}{\partial
x_{3}}\right|_{\mathbf{x}=\mathbf{x}^{(DE)}}=\frac{\delta(k-2)}{2(k-1)}(1-x_{3}^{(DE)})\left[(k+1)\left(-\frac{rc}{k+1}+c+k\alpha\right)-3(\alpha+\beta)\right].$
(S80)

We can see that as long as $0<x_{3}^{(DE)}<1$, we have
$\left.\partial\dot{x}_{3}/\partial
x_{3}\right|_{\mathbf{x}=\mathbf{x}^{(DE)}}>0$. Therefore, if
$\mathbf{x}^{(DE)}$ exists, $\mathbf{x}^{(DE)}$ cannot be stable.

Next, we study the stability of
$\mathbf{x}^{(CE)}=(x_{1}^{(CE)},0,x_{3}^{(CE)})$, which contains infinite
equilibrium points satisfying $x_{1}^{(CE)}+x_{3}^{(CE)}=1$. Similar to well-
mixed populations, this is because when $x_{2}=0$, we have
$\dot{x}_{1}=\dot{x}_{3}$ and strategy 1 and 3 are indistinguishable.
Concerning this, we can treat $x_{1}+x_{3}$ as a whole and study the system
depicted by $\dot{x}_{2}$ and $\dot{x}_{1}+\dot{x}_{3}$. We can cancel
$x_{1}+x_{3}=1-x_{2}$ and have

$\dot{x}_{2}=\frac{\delta(k-2)(k+1)}{2(k-1)}\left\\{-x_{2}(1-x_{2})\left(\frac{rc}{k+1}-c\right)+kx_{2}x_{3}\left(\frac{3-(k-2)(k+3)x_{2}}{k(k+1)}(\alpha+\beta)-\beta\right)\right\\}.$
(S81)

The element in the single-order Jacobian matrix is

$\frac{\mathrm{d}\dot{x}_{2}}{\mathrm{d}x_{2}}=\frac{\delta(k-2)(k+1)}{2(k-1)}\left\\{-(1-2x_{2})\left(\frac{rc}{k+1}-c\right)+kx_{3}\left(\frac{3-2(k-2)(k+3)x_{2}}{k(k+1)}(\alpha+\beta)-\beta\right)\right\\}.$
(S82)

Substituting the value of $\mathbf{x}^{(CE)}$ into
$\mathrm{d}\dot{x}_{2}/\mathrm{d}x_{2}$, we have

$\left.\frac{\mathrm{d}\dot{x}_{2}}{\mathrm{d}x_{2}}\right|_{\mathbf{x}=\mathbf{x}^{(CE)}}=\frac{\delta(k-2)(k+1)}{2(k-1)}\left\\{-\left(\frac{rc}{k+1}-c\right)-kx_{3}^{(CE)}\left(\frac{3}{k(k+1)}(\alpha+\beta)-\beta\right)\right\\}.$
(S83)

Therefore, $\mathbf{x}^{(CE)}$ is stable if $x_{1}^{(CE)}<x_{1,\star}^{(CE)}$
(or $x_{3}^{(CE)}>x_{3,\star}^{(CE)}$), where

$x_{1,\star}^{(CE)}=1+\left(\frac{r}{k+1}-1\right)\frac{c}{k\beta-3(\alpha+\beta)/(k+1)}\equiv\frac{(k-2)(k+3)(\alpha+\beta)}{k(k+1)\beta-3(\alpha+\beta)}x_{2}^{(DE)},$
(S84)

or
$x_{3,\star}^{(CE)}=\left(-\frac{r}{k+1}+1\right)\frac{c}{k\beta-3(\alpha+\beta)/(k+1)}$.
That is, while every point on the $CE$-edge is equilibrium, only the points on
the side of $x_{1}^{(CE)}<x_{1,\star}^{(CE)}$ are stable. Moreover, we note
that $x_{1,\star}^{(CE)}>1$ when $r>k+1$. That is,
$x_{1}^{(CE)}<x_{1,\star}^{(CE)}$ always holds and the points on the $CE$-edge
are stable everywhere for $r>k+1$. Please note that the conclusions here are
based on a regular size of $\beta$. If the value of $\beta$ is very small,
then the conclusion may be reversed, which is a different phenomenon compared
with the well-mixed population.

####  S3.2.3 Discussion

According to stability analysis, we can divide the results into two cases.
First, when $r>k+1$, the fixation of defection at $\mathbf{x}^{(D)}$ is
unstable and the points on the $CE$-edge are stable everywhere. The dilemma is
overcomed through the known principles as in the traditional public goods
game.

Second, when $r<k+1$, the fixation of defection occurs in the traditional
public goods game, but in this 3-strategy system with peer punishment, the
situations can be different. Here, we discuss the effect of peer punishment
when $r<k+1$ and show how the results in Fig. E3 in the main text are
obtained.

As mentioned in the main text, when the punishment strength is intermediate
($\beta_{0}<\beta<\beta^{\star}$ for structured populations and
$\beta>\beta_{0}^{\text{WM}}$ for well-mixed populations, as revealed later),
there are a stable vertex equilibrium $\mathbf{x}^{(D)}$, an unstable edge
equilibrium $\mathbf{x}^{(DE)}$, and a stable equilibrium line
$\mathbf{x}^{(CE)}$ on the $CE$-edge. The vertex equilibrium
$\mathbf{x}^{(D)}$ and the equilibrium line $\mathbf{x}^{(CE)}$ are bi-stable,
depending on the initial state space divided by $\mathbf{x}^{(DE)}$ and
$\mathbf{x}_{\star}^{(CE)}=(x_{1,\star}^{(CE)},0,x_{3,\star}^{(CE)})$. In
well-mixed populations, $x_{2}^{(DE)}\to 1$ as $\beta\to\infty$ according to
Eq. (S58a), and $\mathbf{x}^{(D)}$ is always stable regardless of the value of
$\beta$, according to Eq. (S61). This is the analytical illustration that peer
punishment cannot truly resolve social dilemmas in a well-mixed population.
However, in structured populations, $x_{2}^{(DE)}\to k(k+1)/[(k-2)(k+3)]>1$ as
$\beta\to\infty$ according to Eq. (S71a). In particular, $x_{2}^{(DE)}>1$
happens (i.e., $\mathbf{x}^{(DE)}$ no longer exists) when

$\displaystyle\leavevmode\nobreak\
\frac{1}{(k-2)(k+3)(\alpha+\beta)}\left[(k+1)\left(\frac{rc}{k+1}-c+k\beta\right)-3(\alpha+\beta)\right]>1$
$\displaystyle\Leftrightarrow$ $\displaystyle\leavevmode\nobreak\
\beta>\beta^{\star}\equiv\frac{k+1}{3}\left(-\frac{rc}{k+1}+c+k\alpha\right)-\alpha.$
(S85)

Meanwhile, according to Eq. (S78), $\mathbf{x}^{(D)}$ becomes unstable at the
same time when Eq. (S3.2.3) is satisfied. The bi-stable system state becomes
mono-stable on the $CE$-edge. This is the analytical interpretation that peer
punishment can resolve social dilemmas in structured populations.

On the other hand, structured populations do not always facilitate the
advantage of peer punishment compared to well-mixed populations. This property
can be observed from the critical punishment strength over which peer
punishment starts to play a role. In a structured population, according to Eq.
(S71a), $x_{2}^{(DE)}>0$ means

$\displaystyle\leavevmode\nobreak\
\frac{1}{(k-2)(k+3)(\alpha+\beta)}\left[(k+1)\left(\frac{rc}{k+1}-c+k\beta\right)-3(\alpha+\beta)\right]>0$
$\displaystyle\Leftrightarrow$ $\displaystyle\leavevmode\nobreak\
\beta>\beta_{0}\equiv\frac{k+1}{k^{2}+k-3}\left(-\frac{rc}{k+1}+c\right)+\frac{3\alpha}{k^{2}+k-3}.$
(S86)

And, according to Eq. (S84), $x_{1,\star}^{(CE)}>0\Leftrightarrow
x_{2}^{(DE)}>0$. When $\beta<\beta_{0}$, the system is mono-stable at the
vertex equilibrium $\mathbf{x}^{(D)}$, as it is in the traditional public
goods game. Peer punishment starts functioning and creates bi-stability in
structured populations when $\beta>\beta_{0}$.

In a well-mixed population, according to Eq. (S58a), $x_{2}^{(DE)}>0$ means

$\frac{1}{k(\alpha+\beta)}\left(\frac{rc}{k+1}-c+k\beta\right)>0\Leftrightarrow\beta>\beta_{0}^{\text{WM}}\equiv\frac{1}{k}\left(-\frac{rc}{k+1}+c\right),$
(S87)

and $x_{1,\star}^{(CE)}>0\Leftrightarrow x_{2}^{(DE)}>0$ according to Eq.
(S67). Peer punishment starts working and creates bi-stability in well-mixed
populations when $\beta>\beta_{0}^{\text{WM}}$.

Comparing Eqs. (S3.2.3) and (S87), we see that
$\beta_{0}^{\text{WM}}<\beta_{0}$ always holds. That is, peer punishment first
starts to play a role in well-mixed populations when increasing the punishment
strength.

Actually, there is a considerable interval of $\beta$ that $x_{2}^{(DE)}$ in
structured populations are smaller than the ones in well-mixed populations.
That is, a structured population enlarges the initial state space leading to
full defection, thus less effective to utilize peer punishment. A structured
population becomes more advantages only when its $x_{2}^{(DE)}$ are greater
than the ones in well-mixed populations. According to Eqs. (S58a) and (S71a),
this means

$\displaystyle\leavevmode\nobreak\
\frac{1}{(k-2)(k+3)(\alpha+\beta)}\left[(k+1)\left(\frac{rc}{k+1}-c+k\beta\right)-3(\alpha+\beta)\right]>\frac{1}{k(\alpha+\beta)}\left(\frac{rc}{k+1}-c+k\beta\right)$
$\displaystyle\Leftrightarrow$ $\displaystyle\leavevmode\nobreak\
\beta>\beta_{=}\equiv\frac{2}{k}\left(-\frac{rc}{k+1}+c\right)+\alpha.$ (S88)

When $\beta>\beta_{=}$, a structured population has a smaller initial state
space leading to full defection, thus more effective to utilize peer
punishment than well-mixed populations. In particular, $\beta>\beta^{\star}$
completely eliminates the initial state space leading to full defection as
shown in Eq. (S3.2.3), while a well-mixed population cannot realize this.

###  S3.3 Public goods games with pool punishment ($n=3$)

As mentioned in the main text, there are $n=3$ strategies in the public goods
game with pool punishment [8].

$1=\text{Cooperation}$ ($C$);

$2=\text{Defection}$ ($D$);

$3=\text{Pool punishment}$ ($O$);

Based on the traditional public goods game, the pool punishment strategy is
introduced as an additional strategy. A punishing player always pays a cost
$\alpha$ to establish the institution for pool punishment. A defective player
is charged with a fine $\beta$ if there is at least one punishing co-player.
That is, given $k_{3}$ punishing co-players, a defective player has $\beta
f(k_{3})$ charged, where $f(k_{3})=1$ if $k_{3}>0$ and $f(k_{3})=0$ if
$k_{3}=0$. Meanwhile, we assume punishing players also perform the cooperative
behavior, investing $c$ to the common pool.

Therefore, given the co-player configuration $\mathbf{k}=(k_{1},k_{2},k_{3})$,
we have the following payoff calculation in a single game.

$\displaystyle a_{1|\mathbf{k}}$
$\displaystyle=\frac{r(k_{1}+1+k_{3})c}{k+1}-c=\frac{rc}{k+1}k_{1}+\frac{rc}{k+1}k_{3}+\frac{rc}{k+1}-c,$
(S89a) $\displaystyle a_{2|\mathbf{k}}$
$\displaystyle=\frac{r(k_{1}+k_{3})c}{k+1}-\beta
f(k_{3})=\frac{rc}{k+1}k_{1}+\frac{rc}{k+1}k_{3}-\beta f(k_{3}),$ (S89b)
$\displaystyle a_{3|\mathbf{k}}$
$\displaystyle=\frac{r(k_{1}+k_{3}+1)c}{k+1}k_{1}-c-\alpha=\frac{rc}{k+1}k_{1}+\frac{rc}{k+1}k_{3}+\frac{rc}{k+1}-c-\alpha.$
(S89c)

####  S3.3.1 The well-mixed population

The frequency of strategies 1, 2, 3 is denoted by $x_{1}$, $x_{2}$, and
$x_{3}$ (or $x_{C}$, $x_{D}$, and $x_{O}$ in the main text), respectively. In
a well-mixed population, the expected payoff of the three strategies are
calculated as follows.

$\displaystyle\bar{\pi}_{1}$
$\displaystyle=\sum_{k_{1}+k_{2}+k_{3}=k}\frac{k!}{k_{1}!k_{2}!k_{3}!}{x_{1}}^{k_{1}}{x_{2}}^{k_{2}}{x_{3}}^{k_{3}}a_{1|\mathbf{k}}=\frac{rc}{k+1}kx_{1}+\frac{rc}{k+1}kx_{3}+\frac{rc}{k+1}-c,$
(S90a) $\displaystyle\bar{\pi}_{2}$
$\displaystyle=\sum_{k_{1}+k_{2}+k_{3}=k}\frac{k!}{k_{1}!k_{2}!k_{3}!}{x_{1}}^{k_{1}}{x_{2}}^{k_{2}}{x_{3}}^{k_{3}}a_{2|\mathbf{k}}=\frac{rc}{k+1}kx_{1}+\frac{rc}{k+1}kx_{3}-\beta\left[1-(1-x_{3})^{k}\right],$
(S90b) $\displaystyle\bar{\pi}_{3}$
$\displaystyle=\sum_{k_{1}+k_{2}+k_{3}=k}\frac{k!}{k_{1}!k_{2}!k_{3}!}{x_{1}}^{k_{1}}{x_{2}}^{k_{2}}{x_{3}}^{k_{3}}a_{3|\mathbf{k}}=\frac{rc}{k+1}kx_{1}+\frac{rc}{k+1}kx_{3}+\frac{rc}{k+1}-c-\alpha.$
(S90c)

The expected payoff of the total population is then calculated by

$\bar{\pi}=x_{1}\bar{\pi}_{1}+x_{2}\bar{\pi}_{2}+x_{3}\bar{\pi}_{3}=\frac{rc}{k+1}kx_{1}+\frac{rc}{k+1}kx_{3}+(x_{1}+x_{3})\left(\frac{rc}{k+1}-c\right)-x_{3}\alpha-
x_{2}\left[1-(1-x_{3})^{k}\right].$ (S91)

On this basis, we can write the replicator equations of the well-mixed
population $\dot{x}_{i}=x_{i}(\bar{\pi}_{i}-\bar{\pi})$ as follows.

$\displaystyle\dot{x}_{1}$
$\displaystyle=x_{1}(\bar{\pi}_{1}-\bar{\pi})=x_{1}\left\\{(1-x_{1}-x_{3})\left(\frac{rc}{k+1}-c\right)+x_{3}\alpha+x_{2}\beta\left[1-(1-x_{3})^{k}\right]\right\\},$
(S92a) $\displaystyle\dot{x}_{2}$
$\displaystyle=x_{2}(\bar{\pi}_{2}-\bar{\pi})=x_{2}\left\\{-(x_{1}+x_{3})\left(\frac{rc}{k+1}-c\right)+x_{3}\alpha-(1-x_{2})\beta\left[1-(1-x_{3})^{k}\right]\right\\},$
(S92b) $\displaystyle\dot{x}_{3}$
$\displaystyle=x_{3}(\bar{\pi}_{3}-\bar{\pi})=x_{3}\left\\{(1-x_{1}-x_{3})\left(\frac{rc}{k+1}-c\right)-(1-x_{3})\alpha+x_{2}\beta\left[1-(1-x_{3})^{k}\right]\right\\}.$
(S92c)

We denote the system state $\mathbf{x}=(x_{1},x_{2},x_{3})$. Solving
$\dot{\mathbf{x}}=\mathbf{0}$, we obtain four possible equilibrium points: a
point on the $C$-vertex, $\mathbf{x}^{(C)}=(1,0,0)$, a point on the
$D$-vertex, $\mathbf{x}^{(D)}=(0,1,0)$, a point on the $O$-vertex,
$\mathbf{x}^{(O)}=(0,0,1)$, and a point on the $DO$-edge,
$\mathbf{x}^{(DO)}=(0,x_{2}^{(DO)},x_{3}^{(DO)})$, where

$\displaystyle x_{2}^{(DO)}$
$\displaystyle=\sqrt[k]{1+\frac{1}{\beta}\left(\frac{rc}{k+1}-c-\alpha\right)},$
(S93a) $\displaystyle x_{3}^{(DO)}$
$\displaystyle=1-x_{2}^{(DO)}=1-\sqrt[k]{1+\frac{1}{\beta}\left(\frac{rc}{k+1}-c-\alpha\right)}.$
(S93b)

The stability of these four equilibrium points can be studied by the regular
method. We cancel $x_{2}=1-x_{1}-x_{3}$ and study the dynamics depicted by
$\dot{x}_{1}$ and $\dot{x}_{3}$,

$\displaystyle\dot{x}_{1}$
$\displaystyle=x_{1}\left\\{(1-x_{1}-x_{3})\left(\frac{rc}{k+1}-c\right)+x_{3}\alpha+(1-x_{1}-x_{3})\left[1-(1-x_{3})^{k}\right]\right\\},$
(S94a) $\displaystyle\dot{x}_{3}$
$\displaystyle=x_{3}\left\\{(1-x_{1}-x_{3})\left(\frac{rc}{k+1}-c\right)-(1-x_{3})\alpha+(1-x_{1}-x_{3})\left[1-(1-x_{3})^{k}\right]\right\\}.$
(S94b)

The Jacobian matrix of system (S94) is

$J=\begin{pmatrix}\displaystyle{\frac{\partial\dot{x}_{1}}{\partial
x_{1}}}&\displaystyle{\frac{\partial\dot{x}_{1}}{\partial
x_{3}}}\\\\[10.00002pt] \displaystyle{\frac{\partial\dot{x}_{3}}{\partial
x_{1}}}&\displaystyle{\frac{\partial\dot{x}_{3}}{\partial
x_{3}}}\end{pmatrix},$ (S95)

where

$\displaystyle\frac{\partial\dot{x}_{1}}{\partial x_{1}}=$
$\displaystyle\leavevmode\nobreak\
(1-2x_{1}-x_{3})\left(\frac{rc}{k+1}-c\right)+x_{3}\alpha+(1-2x_{1}-x_{3})\left[1-(1-x_{3})^{k}\right],$
(S96) $\displaystyle\frac{\partial\dot{x}_{1}}{\partial x_{3}}=$
$\displaystyle\leavevmode\nobreak\
x_{1}\left\\{-\frac{rc}{k+1}+c+\alpha-\beta+(1-x_{3})^{k-1}\beta\left[1-x_{3}+k(1-x_{1}-x_{3})\right]\right\\},$
(S97) $\displaystyle\frac{\partial\dot{x}_{3}}{\partial x_{1}}=$
$\displaystyle-
x_{3}\left\\{\frac{rc}{k+1}-c+\beta\left[1-(1-x_{3})^{k}\right]\right\\},$
(S98) $\displaystyle\frac{\partial\dot{x}_{3}}{\partial x_{3}}=$
$\displaystyle\leavevmode\nobreak\
(1-x_{1}-2x_{3})\left(\frac{rc}{k+1}-c\right)-(1-2x_{3})\alpha$
$\displaystyle+\beta\left\\{1-x_{1}-2x_{3}+(1-x_{3})^{k-1}\left[-(1-x_{1}-x_{3})+(k+2)(1-x_{3})x_{3}-(k+1)x_{1}x_{3}\right]\right\\}.$
(S99)

Substituting the value of $\mathbf{x}^{(C)}$ into Eq. (S95), we have

$\left.J\right|_{\mathbf{x}=\mathbf{x}^{(C)}}=\begin{pmatrix}\displaystyle{-\frac{rc}{k+1}+c}&\displaystyle{-\frac{rc}{k+1}+c+\alpha-\beta}\\\\[10.00002pt]
\displaystyle{0}&\displaystyle{-\alpha}\end{pmatrix}.$ (S100)

The condition ensuring $J$ negative-definite is $r>k+1$. Therefore, the
$C$-vertex equilibrium $\mathbf{x}^{(C)}$ is stable if and only if $r>k+1$.

Substituting the value of $\mathbf{x}^{(D)}$ into Eq. (S95), we have

$\left.J\right|_{\mathbf{x}=\mathbf{x}^{(D)}}=\begin{pmatrix}\displaystyle{\frac{rc}{k+1}-c}&\displaystyle{0}\\\\[10.00002pt]
\displaystyle{0}&\displaystyle{\frac{rc}{k+1}-c-\alpha}\end{pmatrix}.$ (S101)

The condition ensuring $J$ negative-definite is $r<k+1$. Therefore, the
$D$-vertex equilibrium $\mathbf{x}^{(D)}$ is stable if and only if $r<k+1$.

Substituting the value of $\mathbf{x}^{(O)}$ into Eq. (S95), we have
$\left.\partial\dot{x}_{1}/\partial
x_{1}\right|_{\mathbf{x}=\mathbf{x}^{(O)}}=\alpha>0$. Therefore, $J$ is not
negative-definite and the $O$-vertex equilibrium $\mathbf{x}^{(O)}$ is not
stable.

Substituting the value of $\mathbf{x}^{(DO)}$ into Eq. (S95), we note that

$\left.\frac{\partial\dot{x}_{1}}{\partial
x_{1}}\right|_{\mathbf{x}=\mathbf{x}^{(DO)}}=\sqrt[k]{1+\frac{1}{\beta}\left(\frac{rc}{k+1}-c-\alpha\right)}\left\\{\left(\frac{rc}{k+1}-c-\alpha\right)+\beta\left[-\frac{1}{\beta}\left(\frac{rc}{k+1}-c-\alpha\right)\right]\right\\}+\alpha=\alpha>0.$
(S102)

Therefore, $J$ is not negative-definite and the equilibrium
$\mathbf{x}^{(DO)}$ on the $DO$-edge is not stable.

####  S3.3.2 The structured population

We notice that the payoff structure given by Eq. (S89) is nonlinear.
Therefore, to obtain the replicator dynamics in structured populations, we
need to utilize Eq. (S26) given by section S2.4.2. The process is to calculate
all elements of the ‘$\langle a_{i|\mathbf{k}_{+j}}\rangle_{i}$ type’ and the
‘$\langle a_{i|\mathbf{k}_{+j}}\rangle_{j}$ type’ manually, which is tedious.

For this $n=3$ system, the $\langle a_{i|\mathbf{k}_{+j}}\rangle_{i}$ type is

$\left[\langle
a_{i|\mathbf{k}_{+j}}\rangle_{i}\right]_{ij}=\displaystyle{\begin{pmatrix}\langle
a_{1|\mathbf{k}_{+1}}\rangle_{1}&\langle
a_{1|\mathbf{k}_{+2}}\rangle_{1}&\langle a_{1|\mathbf{k}_{+3}}\rangle_{1}\\\
\langle a_{2|\mathbf{k}_{+1}}\rangle_{2}&\langle
a_{2|\mathbf{k}_{+2}}\rangle_{2}&\langle a_{2|\mathbf{k}_{+3}}\rangle_{2}\\\
\langle a_{3|\mathbf{k}_{+1}}\rangle_{3}&\langle
a_{3|\mathbf{k}_{+2}}\rangle_{3}&\langle
a_{3|\mathbf{k}_{+3}}\rangle_{3}\end{pmatrix}},$ (S103)

where, by applying Eq. (S89) to Eq. (S21), we have

$\displaystyle\langle a_{1|\mathbf{k}_{+1}}\rangle_{1}$
$\displaystyle=\sum_{k_{1}+k_{2}+k_{3}=k-1}\frac{(k-1)!}{k_{1}!k_{2}!k_{3}!}q_{1|1}^{k_{1}}q_{2|1}^{k_{2}}q_{3|1}^{k_{3}}a_{1|\mathbf{k}_{+1}}$
$\displaystyle=\sum_{k_{1}+k_{2}+k_{3}=k-1}\frac{(k-1)!}{k_{1}!k_{2}!k_{3}!}q_{1|1}^{k_{1}}q_{2|1}^{k_{2}}q_{3|1}^{k_{3}}\left(\frac{rc}{k+1}(k_{1}+1)+\frac{rc}{k+1}k_{3}+\frac{rc}{k+1}-c\right)$
$\displaystyle=\frac{rc}{k+1}(k-1)q_{1|1}+\frac{rc}{k+1}(k-1)q_{3|1}+\frac{2rc}{k+1}-c$
$\displaystyle=\frac{rc}{k+1}(k-2)x_{1}+\frac{rc}{k+1}(k-2)x_{3}+\frac{3rc}{k+1}-c,$
(S104)

and similarly,

$\displaystyle\langle a_{1|\mathbf{k}_{+2}}\rangle_{1}$
$\displaystyle=\frac{rc}{k+1}(k-2)x_{1}+\frac{rc}{k+1}(k-2)x_{3}+\frac{2rc}{k+1}-c,$
(S105a) $\displaystyle\langle a_{1|\mathbf{k}_{+3}}\rangle_{1}$
$\displaystyle=\frac{rc}{k+1}(k-2)x_{1}+\frac{rc}{k+1}(k-2)x_{3}+\frac{3rc}{k+1}-c;$
(S105b) $\displaystyle\langle a_{2|\mathbf{k}_{+1}}\rangle_{2}$
$\displaystyle=\frac{rc}{k+1}(k-2)x_{1}+\frac{rc}{k+1}(k-2)x_{3}+\frac{rc}{k+1}-\beta\left[1-\left(1-\frac{k-2}{k-1}x_{3}\right)^{k-1}\right],$
(S105c) $\displaystyle\langle a_{2|\mathbf{k}_{+2}}\rangle_{2}$
$\displaystyle=\frac{rc}{k+1}(k-2)x_{1}+\frac{rc}{k+1}(k-2)x_{3}-\beta\left[1-\left(1-\frac{k-2}{k-1}x_{3}\right)^{k-1}\right],$
(S105d) $\displaystyle\langle a_{2|\mathbf{k}_{+3}}\rangle_{2}$
$\displaystyle=\frac{rc}{k+1}(k-2)x_{1}+\frac{rc}{k+1}(k-2)x_{3}+\frac{rc}{k+1}-\beta;$
(S105e) $\displaystyle\langle a_{3|\mathbf{k}_{+1}}\rangle_{3}$
$\displaystyle=\frac{rc}{k+1}(k-2)x_{1}+\frac{rc}{k+1}(k-2)x_{3}+\frac{3rc}{k+1}-c-\alpha,$
(S105f) $\displaystyle\langle a_{3|\mathbf{k}_{+2}}\rangle_{3}$
$\displaystyle=\frac{rc}{k+1}(k-2)x_{1}+\frac{rc}{k+1}(k-2)x_{3}+\frac{2rc}{k+1}-c-\alpha,$
(S105g) $\displaystyle\langle a_{3|\mathbf{k}_{+3}}\rangle_{3}$
$\displaystyle=\frac{rc}{k+1}(k-2)x_{1}+\frac{rc}{k+1}(k-2)x_{3}+\frac{3rc}{k+1}-c-\alpha.$
(S105h)

Next, the $\langle a_{i|\mathbf{k}_{+j}}\rangle_{j}$ type is

$\left[\langle
a_{i|\mathbf{k}_{+j}}\rangle_{j}\right]_{ij}=\displaystyle{\begin{pmatrix}\langle
a_{1|\mathbf{k}_{+1}}\rangle_{1}&\langle
a_{1|\mathbf{k}_{+2}}\rangle_{2}&\langle a_{1|\mathbf{k}_{+3}}\rangle_{3}\\\
\langle a_{2|\mathbf{k}_{+1}}\rangle_{1}&\langle
a_{2|\mathbf{k}_{+2}}\rangle_{2}&\langle a_{2|\mathbf{k}_{+3}}\rangle_{3}\\\
\langle a_{3|\mathbf{k}_{+1}}\rangle_{1}&\langle
a_{3|\mathbf{k}_{+2}}\rangle_{2}&\langle
a_{3|\mathbf{k}_{+3}}\rangle_{3}\end{pmatrix}},$ (S106)

where the diagonal elements have been calculated previously, and the remaining
elements are

$\displaystyle\langle a_{1|\mathbf{k}_{+2}}\rangle_{2}$
$\displaystyle=\frac{rc}{k+1}(k-2)x_{1}+\frac{rc}{k+1}(k-2)x_{3}+\frac{2rc}{k+1}-c,$
(S107a) $\displaystyle\langle a_{1|\mathbf{k}_{+3}}\rangle_{3}$
$\displaystyle=\frac{rc}{k+1}(k-2)x_{1}+\frac{rc}{k+1}(k-2)x_{3}+\frac{3rc}{k+1}-c;$
(S107b) $\displaystyle\langle a_{2|\mathbf{k}_{+1}}\rangle_{1}$
$\displaystyle=\frac{rc}{k+1}(k-2)x_{1}+\frac{rc}{k+1}(k-2)x_{3}+\frac{2rc}{k+1}-\beta\left[1-\left(1-\frac{k-2}{k-1}x_{3}\right)^{k-1}\right],$
(S107c) $\displaystyle\langle a_{2|\mathbf{k}_{+3}}\rangle_{3}$
$\displaystyle=\frac{rc}{k+1}(k-2)x_{1}+\frac{rc}{k+1}(k-2)x_{3}+\frac{2rc}{k+1}-\beta;$
(S107d) $\displaystyle\langle a_{3|\mathbf{k}_{+1}}\rangle_{1}$
$\displaystyle=\frac{rc}{k+1}(k-2)x_{1}+\frac{rc}{k+1}(k-2)x_{3}+\frac{3rc}{k+1}-c-\alpha,$
(S107e) $\displaystyle\langle a_{3|\mathbf{k}_{+2}}\rangle_{2}$
$\displaystyle=\frac{rc}{k+1}(k-2)x_{1}+\frac{rc}{k+1}(k-2)x_{3}+\frac{rc}{k+1}-c-\alpha.$
(S107f)

According to Eqs. (S3.3.2), (S105a)–(S105h), and (S107a)–(S107f), there are
totaling $(2n-1)n=15$ elements in the two ‘$\langle
a_{i|\mathbf{k}_{+j}}\rangle_{i}$’ and ‘$\langle
a_{i|\mathbf{k}_{+j}}\rangle_{j}$’ types.

We can apply these elements to Eq. (S26), which is

$\dot{x}_{i}=\frac{\delta(k-2)}{2(k-1)}x_{i}\sum_{j=1}^{3}x_{j}\Bigg{(}\langle
a_{i|\mathbf{k}_{+j}}\rangle_{i}+(k-1)\langle
a_{i|\mathbf{k}_{+j}}\rangle_{j}+\langle
a_{i|\mathbf{k}_{+i}}\rangle_{i}-\langle
a_{j|\mathbf{k}_{+i}}\rangle_{j}-\langle
a_{j|\mathbf{k}_{+i}}\rangle_{i}-(k-2)\sum_{j=1}^{n}x_{l}\langle
a_{j|\mathbf{k}_{+l}}\rangle_{l}-\langle
a_{j|\mathbf{k}_{+j}}\rangle_{j}\Bigg{)}.$ (S108)

Let us do this step by step. For $i=1$, we have

$\displaystyle\sum_{j=1}^{3}x_{j}\langle a_{1|\mathbf{k}_{+j}}\rangle_{1}=$
$\displaystyle\leavevmode\nobreak\
\frac{rc}{k+1}(k-1)x_{1}+\frac{rc}{k+1}(k-1)x_{3}+\frac{2rc}{k+1}-c,$ (S109a)
$\displaystyle\sum_{j=1}^{3}x_{j}\langle a_{1|\mathbf{k}_{+j}}\rangle_{j}=$
$\displaystyle\leavevmode\nobreak\
\frac{rc}{k+1}kx_{1}+\frac{rc}{k+1}kx_{3}+\frac{rc}{k+1}-c,$ (S109b)
$\displaystyle\sum_{j=1}^{3}x_{j}\langle a_{1|\mathbf{k}_{+1}}\rangle_{1}=$
$\displaystyle\leavevmode\nobreak\
\frac{rc}{k+1}(k-2)x_{1}+\frac{rc}{k+1}(k-2)x_{3}+\frac{3rc}{k+1}-c,$ (S109c)
$\displaystyle\sum_{j=1}^{3}x_{j}\langle a_{j|\mathbf{k}_{+1}}\rangle_{1}=$
$\displaystyle\leavevmode\nobreak\
\frac{rc}{k+1}(k-1)x_{1}+\frac{rc}{k+1}(k-1)x_{3}+\frac{2rc}{k+1}-c(x_{1}+x_{3})-\alpha
x_{3}$
$\displaystyle-\beta(1-x_{1}-x_{3})\left[1-\left(1-\frac{k-2}{k-1}x_{3}\right)^{k-1}\right],$
(S109d) $\displaystyle\sum_{j=1}^{3}x_{j}\langle
a_{j|\mathbf{k}_{+1}}\rangle_{j}=$ $\displaystyle\leavevmode\nobreak\
\frac{rc}{k+1}kx_{1}+\frac{rc}{k+1}kx_{3}+\frac{rc}{k+1}-c(x_{1}+x_{3})-\alpha
x_{3}-\beta(1-x_{1}-x_{3})\left[1-\left(1-\frac{k-2}{k-1}x_{3}\right)^{k-1}\right],$
(S109e) $\displaystyle\sum_{j=1}^{3}x_{j}\sum_{l=1}^{3}x_{l}\langle
a_{j|\mathbf{k}_{+l}}\rangle_{l}=$ $\displaystyle\leavevmode\nobreak\
(rc-c)(x_{1}+x_{3})-\alpha
x_{3}-\beta(1-x_{1}-x_{3})(1-x_{3})\left[1-\left(1-\frac{k-2}{k-1}x_{3}\right)^{k-1}\right]-\beta(1-x_{1}-x_{3})x_{3},$
(S109f) $\displaystyle\sum_{j=1}^{3}x_{j}\langle
a_{j|\mathbf{k}_{+j}}\rangle_{j}=$ $\displaystyle\leavevmode\nobreak\
(rc-c)(x_{1}+x_{3})-\alpha
x_{3}-\beta(1-x_{1}-x_{3})\left[1-\left(1-\frac{k-2}{k-1}x_{3}\right)^{k-1}\right],$
(S109g)

which can be used to obtain $\dot{x}_{1}$.

For $i=3$, $\sum_{j=1}^{3}x_{j}\sum_{l=1}^{3}x_{l}\langle
a_{j|\mathbf{k}_{+l}}\rangle_{l}$ and $\sum_{j=1}^{3}x_{j}\langle
a_{j|\mathbf{k}_{+j}}\rangle_{j}$ have been obtained previously, and we have
the remaining terms:

$\displaystyle\sum_{j=1}^{3}x_{j}\langle a_{3|\mathbf{k}_{+j}}\rangle_{3}=$
$\displaystyle\leavevmode\nobreak\
\frac{rc}{k+1}(k-1)x_{1}+\frac{rc}{k+1}(k-1)x_{3}+\frac{2rc}{k+1}-c-\alpha,$
(S110a) $\displaystyle\sum_{j=1}^{3}x_{j}\langle
a_{3|\mathbf{k}_{+j}}\rangle_{j}=$ $\displaystyle\leavevmode\nobreak\
\frac{rc}{k+1}kx_{1}+\frac{rc}{k+1}kx_{3}+\frac{rc}{k+1}-c-\alpha,$ (S110b)
$\displaystyle\sum_{j=1}^{3}x_{j}\langle a_{3|\mathbf{k}_{+3}}\rangle_{3}=$
$\displaystyle\leavevmode\nobreak\
\frac{rc}{k+1}(k-2)x_{1}+\frac{rc}{k+1}(k-2)x_{3}+\frac{3rc}{k+1}-c-\alpha,$
(S110c) $\displaystyle\sum_{j=1}^{3}x_{j}\langle
a_{j|\mathbf{k}_{+3}}\rangle_{3}=$ $\displaystyle\leavevmode\nobreak\
\frac{rc}{k+1}(k-1)x_{1}+\frac{rc}{k+1}(k-1)x_{3}+\frac{2rc}{k+1}-c(x_{1}+x_{3})-\alpha
x_{3}-\beta(1-x_{1}-x_{3}),$ (S110d) $\displaystyle\sum_{j=1}^{3}x_{j}\langle
a_{j|\mathbf{k}_{+3}}\rangle_{j}=$ $\displaystyle\leavevmode\nobreak\
\frac{rc}{k+1}kx_{1}+\frac{rc}{k+1}kx_{3}+\frac{rc}{k+1}-c(x_{1}+x_{3})-\alpha
x_{3}-\beta(1-x_{1}-x_{3}),$ (S110e)

which can be used to obtain $\dot{x}_{3}$.

Applying Eqs. (S109a)–(S109g) and (S110a)–(S110e) to Eq. (S108), we obtain the
following simple replicator equation:

$\displaystyle\dot{x}_{1}$
$\displaystyle=\frac{\delta(k-2)(k+1)}{2(k-1)}x_{1}\left\\{(1-x_{1}-x_{3})\left(\frac{rc}{k+1}-c\right)+x_{3}\alpha+(1-x_{1}-x_{3})\beta\left[1-\left(1-\frac{k-2}{k+1}x_{3}\right)\left(1-\frac{k-2}{k-1}x_{3}\right)^{k-1}\right]\right\\},$
(S111a) $\displaystyle\dot{x}_{3}$
$\displaystyle=\frac{\delta(k-2)(k+1)}{2(k-1)}x_{3}\left\\{(1-x_{1}-x_{3})\left(\frac{rc}{k+1}-c\right)-(1-x_{3})\alpha+(1-x_{1}-x_{3})\beta\left[1-\frac{k-1}{k+1}\left(1-\frac{k-2}{k-1}x_{3}\right)^{k}\right]\right\\}.$
(S111b)

Please note that $\dot{x}_{2}=-\dot{x}_{1}-\dot{x}_{3}$ has been canceled.

Similarly, solving $\dot{\mathbf{x}}=\mathbf{0}$, we obtain four possible
equilibrium points: a point on the $C$-vertex, $\mathbf{x}^{(C)}=(1,0,0)$, a
point on the $D$-vertex, $\mathbf{x}^{(D)}=(0,1,0)$, a point on the
$O$-vertex, $\mathbf{x}^{(O)}=(0,0,1)$, and a point on the $DO$-edge,
$\mathbf{x}^{(DO)}=(0,x_{2}^{(DO)},x_{3}^{(DO)})$, where

$\displaystyle x_{2}^{(DO)}$
$\displaystyle=\frac{k-1}{k-2}\left(-\frac{1}{k-1}+\sqrt[k]{\frac{k+1}{k-1}\left[1+\frac{1}{\beta}\left(\frac{rc}{k+1}-c-\alpha\right)\right]}\right),$
(S112a) $\displaystyle x_{3}^{(DO)}$
$\displaystyle=1-x_{2}^{(DO)}=\frac{k-1}{k-2}\left(1-\sqrt[k]{\frac{k+1}{k-1}\left[1+\frac{1}{\beta}\left(\frac{rc}{k+1}-c-\alpha\right)\right]}\right).$
(S112b)

The Jacobian matrix of system (S111) is

$J=\begin{pmatrix}\displaystyle{\frac{\partial\dot{x}_{1}}{\partial
x_{1}}}&\displaystyle{\frac{\partial\dot{x}_{1}}{\partial
x_{3}}}\\\\[10.00002pt] \displaystyle{\frac{\partial\dot{x}_{3}}{\partial
x_{1}}}&\displaystyle{\frac{\partial\dot{x}_{3}}{\partial
x_{3}}}\end{pmatrix},$ (S113)

where

$\displaystyle\frac{\partial\dot{x}_{1}}{\partial x_{1}}=$
$\displaystyle\leavevmode\nobreak\
\frac{\delta(k-2)(k+1)}{2(k-1)}\Bigg{\\{}(1-2x_{1}-x_{3})\left(\frac{rc}{k+1}-c\right)+x_{3}\alpha$
$\displaystyle+(1-2x_{1}-x_{3})\beta\left[1-\left(1-\frac{k-2}{k+1}x_{3}\right)\left(1-\frac{k-2}{k-1}x_{3}\right)^{k-1}\right]\Bigg{\\}},$
(S114) $\displaystyle\frac{\partial\dot{x}_{1}}{\partial x_{3}}=$
$\displaystyle\leavevmode\nobreak\
\frac{\delta(k-2)(k+1)}{2(k-1)}x_{1}\Bigg{\\{}-\frac{rc}{k+1}+c+\alpha-\beta\left[1-\left(1-\frac{k-2}{k+1}x_{3}\right)\left(1-\frac{k-2}{k-1}x_{3}\right)^{k-1}\right]$
$\displaystyle+(1-x_{1}-x_{3})\beta\left(1-\frac{k-2}{k-1}x_{3}\right)^{k-2}\frac{k-2}{k+1}\left(k+2-\frac{k(k-2)}{k-1}x_{3}\right)\Bigg{\\}},$
(S115) $\displaystyle\frac{\partial\dot{x}_{3}}{\partial x_{1}}=$
$\displaystyle-\frac{\delta(k-2)(k+1)}{2(k-1)}x_{3}\left\\{\frac{rc}{k+1}-c+\beta\left[1-\frac{k-1}{k+1}\left(1-\frac{k-2}{k-1}x_{3}\right)^{k}\right]\right\\},$
(S116) $\displaystyle\frac{\partial\dot{x}_{3}}{\partial x_{3}}=$
$\displaystyle\leavevmode\nobreak\
\frac{\delta(k-2)(k+1)}{2(k-1)}\Bigg{\\{}(1-x_{1}-2x_{3})\left(\frac{rc}{k+1}-c\right)-(1-2x_{3})\alpha+(1-x_{1}-2x_{3})\beta\left[1-\frac{k-1}{k+1}\left(1-\frac{k-2}{k-1}x_{3}\right)^{k}\right]$
$\displaystyle+x_{3}(1-x_{1}-x_{3})\beta\frac{k(k-2)}{k+1}\left(1-\frac{k-2}{k-1}x_{3}\right)^{k-1}\Bigg{\\}}.$
(S117)

Substituting the value of $\mathbf{x}^{(C)}$ into Eq. (S113), we have

$\left.J\right|_{\mathbf{x}=\mathbf{x}^{(C)}}=\frac{\delta(k-2)(k+1)}{2(k-1)}\begin{pmatrix}\displaystyle{-\frac{rc}{k+1}+c}&\displaystyle{-\frac{rc}{k+1}+c+\alpha}\\\\[10.00002pt]
\displaystyle{0}&\displaystyle{-\frac{rc}{k+1}+c-\alpha}\end{pmatrix}.$ (S118)

The condition ensuring $J$ negative-definite is $r>k+1$. Therefore, the
$C$-vertex equilibrium $\mathbf{x}^{(D)}$ is stable if and only if $r>k+1$.

Substituting the value of $\mathbf{x}^{(D)}$ into Eq. (S113), we have

$\left.J\right|_{\mathbf{x}=\mathbf{x}^{(D)}}=\frac{\delta(k-2)(k+1)}{2(k-1)}\begin{pmatrix}\displaystyle{\frac{rc}{k+1}-c}&\displaystyle{0}\\\\[10.00002pt]
\displaystyle{0}&\displaystyle{\frac{rc}{k+1}-c-\alpha+\frac{2\beta}{k+1}}\end{pmatrix}.$
(S119)

The condition ensuring $J$ negative-definite and the $D$-vertex equilibrium
$\mathbf{x}^{(D)}$ stable is $r<k+1$ and $\beta<\beta^{\star}$, where

$\beta^{\star}=\frac{k+1}{2}\left(-\frac{rc}{k+1}+c+\alpha\right).$ (S120)

This indicates that when $\beta>\beta^{\star}$, the $D$-vertex equilibrium
cannot be stable even if $r<k+1$.

Substituting the value of $\mathbf{x}^{(O)}$ into Eq. (S113), we have
$\left.\partial\dot{x}_{1}/\partial
x_{1}\right|_{\mathbf{x}=\mathbf{x}^{(O)}}\sim\alpha>0$. Therefore, $J$ is not
negative-definite and the $O$-vertex equilibrium $\mathbf{x}^{(O)}$ is not
stable.

Substituting the value of $\mathbf{x}^{(DO)}$ into Eq. (S113), we note that
$\left.\partial\dot{x}_{3}/\partial
x_{1}\right|_{\mathbf{x}=\mathbf{x}^{(DO)}}=0$. Therefore, the first- and
second-order sequential principal subformulas’ negativity is equivalent to
$\left.\partial\dot{x}_{1}/\partial
x_{1}\right|_{\mathbf{x}=\mathbf{x}^{(DO)}}<0$ and
$\left.\partial\dot{x}_{3}/\partial
x_{3}\right|_{\mathbf{x}=\mathbf{x}^{(DO)}}<0$. Let us calculate

$\displaystyle\left.\frac{\partial\dot{x}_{3}}{\partial
x_{3}}\right|_{\mathbf{x}=\mathbf{x}^{(DO)}}=$
$\displaystyle\leavevmode\nobreak\
\frac{\delta(k-2)(k+1)}{2(k-1)}\Bigg{\\{}(1-2x_{3}^{(DO)})\left(\frac{rc}{k+1}-c-\alpha+\beta\left[1-\frac{k-1}{k+1}\left(1-\frac{k-2}{k-1}x_{3}^{(DO)}\right)^{k}\right]\right)$
$\displaystyle+x_{3}^{(DO)}(1-x_{3}^{(DO)})\beta\frac{k(k-2)}{k+1}\left(1-\frac{k-2}{k-1}x_{3}^{(DO)}\right)^{k-1}\Bigg{\\}}$
$\displaystyle=$ $\displaystyle\leavevmode\nobreak\
\frac{\delta(k-2)(k+1)}{2(k-1)}x_{3}^{(DO)}(1-x_{3}^{(DO)})\beta\frac{k(k-2)}{k+1}\left(1-\frac{k-2}{k-1}x_{3}^{(DO)}\right)^{k-1}>0.$
(S121)

Therefore, the equilibrium $\mathbf{x}^{(DO)}$ on the $DO$-edge is not stable.

####  S3.3.3 Discussion

According to stability analysis, we can divide the results into two cases.
First, when $r>k+1$, the fixation of defection at $\mathbf{x}^{(D)}$ is
unstable and the fixation of cooperation at $\mathbf{x}^{(C)}$ is stable. The
dilemma is overcomed through the known principles as in the traditional public
goods game.

Second, when $r<k+1$, the fixation of defection occurs in the traditional
public goods game, but in this 3-strategy system with pool punishment, the
situations can be different. Here, we discuss the effect of pool punishment
when $r<k+1$ and show how the results in Fig. E4 in the main text are
obtained.

As mentioned in the main text, when the punishment strength is intermediate
($\beta_{0}<\beta<\beta^{\star}$ for structured populations and
$\beta>\beta_{0}^{\text{WM}}$ for well-mixed populations, as revealed later),
there are a stable vertex equilibrium $\mathbf{x}^{(D)}$, an unstable vertex
equilibrium $\mathbf{x}^{(O)}$, an unstable vertex equilibrium
$\mathbf{x}^{(C)}$, and an unstable edge equilibrium $\mathbf{x}^{(DO)}$. In
the 3-strategy system state space, $\mathbf{x}^{(D)}$ is the only equilibrium
point. However, on the $DO$ edge where $x_{1}=0$, the vertex equilibrium
points $\mathbf{x}^{(D)}$ and $\mathbf{x}^{(O)}$ are bi-stable, depending on
the initial state divided by $\mathbf{x}^{(DO)}$. In well-mixed populations,
$x_{2}^{(DO)}\to 1$ as $\beta\to\infty$ according to Eq. (S93a), and
$\mathbf{x}^{(D)}$ is always stable regardless of the value of $\beta$
according to Eq. (S101). This is the analytical illustration that pool
punishment cannot truly resolve social dilemmas in a well-mixed population.
However, in structured populations,
$x_{2}^{(DO)}\to[(k+1)^{1/k}(k-1)^{1-1/k}-1]/(k-2)>1$ as $\beta\to\infty$
according to Eq. (S112a). In particular, $x_{2}^{(DO)}>1$ happens (i.e.,
$\mathbf{x}^{(DO)}$ no longer exists) when

$\frac{k-1}{k-2}\left(-\frac{1}{k-1}+\sqrt[k]{\frac{k+1}{k-1}\left[1+\frac{1}{\beta}\left(\frac{rc}{k+1}-c-\alpha\right)\right]}\right)>1\Leftrightarrow\beta>\beta^{\star}\equiv\frac{k+1}{2}\left(-\frac{rc}{k+1}+c+\alpha\right)$
(S122)

Meanwhile, according to Eq. (S119), $\mathbf{x}^{(D)}$ becomes unstable at the
same time when Eq. (S122) is satisfied, while $\mathbf{x}^{(O)}$ and
$\mathbf{x}^{(C)}$ remain unstable. In this case, the system enters a rock-
paper-scissors cycle, where defection conquers cooperation
($\partial\dot{x}_{1}/\partial x_{1}<0$ when $x_{3}=0$), cooperation conquers
pool punishment ($\partial\dot{x}_{1}/\partial x_{1}>0$ when $x_{2}=0$), and
pool punishment conquers defection ($\partial\dot{x}_{3}/\partial x_{3}>0$
when $x_{1}=0$). This is the analytical interpretation that pool punishment
can resolve social dilemmas in structured populations.

On the other hand, we also find that structured populations do not always
facilitate the advantage of pool punishment compared to well-mixed
populations. This property can be observed from the critical punishment
strength over which pool punishment starts to play a role. In a structured
population, according to Eq. (S112a), $x_{2}^{(DO)}>0$ means

$\displaystyle\leavevmode\nobreak\
\frac{k-1}{k-2}\left(-\frac{1}{k-1}+\sqrt[k]{\frac{k+1}{k-1}\left[1+\frac{1}{\beta}\left(\frac{rc}{k+1}-c-\alpha\right)\right]}\right)>0$
$\displaystyle\Leftrightarrow$ $\displaystyle\leavevmode\nobreak\
\beta>\beta_{0}\equiv\frac{(k+1)(k-1)^{k-1}}{1-(k+1)(k-1)^{k-1}}\left(\frac{rc}{k+1}-c-\alpha\right).$
(S123)

When $\beta<\beta_{0}$, $\mathbf{x}^{(DO)}$ is not a real vector and the
system is stable at the vertex equilibrium $\mathbf{x}^{(D)}$, even along the
$DO$-edge. Pool punishment starts working and creates bi-stability along the
$DO$-edge in structured populations when $\beta>\beta_{0}$.

In a well-mixed population, according to Eq. (S93a), $x_{2}^{(DO)}>0$ means

$\sqrt[k]{1+\frac{1}{\beta}\left(\frac{rc}{k+1}-c-\alpha\right)}>0\Leftrightarrow\beta>\beta_{0}^{\text{WM}}\equiv-\frac{rc}{k+1}+c+\alpha.$
(S124)

Peer punishment starts functioning and creates bi-stability along the
$DO$-edge in well-mixed populations when $\beta>\beta_{0}^{\text{WM}}$.

Comparing Eqs. (S3.3.3) and (S124), we see that
$\beta_{0}^{\text{WM}}<\beta_{0}$ always holds. That is, peer punishment first
becomes efficient in well-mixed populations when increasing the punishment
strength.

Actually, there is a considerable interval of $\beta$ that $x_{2}^{(DO)}$ in
structured populations are smaller than the ones in well-mixed populations.
That is, a structured population enlarges the initial state interval leading
to full defection along the $DO$-edge, thus less effective to utilize pool
punishment. A structured population becomes more advantages only when its
$x_{2}^{(DO)}$ are greater than the ones in well-mixed populations. According
to Eqs. (S93a) and (S112a), this means

$\displaystyle\leavevmode\nobreak\
\frac{k-1}{k-2}\left(-\frac{1}{k-1}+\sqrt[k]{\frac{k+1}{k-1}\left[1+\frac{1}{\beta}\left(\frac{rc}{k+1}-c-\alpha\right)\right]}\right)>\sqrt[k]{1+\frac{1}{\beta}\left(\frac{rc}{k+1}-c-\alpha\right)}$
$\displaystyle\Leftrightarrow$ $\displaystyle\leavevmode\nobreak\
\beta>\beta_{=}\equiv\frac{\left[(k+1)^{\frac{1}{k}}(k-1)^{1-\frac{1}{k}}-k+2\right]^{k}}{1-\left[(k+1)^{\frac{1}{k}}(k-1)^{1-\frac{1}{k}}-k+2\right]^{k}}\left(\frac{rc}{k+1}-c-\alpha\right).$
(S125)

When $\beta>\beta_{=}$, a structured population has a smaller initial state
interval leading to full defection along the $DO$-edge, thus more effective to
utilize pool punishment than well-mixed populations. In particular,
$\beta>\beta^{\star}$ can even destabilize the full defection equilibrium as
shown in Eq. (S122) and transform the system to a rock-paper-scissors-like
cyclic dominance, while a well-mixed population cannot realize this.

##  Appendix A Death-birth

Here, we supplement the corresponding results under the death-birth rule.

In a unit time, a random focal individual $A$ is selected to die, and $A$’s
neighbors compete for the vacant position proportional to their fitness. The
neighbor $B$ reproduces its strategy to the position of $A$ with the
probability as follows:

$W=\frac{F_{B}}{\sum_{k}{F_{X}}},$ (S126)

where $\sum_{k}{F_{X}}$ represents the total fitness of $A$’s neighbors. If
individual $B$ has a higher fitness among $A$’s neighbors, then $B$ has a
higher probability to reproduce its strategy to the position of $A$. It should
be noted that the focal individual’s fitness is completely ignored under
death-birth.

Below, we analyze the dynamics of the death-birth rule rigorously.

###  A.1 The increase of $i$-players

The increase of $i$-players happens when a focal $j$-player ($j\neq i$) is
selected to update its strategy and an $i$-player takes the position. Given
the focal $j$-player’s neighbor configuration $\mathbf{k}$, the probability
that an $i$-player takes the $j$-player’s position is

$\mathcal{P}(j\leftarrow
i)=\frac{k_{i}F_{i|j}^{\mathbf{k}}}{\sum_{i^{\prime}=1}^{n}k_{i^{\prime}}F_{i^{\prime}|j}^{\mathbf{k}}}=\frac{k_{i}}{k}+\frac{k_{i}}{k}\left(\pi_{i|j}^{\mathbf{k}}-\sum_{i^{\prime}=1}^{n}\frac{k_{i^{\prime}}}{k}\pi_{i^{\prime}|j}^{\mathbf{k}}\right)\delta+\mathcal{O}(\delta^{2}).$
(S127)

Then, we apply it to all possibilities for $j\neq i$ and neighbor
configurations $\mathbf{k}$, obtaining the probability that the number of
$i$-players increases by 1 during a unit time step,

$\displaystyle\mathcal{P}\left(\Delta x_{i}=\frac{1}{N}\right)$
$\displaystyle=\sum_{j=1,j\neq
i}^{n}x_{j}\sum_{\sum_{i^{\prime}=1}^{n}k_{i^{\prime}}=k}\frac{k!}{\prod_{i^{\prime}=1}^{n}k_{i^{\prime}}!}\left(\prod_{i^{\prime}=1}^{n}{q_{i^{\prime}|j}}^{k_{i^{\prime}}}\right)\mathcal{P}(j\leftarrow
i)$ $\displaystyle=\sum_{j=1,j\neq i}^{n}x_{j}q_{i|j}+\sum_{j=1,j\neq
i}^{n}x_{j}\sum_{\sum_{i^{\prime}=1}^{n}k_{i^{\prime}}=k}\frac{k!}{\prod_{i^{\prime}=1}^{n}k_{i^{\prime}}!}\left(\prod_{i^{\prime}=1}^{n}{q_{i^{\prime}|j}}^{k_{i^{\prime}}}\right)\frac{k_{i}}{k}\left(\pi_{i|j}^{\mathbf{k}}-\sum_{i^{\prime}=1}^{n}\frac{k_{i^{\prime}}}{k}\pi_{i^{\prime}|j}^{\mathbf{k}}\right)\delta+\mathcal{O}(\delta^{2}).$
(S128)

###  A.2 The decrease of $i$-players

The decrease of $i$-players happens when a focal $i$-player is selected to
update its strategy and the player who takes the position is not an
$i$-player. Unlike pairwise comparison, where the focal individual has a
probability of keeping its own strategy,
$\sum_{j=1}^{n}\mathcal{P}(i\leftarrow j)\neq 1$, here, under death-birth, the
focal individual must adopt the strategy of a neighbor,
$\sum_{j=1}^{n}\mathcal{P}(i\leftarrow j)=1$. Therefore, given the focal
$i$-player’s neighbor configuration $\mathbf{k}$, the probability that the
player who takes the position is not an $i$-player can be written as

$1-\mathcal{P}(i\leftarrow
i)=1-\frac{k_{i}F_{i|i}^{\mathbf{k}}}{\sum_{i^{\prime}=1}^{n}k_{i^{\prime}}F_{i^{\prime}|i}^{\mathbf{k}}}=\frac{k-k_{i}}{k}-\frac{k_{i}}{k}\left(\pi_{i|i}^{\mathbf{k}}-\sum_{i^{\prime}=1}^{n}\frac{k_{i^{\prime}}}{k}\pi_{i^{\prime}|i}^{\mathbf{k}}\right)\delta+\mathcal{O}(\delta^{2}).$
(S129)

Applying it to all possibilities for the neighbor configuration $\mathbf{k}$
after selecting a focal $i$-player with probability $x_{i}$, we obtain the
probability that the number of $i$-players decreases by 1 during a unit time
step,

$\displaystyle\mathcal{P}\left(\Delta x_{i}=-\frac{1}{N}\right)$
$\displaystyle=x_{i}\sum_{\sum_{i^{\prime}=1}^{n}k_{i^{\prime}}=k}\frac{k!}{\prod_{i^{\prime}=1}^{n}k_{i^{\prime}}!}\left(\prod_{i^{\prime}=1}^{n}{q_{i^{\prime}|i}}^{k_{i^{\prime}}}\right)\left[1-\mathcal{P}(i\leftarrow
i)\right]$
$\displaystyle=x_{i}(1-q_{i|i})-x_{i}\sum_{\sum_{i^{\prime}=1}^{n}k_{i^{\prime}}=k}\frac{k!}{\prod_{i^{\prime}=1}^{n}k_{i^{\prime}}!}\left(\prod_{i^{\prime}=1}^{n}{q_{i^{\prime}|i}}^{k_{i^{\prime}}}\right)\frac{k_{i}}{k}\left(\pi_{i|i}^{\mathbf{k}}-\sum_{i^{\prime}=1}^{n}\frac{k_{i^{\prime}}}{k}\pi_{i^{\prime}|i}^{\mathbf{k}}\right)\delta+\mathcal{O}(\delta^{2}).$
(S130)

###  A.3 The replicator equation

The instant change in the proportion $x_{i}$ of $i$-players consists of the
increase and decrease of $i$-players. Applying Eqs. (A.1) and (A.2), and
considering that a full Monte Carlo step contains $N$ elementary steps, we
have

$\displaystyle\dot{x}_{i}$
$\displaystyle=N\times\left\\{\frac{1}{N}\mathcal{P}\left(\Delta
x_{i}=\frac{1}{N}\right)+\left(-\frac{1}{N}\right)\mathcal{P}\left(\Delta
x_{i}=-\frac{1}{N}\right)\right\\}$
$\displaystyle=\sum_{j=1}^{n}x_{j}\sum_{\sum_{i^{\prime}=1}^{n}k_{i^{\prime}}=k}\frac{k!}{\prod_{i^{\prime}=1}^{n}k_{i^{\prime}}!}\left(\prod_{i^{\prime}=1}^{n}{q_{{i^{\prime}}|j}}^{k_{i^{\prime}}}\right)\frac{k_{i}}{k}\left(\pi_{i|j}^{\mathbf{k}}-\sum_{i^{\prime}=1}^{n}\frac{k_{i^{\prime}}}{k}\pi_{i^{\prime}|j}^{\mathbf{k}}\right)\delta+\mathcal{O}(\delta^{2}).$
(S131)

In Eq. (A.3), the $\delta^{0}$ term has been eliminated: applying Eq. (S1.1),
we know $x_{j}q_{i|j}=p_{ij}=x_{i}q_{j|i}$, such that the $\delta^{0}$ term in
Eq. (A.1) can be expressed as $\sum_{j=1,j\neq
i}^{n}x_{j}q_{i|j}=\sum_{j=1,j\neq i}^{n}x_{i}q_{j|i}$; similarly, we know
$1-q_{i|i}=\sum_{j=1,j\neq i}^{n}q_{j|i}$ by Eqs. (S1.1), which rewrites the
$\delta^{0}$ term in Eq. (A.2) as $x_{i}(1-q_{i|i})=x_{i}\sum_{j=1,j\neq
i}^{n}q_{j|i}$, equal to the one in Eq. (A.1).

The $\delta^{0}$ term being eliminated, the instant change in $x_{i}$ happens
on the order of $\delta^{1}$. Meanwhile, the instant change in $q_{i|j}$ for
$i,j=1,2,\dots,n$ happens on the order of $\delta^{0}$ since the $\delta^{0}$
is non-zero (see Appendix B.2). That is, the change in $q_{i|j}$ is much
faster than $x_{i}$, so that $x_{i}$ changes on the basis of $q_{i|j}$
achieving equilibrium. According to Appendix B.2, we have the following
solution when $q_{i|j}$ achieves stability.

$q_{i|j}=\begin{cases}\displaystyle{\frac{k-2}{k-1}x_{i}},&j\neq i,\\\ \\\
\displaystyle{\frac{k-2}{k-1}x_{i}+\frac{1}{k-1}},&j=i.\end{cases}$ (S132)

Eq. (S132) is independently derived under death-birth but appear the same as
the ones under pairwise comparison.

The primitive replicator equation of $x_{i}$ under death-birth is Eq. (A.3),
where $\mathcal{O}(\delta^{2})=0$, $\pi_{j}^{\mathbf{k}}$ and
$\pi_{i}^{\mathbf{k}}$ are given by Eq. (S3), $\pi_{i|j}^{\mathbf{k}}$ and
$\pi_{j|i}^{\mathbf{k}}$ are given by Eq. (S5), and $q_{{i^{\prime}}|j}$ is
given by Eq. (S132). The degrees of freedom of the replicator dynamics system
are $n-1$, represented by independent variables $x_{i}$ ($i=1,2,\dots,n$,
cancel one of them by $\sum_{i}^{n}x_{i}=1$).

###  A.4 Simplification and discussion

Next, we utilize the useful Theorem 1 to simplify the replicator equation
under death-birth and discuss it.

####  A.4.1 The accumulated payoff level

First, we study the accumulated payoff level by keeping expected accumulated
payoffs in the replicator equation. Using Theorem 1, we can calculate Eq.
(A.3) as follows.

$\displaystyle\dot{x}_{i}=$ $\displaystyle\leavevmode\nobreak\
\delta\sum_{j=1}^{n}x_{j}\sum_{\sum_{i^{\prime}=1}^{n}k_{i^{\prime}}=k}\frac{k!}{\prod_{i^{\prime}=1}^{n}k_{i^{\prime}}!}\left(\prod_{i^{\prime}=1}^{n}{q_{{i^{\prime}}|j}}^{k_{i^{\prime}}}\right)\frac{k_{i}}{k}\left(\pi_{i|j}^{\mathbf{k}}-\sum_{i^{\prime}=1}^{n}\frac{k_{i^{\prime}}}{k}\pi_{i^{\prime}|j}^{\mathbf{k}}\right)$
$\displaystyle=$ $\displaystyle\leavevmode\nobreak\
\delta\sum_{j=1}^{n}x_{j}q_{i|j}\sum_{\sum_{i^{\prime}=1}^{n}k_{i^{\prime}}=k-1}\frac{(k-1)!}{\prod_{i^{\prime}=1}^{n}k_{i^{\prime}}!}\left(\prod_{i^{\prime}=1}^{n}{q_{{i^{\prime}}|j}}^{k_{i^{\prime}}}\right)\pi_{i|j}^{\mathbf{k}_{+i}}$
$\displaystyle\leavevmode\nobreak\
-\delta\sum_{j=1}^{n}x_{j}q_{i|j}\sum_{\sum_{i^{\prime}=1}^{n}k_{i^{\prime}}=k-1}\frac{(k-1)!}{\prod_{i^{\prime}=1}^{n}k_{i^{\prime}}!}\left(\prod_{i^{\prime}=1}^{n}{q_{{i^{\prime}}|j}}^{k_{i^{\prime}}}\right)\left(\frac{k_{i}+1}{k}\pi_{i|j}^{\mathbf{k}_{+i}}+\sum_{i^{\prime}=1,i^{\prime}\neq
i}^{n}\frac{k_{i^{\prime}}}{k}\pi_{i^{\prime}|j}^{\mathbf{k}_{+i}}\right)$
$\displaystyle=$ $\displaystyle\leavevmode\nobreak\
\delta\sum_{j=1}^{n}x_{j}q_{i|j}\sum_{\sum_{i^{\prime}=1}^{n}k_{i^{\prime}}=k-1}\frac{(k-1)!}{\prod_{i^{\prime}=1}^{n}k_{i^{\prime}}!}\left(\prod_{i^{\prime}=1}^{n}{q_{{i^{\prime}}|j}}^{k_{i^{\prime}}}\right)\pi_{i|j}^{\mathbf{k}_{+i}}$
$\displaystyle\leavevmode\nobreak\
-\delta\sum_{j=1}^{n}x_{j}q_{i|j}\sum_{\sum_{i^{\prime}=1}^{n}k_{i^{\prime}}=k-1}\frac{(k-1)!}{\prod_{i^{\prime}=1}^{n}k_{i^{\prime}}!}\left(\prod_{i^{\prime}=1}^{n}{q_{{i^{\prime}}|j}}^{k_{i^{\prime}}}\right)\left(\frac{1}{k}\pi_{i|j}^{\mathbf{k}_{+i}}+\sum_{i^{\prime}=1}^{n}\frac{k_{i^{\prime}}}{k}\pi_{i^{\prime}|j}^{\mathbf{k}_{+i}}\right)$
$\displaystyle=$ $\displaystyle\leavevmode\nobreak\
\frac{\delta(k-1)}{k}x_{i}\sum_{j=1}^{n}q_{j|i}\sum_{\sum_{i^{\prime}=1}^{n}k_{i^{\prime}}=k-1}\frac{(k-1)!}{\prod_{i^{\prime}=1}^{n}k_{i^{\prime}}!}\left(\prod_{i^{\prime}=1}^{n}{q_{{i^{\prime}}|j}}^{k_{i^{\prime}}}\right)\pi_{i|j}^{\mathbf{k}_{+i}}$
$\displaystyle\leavevmode\nobreak\
-\frac{\delta(k-1)}{k}x_{i}\sum_{j=1}^{n}q_{j|i}\sum_{\sum_{i^{\prime}=1}^{n}k_{i^{\prime}}=k-2}\frac{(k-2)!}{\prod_{i^{\prime}=1}^{n}k_{i^{\prime}}!}\left(\prod_{i^{\prime}=1}^{n}{q_{{i^{\prime}}|j}}^{k_{i^{\prime}}}\right)\sum_{i^{\prime}=1}^{n}q_{i^{\prime}|j}\pi_{i^{\prime}|j}^{\mathbf{k}_{+i,+i^{\prime}}}.$
(S133)

Using the aforementioned notations for
$\langle\pi_{i|j}^{\mathbf{k}_{+i}}\rangle$ in Eq. (S16), as well as its
extension for
$\langle\pi_{i^{\prime}|j}^{\mathbf{k}_{+i,+i^{\prime}}}\rangle$,

$\langle\pi_{i^{\prime}|j}^{\mathbf{k}_{+i,+i^{\prime}}}\rangle=\sum_{\sum_{i^{\prime}=1}^{n}k_{i^{\prime}}=k-2}\frac{(k-2)!}{\prod_{i^{\prime}=1}^{n}k_{i^{\prime}}!}\left(\prod_{i^{\prime}=1}^{n}{q_{{i^{\prime}}|j}}^{k_{i^{\prime}}}\right)\pi_{i^{\prime}|j}^{\mathbf{k}_{+i,+i^{\prime}}},$
(S134)

we can write Eq. (A.4.1) as

$\displaystyle\dot{x}_{i}=$ $\displaystyle\leavevmode\nobreak\
\frac{\delta(k-1)}{k}x_{i}\sum_{j=1}^{n}q_{j|i}\left(\langle\pi_{i|j}^{\mathbf{k}_{+i}}\rangle-\sum_{i^{\prime}=1}^{n}q_{i^{\prime}|j}\langle\pi_{i^{\prime}|j}^{\mathbf{k}_{+i,+i^{\prime}}}\rangle\right)$
$\displaystyle=$ $\displaystyle\leavevmode\nobreak\
\frac{\delta(k-1)}{k}x_{i}\left(\langle\pi_{i}^{\mathbf{k}}\rangle-\sum_{j=1}^{n}\sum_{i^{\prime}=1}^{n}q_{i^{\prime}|j}q_{j|i}\langle\pi_{i^{\prime}|j}^{\mathbf{k}_{+i,+i^{\prime}}}\rangle\right).$
(S135)

In Eq. (A.4.1), we have used the relation
$\sum_{j=1}^{n}q_{j|i}\langle\pi_{i|j}^{\mathbf{k}_{+i}}\rangle=\langle\pi_{i}^{\mathbf{k}}\rangle$
as suggested by Theorem 3.

Again, the replicator equation at the accumulated payoff level given by Eq.
(A.4.1) bears an intuitive understanding if we use the following concepts:

  * •

$\pi_{i}^{(0)}=\langle\pi_{i}^{\mathbf{k}}\rangle$, the expected accumulated
payoff of the $i$-player itself (zero-step away on the graph).

  * •

$\pi_{i}^{(2)}=\sum_{j=1}^{n}\sum_{i^{\prime}=1}^{n}q_{i^{\prime}|j}q_{j|i}\langle\pi_{i^{\prime}|j}^{\mathbf{k}_{+i,+i^{\prime}}}\rangle$,
the expected accumulated payoff of the $i$-player’s second-order neighbors
(two-step away on the graph). Intuitively, $q_{i^{\prime}|j}q_{j|i}$ is used
to find two-step away $i^{\prime}$-players away from the $i$-player. The first
step walks to a $j$-player, whose neighbor configuration is
$\mathbf{k}_{+i,+i^{\prime}}$. The second step walks to an
$i^{\prime}$-player.

Using these concepts, we know that $\dot{x}_{i}\sim
x_{i}(\pi_{i}^{(0)}-\pi_{i}^{(2)})$. Under death-birth, the reproduction rate
of $i$-players depends on how much their expected accumulated payoff higher
than the one of neighbors’ neighbors. The essence of replicator dynamics
$\dot{x}_{i}$ is the competition between oneself and its second-order
neighbors. This is also consistent with the result obtained by the identity-
by-descent idea [18, 20], and here we further generalize it to $n$-strategy
multiplayer systems.

####  A.4.2 The single game level

Next, we further divide the replicator equation into the single game level,
stressing the convenience for actual calculation.

Similarly, we use the notation of $\langle a_{i|\mathbf{k}}\rangle_{j}$
defined in Eq. (S21) to represent the payoff in a single game. We know the
expression of $\langle\pi_{i}^{\mathbf{k}}\rangle$ from Eq. (S2.4.2), which is
derived from the payoff calculation stage, independent of strategy update
rules. Furthermore, we expand
$\langle\pi_{i^{\prime}|j}^{\mathbf{k}_{+i,+i^{\prime}}}\rangle$ by Eq. (S5)
and obtain

$\displaystyle\langle\pi_{i^{\prime}|j}^{\mathbf{k}_{+i,+i^{\prime}}}\rangle=$
$\displaystyle\sum_{\sum_{i^{\prime}=1}^{n}k_{i^{\prime}}=k-2}\frac{(k-2)!}{\prod_{i^{\prime}=1}^{n}k_{i^{\prime}}!}\left(\prod_{i^{\prime}=1}^{n}{q_{{i^{\prime}}|j}}^{k_{i^{\prime}}}\right)\Bigg{\\{}a_{i^{\prime}|\mathbf{k}_{+i,+j}}+\sum_{\sum_{l=1}^{n}k^{\prime}_{l}=k-1}\frac{(k-1)!}{\prod_{l=1}^{n}k^{\prime}_{l}!}\left(\prod_{l=1}^{n}{q_{l|i^{\prime}}}^{k^{\prime}_{l}}\right)\Bigg{[}a_{i^{\prime}|\mathbf{k}^{\prime}_{+j}}$
$\displaystyle+\sum_{l=1}^{n}k^{\prime}_{l}\sum_{\sum_{\ell=1}^{n}k^{\prime\prime}_{\ell}=k-1}\frac{(k-1)!}{\prod_{\ell=1}^{n}k^{\prime\prime}_{\ell}!}\left(\prod_{\ell=1}^{n}{q_{\ell|l}}^{k^{\prime\prime}_{\ell}}\right)a_{i^{\prime}|\mathbf{k}^{\prime\prime}_{+l}}\Bigg{]}\Bigg{\\}}$
$\displaystyle=$
$\displaystyle\sum_{\sum_{i^{\prime}=1}^{n}k_{i^{\prime}}=k-2}\frac{(k-2)!}{\prod_{i^{\prime}=1}^{n}k_{i^{\prime}}!}\left(\prod_{i^{\prime}=1}^{n}{q_{{i^{\prime}}|j}}^{k_{i^{\prime}}}\right)a_{i^{\prime}|\mathbf{k}_{+i,+j}}+\sum_{\sum_{l=1}^{n}k^{\prime}_{l}=k-1}\frac{(k-1)!}{\prod_{l=1}^{n}k^{\prime}_{l}!}\left(\prod_{l=1}^{n}{q_{l|i^{\prime}}}^{k^{\prime}_{l}}\right)a_{i^{\prime}|\mathbf{k}^{\prime}_{+j}}$
$\displaystyle+(k-1)\sum_{l=1}^{n}q_{l|i^{\prime}}\sum_{\sum_{\ell=1}^{n}k^{\prime\prime}_{\ell}=k-1}\frac{(k-1)!}{\prod_{\ell=1}^{n}k^{\prime\prime}_{\ell}!}\left(\prod_{\ell=1}^{n}{q_{\ell|l}}^{k^{\prime\prime}_{\ell}}\right)a_{i^{\prime}|\mathbf{k}^{\prime\prime}_{+l}}$
$\displaystyle=$ $\displaystyle\leavevmode\nobreak\ \langle
a_{i^{\prime}|\mathbf{k}_{+i,+j}}\rangle_{j}+\langle
a_{i^{\prime}|\mathbf{k}_{+j}}\rangle_{i^{\prime}}+(k-1)\sum_{l=1}^{n}q_{l|i^{\prime}}\langle
a_{i^{\prime}|\mathbf{k}_{+l}}\rangle_{l}.$ (S136)

Applying Eqs. (S2.4.2) and (A.4.2) to Eq. (A.4.1), we have

$\dot{x}_{i}=\frac{\delta(k-1)}{k}x_{i}\left[\left(\langle
a_{i|\mathbf{k}}\rangle_{i}+k\sum_{j=1}^{n}q_{j|i}\langle
a_{i|\mathbf{k}_{+j}}\rangle_{j}\right)-\sum_{j=1}^{n}\sum_{i^{\prime}=1}^{n}q_{i^{\prime}|j}q_{j|i}\left(\langle
a_{i^{\prime}|\mathbf{k}_{+j}}\rangle_{i^{\prime}}+\langle
a_{i^{\prime}|\mathbf{k}_{+i,+j}}\rangle_{j}+(k-1)\sum_{l=1}^{n}q_{l|i^{\prime}}\langle
a_{i^{\prime}|\mathbf{k}_{+l}}\rangle_{l}\right)\right].$ (S137)

To decrease the number of elements, we transform $\langle
a_{i|\mathbf{k}}\rangle_{i}$ into $\sum_{j=1}^{n}q_{j|i}\langle
a_{i|\mathbf{k}_{+j}}\rangle_{i}$ according to Theorem 4 and rewrite Eq.
(S137) as

$\dot{x}_{i}=\frac{\delta(k-1)}{k}x_{i}\sum_{j=1}^{n}q_{j|i}\left(\langle
a_{i|\mathbf{k}_{+j}}\rangle_{i}+k\langle
a_{i|\mathbf{k}_{+j}}\rangle_{j}-\sum_{i^{\prime}=1}^{n}q_{i^{\prime}|j}\left(\langle
a_{i^{\prime}|\mathbf{k}_{+j}}\rangle_{i^{\prime}}+\langle
a_{i^{\prime}|\mathbf{k}_{+i,+j}}\rangle_{j}+(k-1)\sum_{l=1}^{n}q_{l|i^{\prime}}\langle
a_{i^{\prime}|\mathbf{k}_{+l}}\rangle_{l}\right)\right).$ (S138)

We do not bother giving the general expression of Eq. (S138) with all
$q_{j|i}$ quantities transformed into $x_{j}$, since the results would be
complicated and unruly. Please note that the calculation inside each
$\langle\cdot\rangle$ also contains $q_{j|i}$ quantities and we always need to
transform them manually in applications.

According to Eq. (S138), we can attribute everything about
$\langle\cdot\rangle$ into three types, the ‘$\langle
a_{i|\mathbf{k}_{+j}}\rangle_{i}$ type’, the ‘$\langle
a_{i|\mathbf{k}_{+j}}\rangle_{j}$ type’, and the ‘$\langle
a_{i^{\prime}|\mathbf{k}_{+i,+j}}\rangle_{j}$ type’. The $\langle
a_{i|\mathbf{k}_{+j}}\rangle_{i}$ and $\langle
a_{i|\mathbf{k}_{+j}}\rangle_{j}$ types have been listed as matrices
$\left[\langle a_{i|\mathbf{k}_{+j}}\rangle_{i}\right]_{ij}$ and
$\left[\langle a_{i|\mathbf{k}_{+j}}\rangle_{j}\right]_{ij}$ in Eqs. (S27) and
(S28), respectively. The number of elements to calculate manually in the
$\langle a_{i|\mathbf{k}_{+j}}\rangle_{i}$ and $\langle
a_{i|\mathbf{k}_{+j}}\rangle_{j}$ types is $(2n-1)n$ as we have discussed.

The $\langle a_{i^{\prime}|\mathbf{k}_{+i,+j}}\rangle_{j}$ type, however,
contains three dimensions. It can be expressed by totaling $i^{\prime}$
matrices, each with dimensions $i$ and $j$:

  * •

The $\langle a_{i^{\prime}|\mathbf{k}_{+i,+j}}\rangle_{j}$ type:

$\left[\langle
a_{i^{\prime}|\mathbf{k}_{+i,+j}}\rangle_{j}\right]_{ij}=\displaystyle{\begin{pmatrix}\langle
a_{i^{\prime}|\mathbf{k}_{+1,+1}}\rangle_{1}&\langle
a_{i^{\prime}|\mathbf{k}_{+1,+2}}\rangle_{2}&\cdots&\langle
a_{i^{\prime}|\mathbf{k}_{+1,+n}}\rangle_{n}\\\ \langle
a_{i^{\prime}|\mathbf{k}_{+2,+1}}\rangle_{1}&\langle
a_{i^{\prime}|\mathbf{k}_{+2,+2}}\rangle_{2}&\cdots&\langle
a_{i^{\prime}|\mathbf{k}_{+2,+n}}\rangle_{n}\\\ \vdots&\vdots&\ddots&\vdots\\\
\langle a_{i^{\prime}|\mathbf{k}_{+n,+1}}\rangle_{1}&\langle
a_{i^{\prime}|\mathbf{k}_{+n,+2}}\rangle_{2}&\cdots&\langle
a_{i^{\prime}|\mathbf{k}_{+n,+n}}\rangle_{n}\end{pmatrix}},$ (S139)

for $i^{\prime}=1,2,\dots,n$.

There are $n^{3}$ elements in all $\left[\langle
a_{i^{\prime}|\mathbf{k}_{+i,+j}}\rangle_{j}\right]_{ij}$ matrices. Therefore,
the total computation involves $n^{3}+(2n-1)n=(n^{2}+2n-1)n$ elements. The
computational complexity is $\mathrm{O}(n^{3})$, which is feasible within
polynomial time.

####  A.4.3 Special linear system

We can also further simplify the calculation when faced with special payoff
structures. Given a co-player configuration $\mathbf{k}$, a linear payoff
structure can be determined by the coefficient matrix $\mathbf{b}$ and the
constant vector $\mathbf{c}$, that is,
$a_{i|\mathbf{k}}=\sum_{l=1}^{n}b_{il}k_{l}+c_{i}$, as previously shown in
Eqs. (S29)–(S31).

The expressions of $\langle a_{i|\mathbf{k}_{+j}}\rangle_{i}$ and $\langle
a_{i|\mathbf{k}_{+j}}\rangle_{j}$ as functions of $\mathbf{b}$ and
$\mathbf{c}$ has been obtained in Eqs. (S34a) and (S34b), respectively. We
further calculate the elements of the $\langle
a_{i^{\prime}|\mathbf{k}_{+i,+j}}\rangle_{j}$ type,

$\langle
a_{i^{\prime}|\mathbf{k}_{+i,+j}}\rangle_{j}=(k-2)\sum_{l=1}^{n}b_{i^{\prime}l}q_{l|j}+b_{i^{\prime}i}+b_{i^{\prime}j}+c_{i^{\prime}}=\frac{(k-2)^{2}}{k-1}\sum_{l=1}^{n}b_{i^{\prime}l}x_{l}+\frac{2k-3}{k-1}b_{i^{\prime}j}+b_{i^{\prime}i}+c_{i^{\prime}}.$
(S140)

We substitute Eqs. (S34a), (S34b) and (S140) into Eq. (S138). Through a long
but feasible calculation (transforming all $q_{j|i}$ quantities to $x_{j}$
quantities and organizing the result), we obtain

$\displaystyle\dot{x}_{i}=$ $\displaystyle\leavevmode\nobreak\
\frac{\delta(k-2)}{k(k-1)^{2}}x_{i}\Bigg{\\{}\frac{(k^{2}-2)^{2}}{k}(\bar{\pi}_{i}-\bar{\pi})+\frac{3k^{2}-4}{k}\left((kb_{ii}+c_{i})-\sum_{j=1}^{n}x_{j}(kb_{jj}+c_{j})\right)$
$\displaystyle\leavevmode\nobreak\
-(k^{2}+2k-4)\sum_{j=1}^{n}x_{j}\left(b_{ji}-\sum_{l=1}^{n}x_{l}b_{jl}\right)\Bigg{\\}},$
(S141)

where $\bar{\pi}_{i}=k\sum_{l=1}^{n}x_{l}b_{il}+c_{i}$ is the expected payoff
of $i$-players in a well-mixed population, and
$\bar{\pi}=\sum_{i=1}^{n}x_{i}\pi_{i}$ is the expected payoff of all
individuals in a well-mixed population.

We know that the replicator equation in a well-mixed population is
$\dot{x}_{i}=x_{i}(\bar{\pi}_{i}-\bar{\pi})$. In this way, Eq. (A.4.3) clearly
showed the additional terms brought by death-birth in a structured population
compared to the well-mixed population.

###  A.5 Comparison with two-strategy dynamics

Let us discuss on how our replicator equations under death-birth for
$n$-strategy systems reduce to the 2-strategy system. Similarly, the
configuration $\mathbf{k}=(k_{1},k_{2})=(k_{1},k-k_{1})$ can be represented by
only $k_{1}$, where $k_{1}=0,1,\dots,k$.

The $n$-strategy system simplified to the accumulated payoff level shown in
Eq. (A.4.1) is for comparison with the identity-by-decent idea by Allen and
Nowak [18], that the essence of death-birth is the competition from the
individuals from two-step away. However, in this section, we want to compare
with the 2-strategy system for general multiplayer games in structured
populations, which was studied by Li et al. [44]. In this way, we need to
reorganize the results in Eq. (A.4.1) first.

According to Theorem 5, we have
$\langle\pi_{i}^{\mathbf{k}}\rangle=\sum_{j=1}^{n}\sum_{i^{\prime}=1}^{n}q_{i^{\prime}|j}q_{j|i}\langle\pi_{i|j}^{\mathbf{k}_{+i,+i^{\prime}}}\rangle$.
Substituting this into Eq. (A.4.1), we have

$\dot{x}_{i}=\frac{\delta(k-1)}{k}x_{i}\sum_{j=1}^{n}\sum_{i^{\prime}=1}^{n}q_{i^{\prime}|j}q_{j|i}\left(\langle\pi_{i|j}^{\mathbf{k}_{+i,+i^{\prime}}}\rangle-\langle\pi_{i^{\prime}|j}^{\mathbf{k}_{+i,+i^{\prime}}}\rangle\right).$
(S142)

Further replacing all $q_{j|i}$ quantities with $x_{j}$ quantities by Eq.
(S132) and organizing the results, we have

$\dot{x}_{i}=\frac{\delta(k-2)}{k(k-1)}x_{i}\sum_{j=1}^{n}x_{j}\left[\left(\langle\pi_{i|j}^{\mathbf{k}_{+i,+j}}\rangle-\langle\pi_{j|j}^{\mathbf{k}_{+i,+j}}\rangle\right)+\left(\langle\pi_{i|i}^{\mathbf{k}_{+i,+j}}\rangle-\langle\pi_{j|i}^{\mathbf{k}_{+i,+j}}\rangle\right)+(k-2)\sum_{i^{\prime}=1}^{n}x_{i^{\prime}}\left(\langle\pi_{i|j}^{\mathbf{k}_{+i,+i^{\prime}}}\rangle-\langle\pi_{i^{\prime}|j}^{\mathbf{k}_{+i,+i^{\prime}}}\rangle\right)\right],$
(S143)

which can be reduced to $2$-strategy system depicted by $\dot{x}_{1}$ with the
consideration of $x_{2}=1-x_{1}$,

$\displaystyle\dot{x}_{1}=$ $\displaystyle\leavevmode\nobreak\
\frac{\delta(k-2)}{k(k-1)}x_{1}(1-x_{1})\Big{\\{}k\left(\langle\pi_{1|2}^{\mathbf{k}_{+1,+2}}\rangle-\langle\pi_{2|2}^{\mathbf{k}_{+1,+2}}\rangle\right)$
$\displaystyle+[(k-2)x_{1}+1]\left[\left(\langle\pi_{1|1}^{\mathbf{k}_{+1,+2}}\rangle-\langle\pi_{2|1}^{\mathbf{k}_{+1,+2}}\rangle\right)-\left(\langle\pi_{1|2}^{\mathbf{k}_{+1,+2}}\rangle-\langle\pi_{2|2}^{\mathbf{k}_{+1,+2}}\rangle\right)\right]\Big{\\}}.$
(S144)

Eq. (A.5), which describes general 2-strategy multiplayer games in structured
populations under the death-birth rule, has the same form as the one given by
Li et al. [44]. As a supplement to their work, here we stress a detail as
indicated in Eq. (A.5): the calculation should be made in the remaining $k-2$
individuals given the existence of at least a pair of 1-player and 2-player.

In particular, let us list the quantities in Eq. (A.5) and clarify the
calculation. Referring to Eq. (S16), we have

$\displaystyle\langle\pi_{1|1}^{\mathbf{k}_{+1,+2}}\rangle$
$\displaystyle=\sum_{k_{1}=0}^{k-2}\frac{(k-2)!}{k_{1}!(k-k_{1}-2)!}{q_{1|1}}^{k_{1}}{q_{2|1}}^{k-k_{1}-2}\pi_{1|1}^{\mathbf{k}_{+1,+2}},$
(S145a) $\displaystyle\langle\pi_{2|1}^{\mathbf{k}_{+1,+2}}\rangle$
$\displaystyle=\sum_{k_{1}=0}^{k-2}\frac{(k-2)!}{k_{1}!(k-k_{1}-2)!}{q_{1|1}}^{k_{1}}{q_{2|1}}^{k-k_{1}-2}\pi_{2|1}^{\mathbf{k}_{+1,+2}},$
(S145b) $\displaystyle\langle\pi_{1|2}^{\mathbf{k}_{+1,+2}}\rangle$
$\displaystyle=\sum_{k_{1}=0}^{k-2}\frac{(k-2)!}{k_{1}!(k-k_{1}-2)!}{q_{1|2}}^{k_{1}}{q_{2|2}}^{k-k_{1}-2}\pi_{1|2}^{\mathbf{k}_{+1,+2}},$
(S145c) $\displaystyle\langle\pi_{2|2}^{\mathbf{k}_{+1,+2}}\rangle$
$\displaystyle=\sum_{k_{1}=0}^{k-2}\frac{(k-2)!}{k_{1}!(k-k_{1}-2)!}{q_{1|2}}^{k_{1}}{q_{2|2}}^{k-k_{1}-2}\pi_{2|2}^{\mathbf{k}_{+1,+2}},$
(S145d)

where, according to Eq. (S5), we have

$\displaystyle\pi_{1|1}^{\mathbf{k}_{+1,+2}}=$
$\displaystyle\leavevmode\nobreak\
a_{1|\mathbf{k}_{+1,+2}}+\sum_{k^{\prime}_{1}=0}^{k-1}\frac{(k-1)!}{k^{\prime}_{1}!(k-k^{\prime}_{1}-1)!}{q_{1|1}}^{k^{\prime}_{1}}{q_{2|1}}^{k-k^{\prime}_{1}-1}\Bigg{(}a_{1|\mathbf{k}^{\prime}_{+1}}+k^{\prime}_{1}\sum_{k^{\prime\prime}_{1}=0}^{k-1}\frac{(k-1)!}{k^{\prime\prime}_{1}!(k-k^{\prime\prime}_{1}-1)!}{q_{1|1}}^{k^{\prime\prime}_{1}}{q_{2|1}}^{k-k^{\prime\prime}_{1}-1}a_{1|\mathbf{k}^{\prime\prime}_{+1}}$
$\displaystyle+(k-k^{\prime}_{1}-1)\sum_{k^{\prime\prime}_{1}=0}^{k-1}\frac{(k-1)!}{k^{\prime\prime}_{1}!(k-k^{\prime\prime}_{1}-1)!}{q_{1|2}}^{k^{\prime\prime}_{1}}{q_{2|2}}^{k-k^{\prime\prime}_{1}-1}a_{1|\mathbf{k}^{\prime\prime}_{+2}}\Bigg{)},$
(S146a) $\displaystyle\pi_{2|1}^{\mathbf{k}_{+1,+2}}=$
$\displaystyle\leavevmode\nobreak\
a_{2|\mathbf{k}_{+1,+1}}+\sum_{k^{\prime}_{1}=0}^{k-1}\frac{(k-1)!}{k^{\prime}_{1}!(k-k^{\prime}_{1}-1)!}{q_{1|2}}^{k^{\prime}_{1}}{q_{2|2}}^{k-k^{\prime}_{1}-1}\Bigg{(}a_{2|\mathbf{k}^{\prime}_{+1}}+k^{\prime}_{1}\sum_{k^{\prime\prime}_{1}=0}^{k-1}\frac{(k-1)!}{k^{\prime\prime}_{1}!(k-k^{\prime\prime}_{1}-1)!}{q_{1|1}}^{k^{\prime\prime}_{1}}{q_{2|1}}^{k-k^{\prime\prime}_{1}-1}a_{2|\mathbf{k}^{\prime\prime}_{+1}}$
$\displaystyle+(k-k^{\prime}_{1}-1)\sum_{k^{\prime\prime}_{1}=0}^{k-1}\frac{(k-1)!}{k^{\prime\prime}_{1}!(k-k^{\prime\prime}_{1}-1)!}{q_{1|2}}^{k^{\prime\prime}_{1}}{q_{2|2}}^{k-k^{\prime\prime}_{1}-1}a_{2|\mathbf{k}^{\prime\prime}_{+2}}\Bigg{)},$
(S146b) $\displaystyle\pi_{1|2}^{\mathbf{k}_{+1,+2}}=$
$\displaystyle\leavevmode\nobreak\
a_{1|\mathbf{k}_{+2,+2}}+\sum_{k^{\prime}_{1}=0}^{k-1}\frac{(k-1)!}{k^{\prime}_{1}!(k-k^{\prime}_{1}-1)!}{q_{1|1}}^{k^{\prime}_{1}}{q_{2|1}}^{k-k^{\prime}_{1}-1}\Bigg{(}a_{1|\mathbf{k}^{\prime}_{+2}}+k^{\prime}_{1}\sum_{k^{\prime\prime}_{1}=0}^{k-1}\frac{(k-1)!}{k^{\prime\prime}_{1}!(k-k^{\prime\prime}_{1}-1)!}{q_{1|1}}^{k^{\prime\prime}_{1}}{q_{2|1}}^{k-k^{\prime\prime}_{1}-1}a_{1|\mathbf{k}^{\prime\prime}_{+1}}$
$\displaystyle+(k-k^{\prime}_{1}-1)\sum_{k^{\prime\prime}_{1}=0}^{k-1}\frac{(k-1)!}{k^{\prime\prime}_{1}!(k-k^{\prime\prime}_{1}-1)!}{q_{1|2}}^{k^{\prime\prime}_{1}}{q_{2|2}}^{k-k^{\prime\prime}_{1}-1}a_{1|\mathbf{k}^{\prime\prime}_{+2}}\Bigg{)},$
(S146c) $\displaystyle\pi_{2|2}^{\mathbf{k}_{+1,+2}}=$
$\displaystyle\leavevmode\nobreak\
a_{2|\mathbf{k}_{+1,+2}}+\sum_{k^{\prime}_{1}=0}^{k-1}\frac{(k-1)!}{k^{\prime}_{1}!(k-k^{\prime}_{1}-1)!}{q_{1|2}}^{k^{\prime}_{1}}{q_{2|2}}^{k-k^{\prime}_{1}-1}\Bigg{(}a_{2|\mathbf{k}^{\prime}_{+2}}+k^{\prime}_{1}\sum_{k^{\prime\prime}_{1}=0}^{k-1}\frac{(k-1)!}{k^{\prime\prime}_{1}!(k-k^{\prime\prime}_{1}-1)!}{q_{1|1}}^{k^{\prime\prime}_{1}}{q_{2|1}}^{k-k^{\prime\prime}_{1}-1}a_{2|\mathbf{k}^{\prime\prime}_{+1}}$
$\displaystyle+(k-k^{\prime}_{1}-1)\sum_{k^{\prime\prime}_{1}=0}^{k-1}\frac{(k-1)!}{k^{\prime\prime}_{1}!(k-k^{\prime\prime}_{1}-1)!}{q_{1|2}}^{k^{\prime\prime}_{1}}{q_{2|2}}^{k-k^{\prime\prime}_{1}-1}a_{2|\mathbf{k}^{\prime\prime}_{+2}}\Bigg{)}.$
(S146d)

In this way, the calculation of general 2-strategy multiplayer games given by
Eq. (A.5) is completed.

Similarly, after presenting the 2-strategy case at the accumulated payoff
level, we can further compute the results, which are divided into the single
game level. Under death-birth, however, the results would seem unruly and not
easier than Eq. (A.5). The 2-strategy case at the single game level can be
obtained by applying $n=2$ to Eq. (S138), which leads to an over-complicated
result. In conclusion, we always need to treat the payoff structure of
different models case by case, and using Eq. (A.5) is an acceptable choice
that we have.

In the case of special linear systems, however, we can refer to a well-
organized general equation. Applying $n=2$ to Eq. (A.4.3), we have

$\displaystyle\dot{x}_{1}=$ $\displaystyle\leavevmode\nobreak\
\frac{\delta(k-2)}{k(k-1)^{2}}\Bigg{\\{}\frac{(k^{2}-2)^{2}}{k}x_{1}(\bar{\pi}_{i}-\bar{\pi})+\frac{3k^{2}-4}{k}x_{1}(1-x_{1})\left[k(b_{11}-b_{22})+c_{1}-c_{2}\right]$
$\displaystyle-(k^{2}+2k-4)x_{1}(1-x_{1})\left[(b_{11}-b_{12})x_{1}+(b_{21}-b_{22})(1-x_{1})\right]\Bigg{\\}},$
(S147)

where $b_{ij}$ and $c_{i}$ are given by Eq. (S29) to express the linear payoff
structure $a_{i|\mathbf{k}}=\sum_{l=1}^{n}b_{il}k_{l}+c_{i}$, and
$x_{1}(\bar{\pi}_{1}-\bar{\pi})$ corresponds to the replicator dynamics in a
well-mixed population, which can be obtained by existing knowledge of well-
mixed populations.

##  Appendix B Edge dynamics

As we know from Eq. (S1.1) in section S1.1, the system can be described by
$x_{i}$ and $q_{j|i}$. While we study the dynamics of $x_{i}$ in previous
sections, the dynamics of $q_{j|i}$ is further supplemented here. In
particular, we show how to obtain the relation between $x_{i}$ and $q_{j|i}$
as shown in Eqs. (S13) and (S132), for both pairwise comparison and death-
birth. The same results have been declared in previous literature [33], but
the details of the calculation process were unclear.

###  B.1 Pairwise comparison

First, we study edge dynamics under pairwise comparison to explain the origin
of Eq. (S13) presented in section S2.3.

####  B.1.1 The increase of $il$-edges

The number of $il$-edges in the system increases by $k_{l}$ when a $j$-player
is replaced by an $i$-player ($j\neq i$)—the new $i$-player connects to the
previous $j$-player’s neighbors $l$, thus increasing the number of $il$-edges
by $k_{l}$. Since the total number of edges in the system is $kN/2$, the
frequency of $il$-edges, $p_{il}$, is increased by $k_{l}/(kN/2)=2k_{l}/(kN)$.
For a possible neighbor configuration $\mathbf{k}$ of the $j$-player, we go
through all $j\neq i$, obtaining the probability that the configuration
$\mathbf{k}$ is found around the $j$-player and thus the frequency of
$il$-edges increases by $2k_{l}/(kN)$:

$\mathcal{P}\left(\Delta p_{il}=\frac{2k_{l}}{kN}\right)=\sum_{j=1,j\neq
i}^{n}x_{j}\frac{k!}{\prod_{i^{\prime}=1}^{n}k_{i^{\prime}}!}\left(\prod_{i^{\prime}=1}^{n}{q_{i^{\prime}|j}}^{k_{i^{\prime}}}\right)\mathcal{P}(j\leftarrow
i).$ (S148)

Then, we go through all possible configurations $\mathbf{k}$. We put this step
late here because the change in the frequency of $il$-edges, $2k_{l}/(kN)$,
varies in different $\mathbf{k}$. As we see, the expected change in the
frequency of $il$-edges is calculated as:

$\displaystyle\sum_{\sum_{i^{\prime}=1}^{n}k_{i^{\prime}}=k}\frac{2k_{l}}{kN}\mathcal{P}\left(\Delta
p_{il}=\frac{2k_{l}}{kN}\right)$
$\displaystyle=\sum_{\sum_{i^{\prime}=1}^{n}k_{i^{\prime}}=k}\frac{2k_{l}}{kN}\sum_{j=1,j\neq
i}^{n}x_{j}\frac{k!}{\prod_{i^{\prime}=1}^{n}k_{i^{\prime}}!}\left(\prod_{i^{\prime}=1}^{n}{q_{i^{\prime}|j}}^{k_{i^{\prime}}}\right)\mathcal{P}(j\leftarrow
i)$
$\displaystyle=\sum_{\sum_{i^{\prime}=1}^{n}k_{i^{\prime}}=k}\frac{2k_{l}}{kN}\sum_{j=1,j\neq
i}^{n}x_{j}\frac{k!}{\prod_{i^{\prime}=1}^{n}k_{i^{\prime}}!}\left(\prod_{i^{\prime}=1}^{n}{q_{i^{\prime}|j}}^{k_{i^{\prime}}}\right)\left(\frac{k_{i}}{2k}+\mathcal{O}(\delta)\right)$
$\displaystyle=\begin{cases}\displaystyle{\sum_{j=1,j\neq
i}^{n}x_{j}\frac{q_{i|j}[(k-1)q_{i|j}+1]}{kN}+\mathcal{O}(\delta)},&\mbox{$i=l$,}\\\
\displaystyle{\sum_{j=1,j\neq
i}^{n}x_{j}\frac{(k-1)q_{l|j}q_{i|j}}{kN}+\mathcal{O}(\delta)},&\mbox{$i\neq
l$.}\end{cases}$
$\displaystyle=\frac{x_{i}}{kN}\begin{cases}\displaystyle{\sum_{j=1,j\neq
i}^{n}q_{j|i}[(k-1)q_{i|j}+1]+\mathcal{O}(\delta)},&\mbox{$i=l$,}\\\
\displaystyle{\sum_{j=1,j\neq
i}^{n}(k-1)q_{l|j}q_{j|i}+\mathcal{O}(\delta)},&\mbox{$i\neq l$.}\end{cases}$
(S149)

In Eq. (B.1.1), $\mathcal{P}(j\leftarrow i)$ refers to Eq. (S8). We only take
the first-order Taylor expansion, because the first-order term would be non-
zero as we will see in Eq. (B.1.3), thus not necessary to study the second-
order term.

####  B.1.2 The decrease of $il$-edges

The number of $il$-edges in the system decreases by $k_{l}$ when an $i$-player
is replaced by a $j$-player ($j\neq i$). The died $i$-player disconnects from
the previous neighbors $l$, thus decreasing the number of $il$-edges by
$k_{l}$. In other words, the frequency of $il$-edges $p_{il}$ is decreased by
$k_{l}/(kN/2)=2k_{l}/(kN)$. For a possible neighbor configuration $\mathbf{k}$
of the $i$-player, we go through all possible $j\neq i$, obtaining the
probability that the configuration $\mathbf{k}$ is found around the $i$-player
and thus the frequency of $il$-edges decreases by $2k_{l}/(kN)$:

$\mathcal{P}\left(\Delta p_{il}=-\frac{2k_{l}}{kN}\right)=x_{i}\sum_{j=1,j\neq
i}^{n}\frac{k!}{\prod_{i^{\prime}=1}^{n}k_{i^{\prime}}!}\left(\prod_{i^{\prime}=1}^{n}{q_{i^{\prime}|i}}^{k_{i^{\prime}}}\right)\mathcal{P}(i\leftarrow
j).$ (S150)

Then, we go through all possible configurations $\mathbf{k}$. The expected
change in the frequency of $il$-edges is calculated as:

$\displaystyle\sum_{\sum_{i^{\prime}=1}^{n}k_{i^{\prime}}=k}\left(-\frac{2k_{l}}{kN}\right)\mathcal{P}\left(\Delta
p_{il}=-\frac{2k_{l}}{kN}\right)$
$\displaystyle=\sum_{\sum_{i^{\prime}=1}^{n}k_{i^{\prime}}=k}\left(-\frac{2k_{l}}{kN}\right)x_{i}\sum_{j=1,j\neq
i}^{n}\frac{k!}{\prod_{i^{\prime}=1}^{n}k_{i^{\prime}}!}\left(\prod_{i^{\prime}=1}^{n}{q_{i^{\prime}|i}}^{k_{i^{\prime}}}\right)\mathcal{P}(i\leftarrow
j)$
$\displaystyle=\sum_{\sum_{i^{\prime}=1}^{n}k_{i^{\prime}}=k}\left(-\frac{2k_{l}}{kN}\right)x_{i}\sum_{j=1,j\neq
i}^{n}\frac{k!}{\prod_{i^{\prime}=1}^{n}k_{i^{\prime}}!}\left(\prod_{i^{\prime}=1}^{n}{q_{i^{\prime}|i}}^{k_{i^{\prime}}}\right)\left(\frac{k_{j}}{2k}+\mathcal{O}(\delta)\right)$
$\displaystyle=\sum_{\sum_{i^{\prime}=1}^{n}k_{i^{\prime}}=k}\left(-\frac{2k_{l}}{kN}\right)x_{i}\frac{k!}{\prod_{i^{\prime}=1}^{n}k_{i^{\prime}}!}\left(\prod_{i^{\prime}=1}^{n}{q_{i^{\prime}|i}}^{k_{i^{\prime}}}\right)\left(\frac{k-k_{i}}{2k}+\mathcal{O}(\delta)\right)$
$\displaystyle=\begin{cases}\displaystyle{x_{i}\frac{q_{i|i}}{N}-x_{i}\frac{q_{i|i}[(k-1)q_{i|i}+1]}{kN}+\mathcal{O}(\delta)},&\mbox{$i=l$,}\\\
\displaystyle{x_{i}\frac{q_{l|i}}{N}-x_{i}\frac{(k-1)q_{l|j}q_{i|j}}{kN}+\mathcal{O}(\delta)},&\mbox{$i\neq
l$.}\end{cases}$
$\displaystyle=\frac{x_{i}}{kN}\begin{cases}\displaystyle{kq_{i|i}-q_{i|i}[(k-1)q_{i|i}+1]+\mathcal{O}(\delta)},&\mbox{$i=l$,}\\\
\displaystyle{kq_{l|i}-(k-1)q_{l|i}q_{i|i}+\mathcal{O}(\delta)},&\mbox{$i\neq
l$.}\end{cases}$ (S151)

In Eq. (B.1.2), $\mathcal{P}(i\leftarrow j)$ refers to Eq. (S10). Again, we
only take the first-order Taylor expansion, because the first-order term would
be non-zero as we will see in Eq. (B.1.3), thus not necessary to study the
second-order term.

####  B.1.3 Separation of different time scales

We combine the increase and decrease of $il$-edges given by Eqs. (B.1.1) and
(B.1.2), and consider that a full MC step contains $N$ elementary steps. In
this way, we have the master equation of the change in the $il$-edge
frequency:

$\displaystyle\dot{p}_{il}$
$\displaystyle=N\times\left(\sum_{\sum_{i^{\prime}=1}^{n}k_{i^{\prime}}=k}\frac{2k_{l}}{kN}\mathcal{P}\left(\Delta
p_{il}=\frac{2k_{l}}{kN}\right)+\sum_{\sum_{i^{\prime}=1}^{n}k_{i^{\prime}}=k}\left(-\frac{2k_{l}}{kN}\right)\mathcal{P}\left(\Delta
p_{il}=-\frac{2k_{l}}{kN}\right)\right)$
$\displaystyle=\frac{x_{i}}{k}\begin{cases}\displaystyle{\sum_{j=1,j\neq
i}^{n}q_{j|i}[(k-1)q_{i|j}+1]-kq_{i|i}+q_{i|i}[(k-1)q_{i|i}+1]+\mathcal{O}(\delta)},&\mbox{$i=l$,}\\\
\displaystyle{\sum_{j=1,j\neq
i}^{n}(k-1)q_{l|j}q_{j|i}-kq_{l|i}+(k-1)q_{l|i}q_{i|i}+\mathcal{O}(\delta)},&\mbox{$i\neq
l$.}\end{cases}$
$\displaystyle=\frac{x_{i}}{k}\begin{cases}\displaystyle{\sum_{j=1}^{n}q_{j|i}+(k-1)\sum_{j=1}^{n}q_{l|j}q_{j|i}-kq_{l|i}+\mathcal{O}(\delta)},&\mbox{$i=l$,}\\\
\displaystyle{(k-1)\sum_{j=1}^{n}q_{l|j}q_{j|i}-kq_{l|i}+\mathcal{O}(\delta)},&\mbox{$i\neq
l$.}\end{cases}$
$\displaystyle=\frac{x_{i}}{k}\left(\theta_{il}+(k-1)\sum_{j=1}^{n}q_{l|j}q_{j|i}-kq_{l|i}\right)+\mathcal{O}(\delta),$
(S152)

where

$\theta_{il}=\begin{cases}\displaystyle{1},&\mbox{$i=l$,}\\\
\displaystyle{0},&\mbox{$i\neq l$.}\end{cases}$ (S153)

Previously, we conclude from Eq. (S1.1) in section S1.1 that the system can be
expressed by $x_{i}$ and $q_{j|i}$, eliminating the need to use $p_{ij}$.
Therefore, we consider $q_{l|i}=p_{il}/x_{i}$ and write $\dot{q}_{l|i}$ by
$\dot{p}_{il}$ as follows.

$\dot{q}_{l|i}=\frac{\mathrm{d}}{\mathrm{d}t}\left(\frac{p_{il}}{x_{i}}\right)=\frac{\dot{p}_{il}x_{i}-\dot{x}_{i}p_{il}}{{x_{i}}^{2}}=\frac{\dot{p}_{il}}{x_{i}}=\frac{1}{k}\left(\theta_{il}+(k-1)\sum_{j=1}^{n}q_{l|j}q_{j|i}-kq_{l|i}\right)+\mathcal{O}(\delta).$
(S154)

In Eq. (S154), we only take the first-order Taylor expansion, as the first-
order term is non-zero and one can ignore the remaining higher order
infinitesimals. It is noted that we have taken $\dot{x}_{i}=0$ in the
calculation, because as seen in Eq. (S2.3), $\dot{x}_{i}$ is zero in its
$\delta^{0}$ term and only becomes non-zero since the $\delta^{1}$ term, which
is the ignored higher order infinitesimals when calculating $\dot{q}_{l|i}$.

Here, an important characteristic arises: the time scales of the changes in
$\dot{x}_{i}$ and $\dot{q}_{l|i}$ are different. Comparing Eqs. (S2.3) and
(S154), we observe that the change in $x_{i}$ happens at the $\delta^{1}$
scale, while the change in $q_{l|i}$ happens at the $\delta^{0}$ scale. In
other words, the change of $q_{l|i}$ is much faster than in $x_{i}$.

Therefore, we can study the dynamics of $x_{i}$ on the basis of $q_{l|i}$
achieving equilibrium (rapidly). That is, we can solve $\dot{q}_{l|i}=0$,
which is

$\theta_{il}+(k-1)\sum_{j=1}^{n}q_{l|j}q_{j|i}-kq_{l|i}=0,$ (S155)

and represent all $q_{l|i}$ quantities by $x_{i}$ quantities. Considering
$x_{i}q_{l|i}=x_{l}q_{i|l}$, the solution of Eq. (S155) is

$q_{l|i}=\frac{(k-2)x_{l}+\theta_{il}}{k-1},$ (S156)

which is Eq. (S13) presented in section S2.3, with the consideration of
$\theta_{il}=1$ if $i=l$ and $\theta_{il}=0$ otherwise.

###  B.2 Death-birth

Similarly, we supplement edge dynamics under pairwise comparison to explain
the origin of Eq. (S132) presented in section A.3.

####  B.2.1 The increase of $il$-edges

The number of $il$-edges in the system increases by $k_{l}$ when a $j$-player
is replaced by an $i$-player ($j\neq i$). The frequency of $il$-edges,
$p_{il}$, is increased by $2k_{l}/(kN)$. For a possible neighbor configuration
$\mathbf{k}$ of the $j$-player, we go through all $j\neq i$, obtaining the
probability that the configuration $\mathbf{k}$ is found around the $j$-player
and thus the frequency of $il$-edges increases by $2k_{l}/(kN)$:

$\mathcal{P}\left(\Delta p_{il}=\frac{2k_{l}}{kN}\right)=\sum_{j=1,j\neq
i}^{n}x_{j}\frac{k!}{\prod_{i^{\prime}=1}^{n}k_{i^{\prime}}!}\left(\prod_{i^{\prime}=1}^{n}{q_{i^{\prime}|j}}^{k_{i^{\prime}}}\right)\mathcal{P}(j\leftarrow
i).$ (S157)

Then, we go through all possible configurations $\mathbf{k}$. The expected
change in the frequency of $il$-edges is calculated as:

$\displaystyle\sum_{\sum_{i^{\prime}=1}^{n}k_{i^{\prime}}=k}\frac{2k_{l}}{kN}\mathcal{P}\left(\Delta
p_{il}=\frac{2k_{l}}{kN}\right)$
$\displaystyle=\sum_{\sum_{i^{\prime}=1}^{n}k_{i^{\prime}}=k}\frac{2k_{l}}{kN}\sum_{j=1,j\neq
i}^{n}x_{j}\frac{k!}{\prod_{i^{\prime}=1}^{n}k_{i^{\prime}}!}\left(\prod_{i^{\prime}=1}^{n}{q_{i^{\prime}|j}}^{k_{i^{\prime}}}\right)\mathcal{P}(j\leftarrow
i)$
$\displaystyle=\sum_{\sum_{i^{\prime}=1}^{n}k_{i^{\prime}}=k}\frac{2k_{l}}{kN}\sum_{j=1,j\neq
i}^{n}x_{j}\frac{k!}{\prod_{i^{\prime}=1}^{n}k_{i^{\prime}}!}\left(\prod_{i^{\prime}=1}^{n}{q_{i^{\prime}|j}}^{k_{i^{\prime}}}\right)\left(\frac{k_{i}}{k}+\mathcal{O}(\delta)\right)$
$\displaystyle=\frac{2x_{i}}{kN}\begin{cases}\displaystyle{\sum_{j=1,j\neq
i}^{n}q_{j|i}[(k-1)q_{i|j}+1]+\mathcal{O}(\delta)},&\mbox{$i=l$,}\\\
\displaystyle{\sum_{j=1,j\neq
i}^{n}(k-1)q_{l|j}q_{j|i}+\mathcal{O}(\delta)},&\mbox{$i\neq l$.}\end{cases}$
(S158)

In Eq. (B.2.1), $\mathcal{P}(j\leftarrow i)$ refers to Eq. (S127). We only
take the first-order Taylor expansion, because the first-order term would be
non-zero as we will see in Eq. (B.2.3).

####  B.2.2 The decrease of $il$-edges

The number of $il$-edges in the system decreases by $k_{l}$ when an $i$-player
is replaced by a $j$-player ($j\neq i$). The frequency of $il$-edges $p_{il}$
is decreased by $2k_{l}/(kN)$. For a possible neighbor configuration
$\mathbf{k}$ of the $i$-player, we go through all possible $j\neq i$,
obtaining the probability that the configuration $\mathbf{k}$ is found around
the $i$-player and thus the frequency of $il$-edges decreases by
$2k_{l}/(kN)$:

$\mathcal{P}\left(\Delta
p_{il}=-\frac{2k_{l}}{kN}\right)=x_{i}\frac{k!}{\prod_{i^{\prime}=1}^{n}k_{i^{\prime}}!}\left(\prod_{i^{\prime}=1}^{n}{q_{i^{\prime}|i}}^{k_{i^{\prime}}}\right)[1-\mathcal{P}(i\leftarrow
i)].$ (S159)

Then, we go through all possible configurations $\mathbf{k}$. The expected
change in the frequency of $il$-edges is calculated as:

$\displaystyle\sum_{\sum_{i^{\prime}=1}^{n}k_{i^{\prime}}=k}\left(-\frac{2k_{l}}{kN}\right)\mathcal{P}\left(\Delta
p_{il}=-\frac{2k_{l}}{kN}\right)$
$\displaystyle=\sum_{\sum_{i^{\prime}=1}^{n}k_{i^{\prime}}=k}\left(-\frac{2k_{l}}{kN}\right)x_{i}\frac{k!}{\prod_{i^{\prime}=1}^{n}k_{i^{\prime}}!}\left(\prod_{i^{\prime}=1}^{n}{q_{i^{\prime}|i}}^{k_{i^{\prime}}}\right)[1-\mathcal{P}(i\leftarrow
i)]$
$\displaystyle=\sum_{\sum_{i^{\prime}=1}^{n}k_{i^{\prime}}=k}\left(-\frac{2k_{l}}{kN}\right)x_{i}\frac{k!}{\prod_{i^{\prime}=1}^{n}k_{i^{\prime}}!}\left(\prod_{i^{\prime}=1}^{n}{q_{i^{\prime}|i}}^{k_{i^{\prime}}}\right)\left(1-\frac{k_{i}}{k}+\mathcal{O}(\delta)\right)$
$\displaystyle=\frac{2x_{i}}{kN}\begin{cases}\displaystyle{kq_{i|i}-q_{i|i}[(k-1)q_{i|i}+1]+\mathcal{O}(\delta)},&\mbox{$i=l$,}\\\
\displaystyle{kq_{l|i}-(k-1)q_{l|i}q_{i|i}+\mathcal{O}(\delta)},&\mbox{$i\neq
l$.}\end{cases}$ (S160)

In Eq. (B.2.2), $[1-\mathcal{P}(i\leftarrow i)]$ refers to Eq. (S129). Again,
we only take the first-order Taylor expansion, because the first-order term
would be non-zero as we will see in Eq. (B.2.3).

####  B.2.3 Separation of different time scales

We combine the increase and decrease of $il$-edges given by Eqs. (B.2.1) and
(B.2.2), and consider that a full MC step contains $N$ elementary steps. In
this way, we have the master equation of the change in the $il$-edge frequency
under death-birth:

$\displaystyle\dot{p}_{il}$
$\displaystyle=N\times\left(\sum_{\sum_{i^{\prime}=1}^{n}k_{i^{\prime}}=k}\frac{2k_{l}}{kN}\mathcal{P}\left(\Delta
p_{il}=\frac{2k_{l}}{kN}\right)+\sum_{\sum_{i^{\prime}=1}^{n}k_{i^{\prime}}=k}\left(-\frac{2k_{l}}{kN}\right)\mathcal{P}\left(\Delta
p_{il}=-\frac{2k_{l}}{kN}\right)\right)$
$\displaystyle=\frac{2x_{i}}{k}\left(\theta_{il}+(k-1)\sum_{j=1}^{n}q_{l|j}q_{j|i}-kq_{l|i}\right)+\mathcal{O}(\delta),$
(S161)

where $\theta_{il}=1$ if $i=l$ and $\theta_{il}=0$ otherwise.

According to Eq. (S1.1), we consider $q_{l|i}=p_{il}/x_{i}$ and write
$\dot{q}_{l|i}$ by $\dot{p}_{il}$ as follows.

$\dot{q}_{l|i}=\frac{\mathrm{d}}{\mathrm{d}t}\left(\frac{p_{il}}{x_{i}}\right)=\frac{2}{k}\left(\theta_{il}+(k-1)\sum_{j=1}^{n}q_{l|j}q_{j|i}-kq_{l|i}\right)+\mathcal{O}(\delta).$
(S162)

In Eq. (S162), we only take the first-order Taylor expansion, as the first-
order term is non-zero and one can ignore the remaining higher order
infinitesimals. Comparing Eqs. (A.3) and (S162), we observe that the change in
$x_{i}$ happens at the $\delta^{1}$ scale, while the change in $q_{l|i}$
happens at the $\delta^{0}$ scale. The change of $q_{l|i}$ are much faster
than $x_{i}$.

Therefore, we can study the dynamics of $x_{i}$ on the basis of $q_{l|i}$
achieving equilibrium (rapidly). That is, we can solve $\dot{q}_{l|i}=0$,
which is $\theta_{il}+(k-1)\sum_{j=1}^{n}q_{l|j}q_{j|i}-kq_{l|i}=0$, and
represent all $q_{l|i}$ quantities by $x_{i}$ quantities. Considering
$x_{i}q_{l|i}=x_{l}q_{i|l}$, the solution is

$q_{l|i}=\frac{(k-2)x_{l}+\theta_{il}}{k-1},$ (S163)

which is Eq. (S132) presented in section A.3, with the consideration of
$\theta_{il}=1$ if $i=l$ and $\theta_{il}=0$ otherwise.

##  Appendix C Theorem on computation and operators

A summation relation central to the theoretical derivation in this Article is
given by Theorem 1.

######  Theorem 1.

For any real number $z_{j}$, where $j=1,2,\dots,n$, satisfying $0\leq
z_{j}\leq 1$ and $\sum_{j=1}^{n}z_{j}=1$, and any scalar multivariate function
$g(\mathbf{k})$ of $\mathbf{k}$, we have the following relation:

$\sum_{\sum_{j=1}^{n}k_{j}=k}\frac{k!}{\prod_{j=1}^{n}k_{j}!}\left(\prod_{j=1}^{n}{z_{j}}^{k_{j}}\right)k_{i}g(\mathbf{k})=kz_{i}\sum_{\sum_{j=1}^{n}k_{j}=k-1}\frac{(k-1)!}{\prod_{j=1}^{n}k_{j}!}\left(\prod_{j=1}^{n}{z_{j}}^{k_{j}}\right)g(\mathbf{k}_{+i}).$
(S164)

###### Proof.

Recall the definitions given in the main text or section S1:
$\mathbf{k}=(k_{1},k_{2},\dots,k_{n})$, where $\sum_{j=1}^{n}k_{j}=k$; and
$\mathbf{k}_{+i}=(k_{1},k_{2},\dots,k_{i}+1,\dots,k_{n})$, where
$\sum_{j=1}^{n}k_{j}=k-1$. In this way, we have

$\displaystyle\sum_{\sum_{j=1}^{n}k_{j}=k}\frac{k!}{\prod_{j=1}^{n}k_{j}!}\left(\prod_{j=1}^{n}{z_{j}}^{k_{j}}\right)k_{i}g(\mathbf{k})$
$\displaystyle=kz_{i}\sum_{\sum_{j=1}^{n}k_{j}=k,k_{i}\neq
0}\frac{(k-1)!}{(k_{i}-1)!\prod_{j=1,j\neq
i}^{n}k_{j}!}\left({z_{i}}^{k_{i}-1}\prod_{j=1,j\neq
i}^{n}{z_{j}}^{k_{j}}\right)g(\mathbf{k})$
$\displaystyle=kz_{i}\sum_{\sum_{j=1}^{n}k_{j}=k-1}\frac{(k-1)!}{\prod_{j=1}^{n}k_{j}!}\left(\prod_{j=1}^{n}{z_{j}}^{k_{j}}\right)g(\mathbf{k}_{+i}).$
(S165)

The first step of Eq. (1) uses $k!=k(k-1)!$,
${z_{i}}^{k_{i}}=z_{i}{z_{i}}^{k_{i}-1}$, $k_{i}/(k_{i}!)=1/[(k_{i}-1)!]$ and
considers that the entire formula equals to zero when $k_{i}=0$. The second
step of Eq. (1) is a straightforward application of the $\mathbf{k}_{+i}$
concept. ∎

Theorem 1 gives a general relation which helps us eliminate $k_{i}$ inside
summations at the cost of transforming $g(\mathbf{k})$ into
$g(\mathbf{k}_{+i})$. For convenience of our readers, there are several useful
applications of Theorem 1.

For example, taking $g(\mathbf{k})=1$, we have

$\sum_{\sum_{j=1}^{n}k_{j}=k}\frac{k!}{\prod_{j=1}^{n}k_{j}!}\left(\prod_{j=1}^{n}{z_{j}}^{k_{j}}\right)k_{i}=kz_{i}.$
(S166)

Another example is $g(\mathbf{k})=k_{l}$. If $l=i$, we have

$\sum_{\sum_{j=1}^{n}k_{j}=k}\frac{k!}{\prod_{j=1}^{n}k_{j}!}\left(\prod_{j=1}^{n}{z_{j}}^{k_{j}}\right)k_{i}k_{l}=kz_{i}\sum_{\sum_{j=1}^{n}k_{j}=k-1}\frac{(k-1)!}{\prod_{j=1}^{n}k_{j}!}\left(\prod_{j=1}^{n}{z_{j}}^{k_{j}}\right)(k_{l}+1)=kz_{i}[(k-1)z_{l}+1],$
(S167)

and if $l\neq i$, it becomes

$\sum_{\sum_{j=1}^{n}k_{j}=k}\frac{k!}{\prod_{j=1}^{n}k_{j}!}\left(\prod_{j=1}^{n}{z_{j}}^{k_{j}}\right)k_{i}k_{l}=kz_{i}\sum_{\sum_{j=1}^{n}k_{j}=k-1}\frac{(k-1)!}{\prod_{j=1}^{n}k_{j}!}\left(\prod_{j=1}^{n}{z_{j}}^{k_{j}}\right)k_{l}=k(k-1)z_{i}z_{l}.$
(S168)

Finally, a potentially misleading example is $g(\mathbf{k})=k^{\prime}_{l}$,
where $k^{\prime}_{l}\in\mathbf{k}^{\prime}\neq\mathbf{k}$ (recall that
$\mathbf{k}^{\prime}$ is an independent variable of $\mathbf{k}$). In this
way, elements in $\mathbf{k}^{\prime}$ are essentially constants when studying
$g(\mathbf{k})$,

$\sum_{\sum_{j=1}^{n}k_{j}=k}\frac{k!}{\prod_{j=1}^{n}k_{j}!}\left(\prod_{j=1}^{n}{z_{j}}^{k_{j}}\right)k_{i}k^{\prime}_{l}=kz_{i}k^{\prime}_{l}.$
(S169)

In this Article, we have created many operators to convey various concepts.
The main relations between these operators are summarized in Theorems 2–5.

######  Theorem 2.

The following equation holds:

$\langle\pi_{i|j}^{\mathbf{k}_{+i}}\rangle=\langle\pi_{i}^{\mathbf{k}_{+j}}\rangle.$
(S170)

Intuitively speaking, these two concepts are equivalent: (1) The statistical
mean payoff of an $i$-player neighboring a $j$-player, where the $j$-player
has at least one $i$-player neighbor. (2) The statistical mean payoff of an
$i$-player, where the $i$-player has at least one $j$-player neighbor.

###### Proof.

Let us expand $\pi_{i|j}^{\mathbf{k}_{+i}}$ according to Eq. (S5),

$\pi_{i|j}^{\mathbf{k}_{+i}}=a_{i|\mathbf{k}_{+j}}+\sum_{\sum_{l=1}^{n}k^{\prime}_{l}=k-1}\frac{(k-1)!}{\prod_{l=1}^{n}k^{\prime}_{l}!}\left(\prod_{l=1}^{n}{q_{l|i}}^{k^{\prime}_{l}}\right)\left(a_{i|\mathbf{k}^{\prime}_{+j}}+\sum_{l=1}^{n}k^{\prime}_{l}\sum_{\sum_{\ell=1}^{n}k^{\prime\prime}_{\ell}=k-1}\frac{(k-1)!}{\prod_{\ell=1}^{n}k^{\prime\prime}_{\ell}!}\left(\prod_{\ell=1}^{n}{q_{\ell|l}}^{k^{\prime\prime}_{\ell}}\right)a_{i|\mathbf{k}^{\prime\prime}_{+l}}\right).$
(S171)

By substituting Eq. (S171) into Eq. (S16), we calculate

$\displaystyle\langle\pi_{i|j}^{\mathbf{k}_{+i}}\rangle=$
$\displaystyle\leavevmode\nobreak\
\sum_{\sum_{i^{\prime}=1}^{n}k_{i^{\prime}}=k-1}\frac{(k-1)!}{\prod_{i^{\prime}=1}^{n}k_{i^{\prime}}!}\left(\prod_{i^{\prime}=1}^{n}{q_{{i^{\prime}}|j}}^{k_{i^{\prime}}}\right)a_{i|\mathbf{k}_{+j}}$
$\displaystyle\leavevmode\nobreak\
+\sum_{\sum_{i^{\prime}=1}^{n}k_{i^{\prime}}=k-1}\frac{(k-1)!}{\prod_{i^{\prime}=1}^{n}k_{i^{\prime}}!}\left(\prod_{i^{\prime}=1}^{n}{q_{{i^{\prime}}|j}}^{k_{i^{\prime}}}\right)\sum_{\sum_{l=1}^{n}k^{\prime}_{l}=k-1}\frac{(k-1)!}{\prod_{l=1}^{n}k^{\prime}_{l}!}\left(\prod_{l=1}^{n}{q_{l|i}}^{k^{\prime}_{l}}\right)\Bigg{[}a_{i|\mathbf{k}^{\prime}_{+j}}$
$\displaystyle\leavevmode\nobreak\
+\sum_{l=1}^{n}k^{\prime}_{l}\sum_{\sum_{\ell=1}^{n}k^{\prime\prime}_{\ell}=k-1}\frac{(k-1)!}{\prod_{\ell=1}^{n}k^{\prime\prime}_{\ell}!}\left(\prod_{\ell=1}^{n}{q_{\ell|l}}^{k^{\prime\prime}_{\ell}}\right)a_{i|\mathbf{k}^{\prime\prime}_{+l}}\Bigg{]}$
$\displaystyle=$ $\displaystyle\leavevmode\nobreak\
\sum_{\sum_{i^{\prime}=1}^{n}k_{i^{\prime}}=k-1}\frac{(k-1)!}{\prod_{i^{\prime}=1}^{n}k_{i^{\prime}}!}\left(\prod_{i^{\prime}=1}^{n}{q_{{i^{\prime}}|j}}^{k_{i^{\prime}}}\right)a_{i|\mathbf{k}_{+j}}+\sum_{\sum_{l=1}^{n}k^{\prime}_{l}=k-1}\frac{(k-1)!}{\prod_{l=1}^{n}k^{\prime}_{l}!}\left(\prod_{l=1}^{n}{q_{l|i}}^{k^{\prime}_{l}}\right)a_{i|\mathbf{k}^{\prime}_{+j}}$
$\displaystyle\leavevmode\nobreak\
+\sum_{\sum_{l=1}^{n}k^{\prime}_{l}=k-1}\frac{(k-1)!}{\prod_{l=1}^{n}k^{\prime}_{l}!}\left(\prod_{l=1}^{n}{q_{l|i}}^{k^{\prime}_{l}}\right)\sum_{l=1}^{n}k^{\prime}_{l}\sum_{\sum_{\ell=1}^{n}k^{\prime\prime}_{\ell}=k-1}\frac{(k-1)!}{\prod_{\ell=1}^{n}k^{\prime\prime}_{\ell}!}\left(\prod_{\ell=1}^{n}{q_{\ell|l}}^{k^{\prime\prime}_{\ell}}\right)a_{i|\mathbf{k}^{\prime\prime}_{+l}}.$
(S172)

Similarly, we expand $\pi_{i}^{\mathbf{k}_{+j}}$ according to Eq. (S3),

$\pi_{i}^{\mathbf{k}_{+j}}=a_{i|\mathbf{k}_{+j}}+\sum_{l=1}^{n}k_{l}\sum_{\sum_{\ell=1}^{n}k^{\prime}_{\ell}=k-1}\frac{(k-1)!}{\prod_{\ell=1}^{n}k^{\prime}_{\ell}!}\left(\prod_{\ell=1}^{n}{q_{\ell|l}}^{k^{\prime}_{\ell}}\right)a_{i|\mathbf{k}^{\prime}_{+l}}.$
(S173)

By substituting Eq. (S173) into Eq. (S17), we calculate

$\displaystyle\langle\pi_{i}^{\mathbf{k}_{+j}}\rangle=$
$\displaystyle\leavevmode\nobreak\
\sum_{\sum_{i^{\prime}=1}^{n}k_{i^{\prime}}=k-1}\frac{(k-1)!}{\prod_{i^{\prime}=1}^{n}k_{i^{\prime}}!}\left(\prod_{i^{\prime}=1}^{n}{q_{{i^{\prime}}|i}}^{k_{i^{\prime}}}\right)\pi_{i}^{\mathbf{k}_{+j}}$
$\displaystyle=$ $\displaystyle\leavevmode\nobreak\
\sum_{\sum_{i^{\prime}=1}^{n}k_{i^{\prime}}=k-1}\frac{(k-1)!}{\prod_{i^{\prime}=1}^{n}k_{i^{\prime}}!}\left(\prod_{i^{\prime}=1}^{n}{q_{{i^{\prime}}|i}}^{k_{i^{\prime}}}\right)a_{i|\mathbf{k}_{+j}}+\sum_{\sum_{\ell=1}^{n}k^{\prime}_{\ell}=k-1}\frac{(k-1)!}{\prod_{\ell=1}^{n}k^{\prime}_{\ell}!}\left(\prod_{\ell=1}^{n}{q_{\ell|j}}^{k^{\prime}_{\ell}}\right)a_{i|\mathbf{k}^{\prime}_{+j}}$
$\displaystyle\leavevmode\nobreak\
+\sum_{\sum_{i^{\prime}=1}^{n}k_{i^{\prime}}=k-1}\frac{(k-1)!}{\prod_{i^{\prime}=1}^{n}k_{i^{\prime}}!}\left(\prod_{i^{\prime}=1}^{n}{q_{{i^{\prime}}|i}}^{k_{i^{\prime}}}\right)\sum_{l=1}^{n}k_{l}\sum_{\sum_{\ell=1}^{n}k^{\prime}_{\ell}=k-1}\frac{(k-1)!}{\prod_{\ell=1}^{n}k^{\prime}_{\ell}!}\left(\prod_{\ell=1}^{n}{q_{\ell|l}}^{k^{\prime}_{\ell}}\right)a_{i|\mathbf{k}^{\prime}_{+l}}.$
(S174)

Comparing Eqs. (2) and (2) and realizing the equivalence of auxiliary
variables, we have
$\langle\pi_{i|j}^{\mathbf{k}_{+i}}\rangle=\langle\pi_{i}^{\mathbf{k}_{+j}}\rangle$.
∎

######  Theorem 3.

The following equation holds:

$\sum_{j=1}^{n}q_{j|i}\langle\pi_{i|j}^{\mathbf{k}_{+i}}\rangle=\langle\pi_{i}^{\mathbf{k}}\rangle.$
(S175)

Intuitive interpretation: (1) the statistical mean payoff of an $i$-player
over all possible neighboring $j$-players found near an $i$-player, where the
$j$-player has at least one $i$-player neighbor, is equivalent to (2) the
statistical mean payoff of an $i$-player.

###### Proof.

Further developing Eq. (2) results in

$\displaystyle\sum_{j=1}^{n}q_{j|i}\langle\pi_{i|j}^{\mathbf{k}_{+i}}\rangle=$
$\displaystyle\sum_{j=1}^{n}q_{j|i}\sum_{\sum_{i^{\prime}=1}^{n}k_{i^{\prime}}=k-1}\frac{(k-1)!}{\prod_{i^{\prime}=1}^{n}k_{i^{\prime}}!}\left(\prod_{i^{\prime}=1}^{n}{q_{{i^{\prime}}|j}}^{k_{i^{\prime}}}\right)a_{i|\mathbf{k}_{+j}}+\sum_{j=1}^{n}q_{j|i}\sum_{\sum_{l=1}^{n}k^{\prime}_{l}=k-1}\frac{(k-1)!}{\prod_{l=1}^{n}k^{\prime}_{l}!}\left(\prod_{l=1}^{n}{q_{l|i}}^{k^{\prime}_{l}}\right)a_{i|\mathbf{k}^{\prime}_{+j}}$
$\displaystyle+\sum_{j=1}^{n}q_{j|i}\sum_{\sum_{l=1}^{n}k^{\prime}_{l}=k-1}\frac{(k-1)!}{\prod_{l=1}^{n}k^{\prime}_{l}!}\left(\prod_{l=1}^{n}{q_{l|i}}^{k^{\prime}_{l}}\right)\sum_{l=1}^{n}k^{\prime}_{l}\sum_{\sum_{\ell=1}^{n}k^{\prime\prime}_{\ell}=k-1}\frac{(k-1)!}{\prod_{\ell=1}^{n}k^{\prime\prime}_{\ell}!}\left(\prod_{\ell=1}^{n}{q_{\ell|l}}^{k^{\prime\prime}_{\ell}}\right)a_{i|\mathbf{k}^{\prime\prime}_{+l}}$
$\displaystyle=$
$\displaystyle\sum_{j=1}^{n}q_{j|i}\sum_{\sum_{i^{\prime}=1}^{n}k_{i^{\prime}}=k-1}\frac{(k-1)!}{\prod_{i^{\prime}=1}^{n}k_{i^{\prime}}!}\left(\prod_{i^{\prime}=1}^{n}{q_{{i^{\prime}}|j}}^{k_{i^{\prime}}}\right)a_{i|\mathbf{k}_{+j}}+\sum_{j=1}^{n}\frac{k_{j}}{k}\sum_{\sum_{l=1}^{n}k^{\prime}_{l}=k}\frac{k!}{\prod_{l=1}^{n}k^{\prime}_{l}!}\left(\prod_{l=1}^{n}{q_{l|i}}^{k^{\prime}_{l}}\right)a_{i|\mathbf{k}^{\prime}}$
$\displaystyle+(k-1)\sum_{\sum_{l=1}^{n}k^{\prime}_{l}=k-2}\frac{(k-2)!}{\prod_{l=1}^{n}k^{\prime}_{l}!}\left(\prod_{l=1}^{n}{q_{l|i}}^{k^{\prime}_{l}}\right)\sum_{l=1}^{n}q_{l|i}\sum_{\sum_{\ell=1}^{n}k^{\prime\prime}_{\ell}=k-1}\frac{(k-1)!}{\prod_{\ell=1}^{n}k^{\prime\prime}_{\ell}!}\left(\prod_{\ell=1}^{n}{q_{\ell|l}}^{k^{\prime\prime}_{\ell}}\right)a_{i|\mathbf{k}^{\prime\prime}_{+l}}$
$\displaystyle=$
$\displaystyle\sum_{j=1}^{n}q_{j|i}\sum_{\sum_{i^{\prime}=1}^{n}k_{i^{\prime}}=k-1}\frac{(k-1)!}{\prod_{i^{\prime}=1}^{n}k_{i^{\prime}}!}\left(\prod_{i^{\prime}=1}^{n}{q_{{i^{\prime}}|j}}^{k_{i^{\prime}}}\right)a_{i|\mathbf{k}_{+j}}+\sum_{\sum_{l=1}^{n}k^{\prime}_{l}=k}\frac{k!}{\prod_{l=1}^{n}k^{\prime}_{l}!}\left(\prod_{l=1}^{n}{q_{l|i}}^{k^{\prime}_{l}}\right)a_{i|\mathbf{k}^{\prime}}$
$\displaystyle+(k-1)\sum_{l=1}^{n}q_{l|i}\sum_{\sum_{\ell=1}^{n}k^{\prime\prime}_{\ell}=k-1}\frac{(k-1)!}{\prod_{\ell=1}^{n}k^{\prime\prime}_{\ell}!}\left(\prod_{\ell=1}^{n}{q_{\ell|l}}^{k^{\prime\prime}_{\ell}}\right)a_{i|\mathbf{k}^{\prime\prime}_{+l}}$
$\displaystyle=$
$\displaystyle\sum_{\sum_{l=1}^{n}k^{\prime}_{l}=k}\frac{k!}{\prod_{l=1}^{n}k^{\prime}_{l}!}\left(\prod_{l=1}^{n}{q_{l|i}}^{k^{\prime}_{l}}\right)a_{i|\mathbf{k}^{\prime}}+k\sum_{l=1}^{n}q_{l|i}\sum_{\sum_{\ell=1}^{n}k^{\prime\prime}_{\ell}=k-1}\frac{(k-1)!}{\prod_{\ell=1}^{n}k^{\prime\prime}_{\ell}!}\left(\prod_{\ell=1}^{n}{q_{\ell|l}}^{k^{\prime\prime}_{\ell}}\right)a_{i|\mathbf{k}^{\prime\prime}_{+l}}.$
(S176)

On the other hand, we apply the expression of $\pi_{i}^{\mathbf{k}}$ in Eq.
(S3) to Eq. (S17) and calculate

$\displaystyle\langle\pi_{i}^{\mathbf{k}}\rangle=$
$\displaystyle\leavevmode\nobreak\
\sum_{\sum_{i^{\prime}=1}^{n}k_{i^{\prime}}=k}\frac{k!}{\prod_{i^{\prime}=1}^{n}k_{i^{\prime}}!}\left(\prod_{i^{\prime}=1}^{n}{q_{{i^{\prime}}|i}}^{k_{i^{\prime}}}\right)\pi_{i}^{\mathbf{k}}$
$\displaystyle=$ $\displaystyle\leavevmode\nobreak\
\sum_{\sum_{i^{\prime}=1}^{n}k_{i^{\prime}}=k}\frac{k!}{\prod_{i^{\prime}=1}^{n}k_{i^{\prime}}!}\left(\prod_{i^{\prime}=1}^{n}{q_{{i^{\prime}}|i}}^{k_{i^{\prime}}}\right)a_{i|\mathbf{k}}$
$\displaystyle\leavevmode\nobreak\
+\sum_{\sum_{i^{\prime}=1}^{n}k_{i^{\prime}}=k}\frac{k!}{\prod_{i^{\prime}=1}^{n}k_{i^{\prime}}!}\left(\prod_{i^{\prime}=1}^{n}{q_{{i^{\prime}}|i}}^{k_{i^{\prime}}}\right)\sum_{l=1}^{n}k_{l}\sum_{\sum_{\ell=1}^{n}k^{\prime}_{\ell}=k-1}\frac{(k-1)!}{\prod_{\ell=1}^{n}k^{\prime}_{\ell}!}\left(\prod_{\ell=1}^{n}{q_{\ell|l}}^{k^{\prime}_{\ell}}\right)a_{i|\mathbf{k}^{\prime}_{+l}}$
$\displaystyle=$ $\displaystyle\leavevmode\nobreak\
\sum_{\sum_{i^{\prime}=1}^{n}k_{i^{\prime}}=k}\frac{k!}{\prod_{i^{\prime}=1}^{n}k_{i^{\prime}}!}\left(\prod_{i^{\prime}=1}^{n}{q_{{i^{\prime}}|i}}^{k_{i^{\prime}}}\right)a_{i|\mathbf{k}}+k\sum_{l=1}^{n}q_{l|i}\sum_{\sum_{\ell=1}^{n}k^{\prime}_{\ell}=k-1}\frac{(k-1)!}{\prod_{\ell=1}^{n}k^{\prime}_{\ell}!}\left(\prod_{\ell=1}^{n}{q_{\ell|l}}^{k^{\prime}_{\ell}}\right)a_{i|\mathbf{k}^{\prime}_{+l}}.$
(S177)

Comparing Eqs. (3) and (3) and noting that $\mathbf{k}^{\prime}$,
$\mathbf{k}^{\prime\prime}_{+l}$, $\mathbf{k}$, and $\mathbf{k}^{\prime}_{+l}$
are calculated separately here (i.e., we can treat $\mathbf{k}^{\prime}$,
$\mathbf{k}$ as $\mathbf{k}$ and treat $\mathbf{k}^{\prime\prime}_{+l}$,
$\mathbf{k}^{\prime}_{+l}$ as $\mathbf{k}_{+l}$), we have
$\sum_{j=1}^{n}q_{j|i}\langle\pi_{i|j}^{\mathbf{k}_{+i}}\rangle=\langle\pi_{i}^{\mathbf{k}}\rangle$.
∎

######  Theorem 4.

The following equation holds:

$\langle a_{i|\mathbf{k}}\rangle_{i}=\sum_{j=1}^{n}q_{j|i}\langle
a_{i|\mathbf{k}_{+j}}\rangle_{i}.$ (S178)

Intuitive interpretation: (1) the expected payoff of an $i$-player in a single
game is equivalent to (2) the expected payoff of an $i$-player in a single
game conditioning on all possible neighboring $j$-players, where the
$i$-player has at least one $j$-player neighbor.

###### Proof.

According to the definition of $\langle a_{i|\mathbf{k}}\rangle_{j}$ in Eq.
(S21), we have

$\displaystyle\langle a_{i|\mathbf{k}}\rangle_{i}$
$\displaystyle=\sum_{\sum_{i^{\prime}=1}^{n}k_{i^{\prime}}=k}\frac{k!}{\prod_{i^{\prime}=1}^{n}k_{i^{\prime}}!}\left(\prod_{i^{\prime}=1}^{n}{q_{{i^{\prime}}|i}}^{k_{i^{\prime}}}\right)a_{i|\mathbf{k}}$
$\displaystyle=\sum_{j=1}^{n}\sum_{\sum_{i^{\prime}=1}^{n}k_{i^{\prime}}=k}\frac{k!}{\prod_{i^{\prime}=1}^{n}k_{i^{\prime}}!}\left(\prod_{i^{\prime}=1}^{n}{q_{{i^{\prime}}|i}}^{k_{i^{\prime}}}\right)\frac{k_{j}}{k}a_{i|\mathbf{k}}$
$\displaystyle=\sum_{j=1}^{n}q_{j|i}\sum_{\sum_{i^{\prime}=1}^{n}k_{i^{\prime}}=k-1}\frac{(k-1)!}{\prod_{i^{\prime}=1}^{n}k_{i^{\prime}}!}\left(\prod_{i^{\prime}=1}^{n}{q_{{i^{\prime}}|i}}^{k_{i^{\prime}}}\right)a_{i|\mathbf{k}_{+j}}$
$\displaystyle=\sum_{j=1}^{n}q_{j|i}\langle a_{i|\mathbf{k}_{+j}}\rangle_{i},$
(S179)

which completes the proof. ∎

######  Theorem 5.

The following equation holds:

$\sum_{j=1}^{n}\sum_{i^{\prime}=1}^{n}q_{i^{\prime}|j}q_{j|i}\langle\pi_{i|j}^{\mathbf{k}_{+i,+i^{\prime}}}\rangle=\langle\pi_{i}^{\mathbf{k}}\rangle.$
(S180)

Intuitive interpretation: (1) the statistical mean payoff of an $i$-player
over all possible neighboring $j$-players found near an $i$-player and all
possible neighboring $i^{\prime}$-players found near an $j$-player, where the
$j$-player has at least one $i$-player neighbor and one $i^{\prime}$-player
neighbor, is equivalent to (2) the statistical mean payoff of an $i$-player.

###### Proof.

According to Theorem 3, we have
$\langle\pi_{i}^{\mathbf{k}}\rangle=\sum_{j=1}^{n}q_{j|i}\langle\pi_{i|j}^{\mathbf{k}_{+i}}\rangle$.
Furthermore, we calculate

$\displaystyle\langle\pi_{i|j}^{\mathbf{k}_{+i}}\rangle=$
$\displaystyle\sum_{\sum_{i^{\prime}=1}^{n}k_{i^{\prime}}=k-1}\frac{(k-1)!}{\prod_{i^{\prime}=1}^{n}k_{i^{\prime}}!}\left(\prod_{i^{\prime}=1}^{n}{q_{{i^{\prime}}|j}}^{k_{i^{\prime}}}\right)\pi_{i|j}^{\mathbf{k}_{+i}}$
$\displaystyle=$
$\displaystyle\sum_{\sum_{i^{\prime}=1}^{n}k_{i^{\prime}}=k-1}\frac{(k-1)!}{\prod_{i^{\prime}=1}^{n}k_{i^{\prime}}!}\left(\prod_{i^{\prime}=1}^{n}{q_{{i^{\prime}}|j}}^{k_{i^{\prime}}}\right)\sum_{i^{\prime}=1}^{n}\frac{k_{i^{\prime}}}{k-1}\pi_{i|j}^{\mathbf{k}_{+i}}$
$\displaystyle=$
$\displaystyle\sum_{\sum_{i^{\prime}=1}^{n}k_{i^{\prime}}=k-2}\frac{(k-2)!}{\prod_{i^{\prime}=1}^{n}k_{i^{\prime}}!}\left(\prod_{i^{\prime}=1}^{n}{q_{{i^{\prime}}|j}}^{k_{i^{\prime}}}\right)\sum_{i^{\prime}=1}^{n}q_{i^{\prime}|j}\pi_{i|j}^{\mathbf{k}_{+i,+i^{\prime}}}$
$\displaystyle=$
$\displaystyle\sum_{i^{\prime}=1}^{n}q_{i^{\prime}|j}\langle\pi_{i|j}^{\mathbf{k}_{+i,+i^{\prime}}}\rangle.$
(S181)

Applying Eq. (5) to Theorem 3, we obtain
$\langle\pi_{i}^{\mathbf{k}}\rangle=\sum_{j=1}^{n}\sum_{i^{\prime}=1}^{n}q_{i^{\prime}|j}q_{j|i}\langle\pi_{i|j}^{\mathbf{k}_{+i,+i^{\prime}}}\rangle$,
which completes the proof. ∎

## References

  * [1] Szolnoki, A. _et al._ Cyclic dominance in evolutionary games: a review.  _Journal of the Royal Society Interface_ 11, 20140735, DOI:10.1098/rsif.2014.0735 (2014).
  * [2] Sinervo, B. & Lively, C. M.  The rock-paper-scissors game and the evolution of alternative male strategies.  _Nature_ 380, 240–243, DOI:10.1038/380240a0 (1996).
  * [3] Kerr, B., Riley, M. A., Feldman, M. W. & Bohannan, B. J.  Local dispersal promotes biodiversity in a real-life game of rock-paper-scissors.  _Nature_ 418, 171–174, DOI:10.1038/nature00823 (2002).
  * [4] Hauert, C., De Monte, S., Hofbauer, J. & Sigmund, K.  Volunteering as red queen mechanism for cooperation in public goods games.  _Science_ 296, 1129–1132, DOI:10.1126/science.1070582 (2002).
  * [5] Hofbauer, J. & Sigmund, K.  _Evolutionary Games and Population Dynamics_ (Cambridge University Press, 1998).
  * [6] Semmann, D., Krambeck, H.-J. & Milinski, M.  Volunteering leads to rock-paper-scissors dynamics in a public goods game.  _Nature_ 425, 390–393, DOI:10.1038/nature01986 (2003).
  * [7] Fehr, E. & Gächter, S.  Altruistic punishment in humans.  _Nature_ 415, 137–140, DOI:10.1038/415137a (2002).
  * [8] Szolnoki, A., Szabó, G. & Perc, M.  Phase diagrams for the spatial public goods game with pool punishment.  _Physical Review E_ 83, 036101, DOI:10.1103/PhysRevE.83.036101 (2011).
  * [9] Sigmund, K., Hauert, C. & Nowak, M. A.  Reward and punishment.  _Proceedings of the National Academy of Sciences_ 98, 10757–10762, DOI:10.1073/pnas.161155698 (2001).
  * [10] Wang, S., Liu, L. & Chen, X.  Tax-based pure punishment and reward in the public goods game.  _Physics Letters A_ 386, 126965, DOI:10.1016/j.physleta.2020.126965 (2021).
  * [11] Sigmund, K.  _The calculus of selfishness_ (Princeton University Press, 2010).
  * [12] Nowak, M. A. & May, R. M.  Evolutionary games and spatial chaos.  _Nature_ 359, 826–829, DOI:10.1038/359826a0 (1992).
  * [13] Nowak, M. A.  Five rules for the evolution of cooperation.  _Science_ 314, 1560–1563, DOI:10.1126/science.1133755 (2006).
  * [14] Nowak, M. A., Tarnita, C. E. & Antal, T.  Evolutionary dynamics in structured populations.  _Philosophical Transactions of the Royal Society B: Biological Sciences_ 365, 19–30, DOI:10.1098/rstb.2009.0215 (2010).
  * [15] Ibsen-Jensen, R., Chatterjee, K. & Nowak, M. A.  Computational complexity of ecological and evolutionary spatial dynamics.  _Proceedings of the National Academy of Sciences_ 112, 15636–15641, DOI:10.1073/pnas.1511366112 (2015).
  * [16] Lieberman, E., Hauert, C. & Nowak, M. A.  Evolutionary dynamics on graphs.  _Nature_ 433, 312–316, DOI:10.1038/nature03204 (2005).
  * [17] Taylor, P. D., Day, T. & Wild, G.  Evolution of cooperation in a finite homogeneous graph.  _Nature_ 447, 469–472, DOI:10.1038/nature05784 (2007).
  * [18] Allen, B. & Nowak, M. A.  Games on graphs.  _EMS Surveys in Mathematical Sciences_ 1, 113–151, DOI:10.4171/EMSS/3 (2014).
  * [19] Débarre, F., Hauert, C. & Doebeli, M.  Social evolution in structured populations.  _Nature Communications_ 5, 3409, DOI:10.1038/ncomms4409 (2014).
  * [20] Allen, B. _et al._ Evolutionary dynamics on any population structure.  _Nature_ 544, 227–230, DOI:10.1038/nature21723 (2017).
  * [21] McAvoy, A., Allen, B. & Nowak, M. A.  Social goods dilemmas in heterogeneous societies.  _Nature Human Behaviour_ 4, 819–831, DOI:10.1038/s41562-020-0881-2 (2020).
  * [22] Su, Q., McAvoy, A., Mori, Y. & Plotkin, J. B.  Evolution of prosocial behaviours in multilayer populations.  _Nature Human Behaviour_ 6, 338–348, DOI:10.1038/s41562-021-01241-2 (2022).
  * [23] Su, Q., McAvoy, A. & Plotkin, J. B.  Strategy evolution on dynamic networks.  _Nature Computational Science_ 3, 763–776, DOI:10.1038/s43588-023-00509-z (2023).
  * [24] Gutowitz, H. A., Victor, J. D. & Knight, B. W.  Local structure theory for cellular automata.  _Physica D_ 28, 18–48, DOI:10.1016/0167-2789(87)90120-5 (1987).
  * [25] Matsuda, H., Tamachi, N., Sasaki, A. & Ogita, N.  A lattice model for population biology.  In _Mathematical Topics in Population Biology, Morphogenesis and Neurosciences_ , vol. 71, 154–161, DOI:10.1007/978-3-642-93360-8_15 (Springer, 1987).
  * [26] Szabó, G., Szolnoki, A. & Bodócs, L.  Correlations induced by transport in one-dimensional lattice gas.  _Physical Review A_ 44, 6375, DOI:10.1103/PhysRevA.44.6375 (1991).
  * [27] Matsuda, H., Ogita, N., Sasaki, A. & Satō, K.  Statistical mechanics of population: the lattice Lotka-Volterra model.  _Progress of Theoretical Physics_ 88, 1035–1049, DOI:10.1143/ptp/88.6.1035 (1992).
  * [28] Szabó, G. & Fath, G.  Evolutionary games on graphs.  _Physics Reports_ 446, 97–216, DOI:10.1016/j.physrep.2007.04.004 (2007).
  * [29] Ohtsuki, H., Hauert, C., Lieberman, E. & Nowak, M. A.  A simple rule for the evolution of cooperation on graphs and social networks.  _Nature_ 441, 502–505, DOI:10.1038/nature04605 (2006).
  * [30] Ohtsuki, H., Nowak, M. A. & Pacheco, J. M.  Breaking the symmetry between interaction and replacement in evolutionary dynamics on graphs.  _Physical Review Letters_ 98, 108106, DOI:10.1103/PhysRevLett.98.108106 (2007).
  * [31] Su, Q., Allen, B. & Plotkin, J. B.  Evolution of cooperation with asymmetric social interactions.  _Proceedings of the National Academy of Sciences_ 119, e2113468118, DOI:10.1073/pnas.2113468118 (2022).
  * [32] Su, Q., McAvoy, A., Wang, L. & Nowak, M. A.  Evolutionary dynamics with game transitions.  _Proceedings of the National Academy of Sciences_ 116, 25398–25404, DOI:10.1073/pnas.1908936116 (2019).
  * [33] Ohtsuki, H. & Nowak, M. A.  The replicator equation on graphs.  _Journal of Theoretical Biology_ 243, 86–97, DOI:10.1016/j.jtbi.2006.06.004 (2006).
  * [34] Taylor, P. D. & Jonker, L. B.  Evolutionary stable strategies and game dynamics.  _Mathematical Biosciences_ 40, 145–156, DOI:10.1016/0025-5564(78)90077-9 (1978).
  * [35] Perc, M., Gómez-Gardenes, J., Szolnoki, A., Floría, L. M. & Moreno, Y.  Evolutionary dynamics of group interactions on structured populations: a review.  _Journal of the Royal Society Interface_ 10, 20120997, DOI:10.1098/rsif.2012.0997 (2013).
  * [36] Alvarez-Rodriguez, U. _et al._ Evolutionary dynamics of higher-order interactions in social networks.  _Nature Human Behaviour_ 5, 586–595, DOI:10.1038/s41562-020-01024-1 (2021).
  * [37] Battiston, F. _et al._ The physics of higher-order interactions in complex systems.  _Nature Physics_ 17, 1093–1098, DOI:10.1038/s41567-021-01371-4 (2021).
  * [38] Tarnita, C. E., Ohtsuki, H., Antal, T., Fu, F. & Nowak, M. A.  Strategy selection in structured populations.  _Journal of Theoretical Biology_ 259, 570–581, DOI:10.1016/j.jtbi.2009.03.035 (2009).
  * [39] Tarnita, C. E., Wage, N. & Nowak, M. A.  Multiple strategies in structured populations.  _Proceedings of the National Academy of Sciences_ 108, 2334–2337, DOI:10.1073/pnas.1016008108 (2011).
  * [40] McAvoy, A. & Wakeley, J.  Evaluating the structure-coefficient theorem of evolutionary game theory.  _Proceedings of the National Academy of Sciences_ 119, e2119656119, DOI:10.1073/pnas.2119656119 (2022).
  * [41] Wu, B., Traulsen, A. & Gokhale, C. S.  Dynamic properties of evolutionary multi-player games in finite populations.  _Games_ 4, 182–199, DOI:10.3390/g4020182 (2013).
  * [42] McAvoy, A. & Hauert, C.  Structure coefficients and strategy selection in multiplayer games.  _Journal of Mathematical Biology_ 72, 203–238, DOI:10.1007/s00285-015-0882-3 (2016).
  * [43] Li, A., Wu, B. & Wang, L.  Cooperation with both synergistic and local interactions can be worse than each alone.  _Scientific Reports_ 4, 5536, DOI:10.1038/srep05536 (2014).
  * [44] Li, A., Broom, M., Du, J. & Wang, L.  Evolutionary dynamics of general group interactions in structured populations.  _Physical Review E_ 93, 022407, DOI:10.1103/PhysRevE.93.022407 (2016).
  * [45] Su, Q., Li, A., Wang, L. & Eugene Stanley, H.  Spatial reciprocity in the evolution of cooperation.  _Proceedings of the Royal Society B_ 286, 20190041, DOI:10.1098/rspb.2019.0041 (2019).
  * [46] Wang, C. & Szolnoki, A.  Inertia in spatial public goods games under weak selection.  _Applied Mathematics and Computation_ 449, 127941, DOI:10.1016/j.amc.2023.127941 (2023).
  * [47] Sigmund, K.  Punish or perish? retaliation and collaboration among humans.  _Trends in ecology & evolution_ 22, 593–600, DOI:10.1016/j.tree.2007.06.012 (2007).
  * [48] Helbing, D., Szolnoki, A., Perc, M. & Szabó, G.  Punish, but not too hard: how costly punishment spreads in the spatial public goods game.  _New Journal of Physics_ 12, 083005, DOI:10.1088/1367-2630/12/8/083005 (2010).
  * [49] Szolnoki, A., Szabó, G. & Czakó, L.  Competition of individual and institutional punishments in spatial public goods games.  _Physical Review E_ 84, 046106, DOI:10.1103/PhysRevE.84.046106 (2011).
  * [50] Rand, D. G., Dreber, A., Ellingsen, T., Fudenberg, D. & Nowak, M. A.  Positive interactions promote public cooperation.  _Science_ 325, 1272–1275, DOI:10.1126/science.1177418 (2009).
  * [51] Hilbe, C. & Sigmund, K.  Incentives and opportunism: from the carrot to the stick.  _Proceedings of the Royal Society B: Biological Sciences_ 277, 2427–2433, DOI:10.1098/rspb.2010.0065 (2010).
  * [52] Szabó, G. & Hauert, C.  Phase transitions and volunteering in spatial public goods games.  _Physical Review Letters_ 89, 118101, DOI:10.1103/PhysRevLett.89.118101 (2002).
  * [53] Wang, S., Chen, X., Xiao, Z. & Szolnoki, A.  Decentralized incentives for general well-being in networked public goods game.  _Applied Mathematics and Computation_ 431, 127308, DOI:10.1016/j.amc.2022.127308 (2022).
  * [54] Sun, Z., Chen, X. & Szolnoki, A.  State-dependent optimal incentive allocation protocols for cooperation in public goods games on regular networks.  _IEEE Transactions on Network Science and Engineering_ 10, 3975–3988, DOI:10.1109/TNSE.2023.3279094 (2023).
  * [55] Szabó, G. & Tőke, C.  Evolutionary prisoner’s dilemma game on a square lattice.  _Physical Review E_ 58, 69, DOI:10.1103/PhysRevE.58.69 (1998).
  * [56] Ohtsuki, H., Iwasa, Y. & Nowak, M. A.  Indirect reciprocity provides only a narrow margin of efficiency for costly punishment.  _Nature_ 457, 79–82, DOI:10.1038/nature07601 (2009).
  * [57] Sigmund, K., De Silva, H., Traulsen, A. & Hauert, C.  Social learning promotes institutions for governing the commons.  _Nature_ 466, 861–863, DOI:10.1038/nature09203 (2010).
  * [58] Hardin, G.  The tragedy of the commons.  _Science_ 162, 1243–1248, DOI:10.1126/science.162.3859.1243 (1968).
  * [59] Wang, C., Zhu, W. & Szolnoki, A.  The conflict between self-interaction and updating passivity in the evolution of cooperation.  _Chaos, Solitons & Fractals_ 173, 113667, DOI:10.1016/j.chaos.2023.113667 (2023).
  * [60] Wang, C., Zhu, W. & Szolnoki, A.  When greediness and self-confidence meet in a social dilemma.  _Physica A_ 625, 129033, DOI:10.1016/j.physa.2023.129033 (2023).
  * [61] Brandt, H., Hauert, C. & Sigmund, K.  Punishment and reputation in spatial public goods games.  _Proceedings of the Royal Society B_ 270, 1099–1104, DOI:10.1098/rspb.2003.2336 (2003).
  * [62] Helbing, D., Szolnoki, A., Perc, M. & Szabó, G.  Evolutionary establishment of moral and double moral standards through spatial interactions.  _PLoS Computational Biology_ 6, e1000758, DOI:10.1371/journal.pcbi.1000758 (2010).
  * [63] Clifford, P. & Sudbury, A.  A model for spatial conflict.  _Biometrika_ 60, 581–588, DOI:10.1093/biomet/60.3.581 (1973).
  * [64] Liggett, T. M.  _Interacting Particle Systems_ (Springer, 1985).
  * [65] Helbing, D., Szolnoki, A., Perc, M. & Szabó, G.  Defector-accelerated cooperativeness and punishment in public goods games with mutations.  _Physical Review E_ 81, 057104, DOI:10.1103/PhysRevE.81.057104 (2010).
  * [66] Sasaki, T., Uchida, S. & Chen, X.  Voluntary rewards mediate the evolution of pool punishment for maintaining public goods in large populations.  _Scientific Reports_ 5, 8917, DOI:10.1038/srep08917 (2015).
  * [67] Luo, Q., Liu, L. & Chen, X.  Evolutionary dynamics of cooperation in the $N$-person stag hunt game.  _Physica D_ 424, 132943, DOI:10.1016/j.physd.2021.132943 (2021).
  * [68] Szolnoki, A. & Chen, X.  Tactical cooperation of defectors in a multi-stage public goods game.  _Chaos, Solitons & Fractals_ 155, 111696, DOI:10.1016/j.chaos.2021.111696 (2022).