Average abundancy of cooperation in multi-player games with random payoffs

We consider interactions between players in groups of size $d\geq 2$ with payoffs that not only depend on the strategies used in the group but also fluctuate at random over time. An individual can adopt either cooperation or defection as strategy and the population is updated from one-time step to the next by a birth-death event according to a Moran model. Assuming recurrent symmetric mutation and payoffs with expected values, variances, and covariances of the same small order, we derive a first-order approximation of the average abundance of cooperation in the selection-mutation equilibrium. We show that increasing the variance of any payoff for defection or decreasing the variance of any payoff for cooperation increases the average abundance of cooperation. As for the effect of the covariance between any payoff for cooperation and any payoff for defection, we show that it depends on the number of cooperators in the group associated with these payoffs. We study in particular the public goods game, the stag hunt game, and the snowdrift game, all social dilemmas based on random benefit $b$ and cost $c$ for cooperation. We show that a decrease in the scaled variance of $b$ or $c$, or an increase in their scaled covariance, makes it easier for weak selection to favor the abundance of cooperation in the stag hunt game and the snowdrift game. The same conclusion holds for the public goods game except that the covariance of $b$ has no effect on the average abundance of $C$. On the other hand, increasing the scaled mutation rate or the group size can enhance or lessen the condition for weak selection to favor the abundance of $C$.

Evolutionary game theory assumes that the reproductive success of a strategy is not constant but depends on the frequencies of the different strategies (Maynard Smith and Price [40], Maynard Smith [39], Hofbauer and Sigmund [19]). This reproductive success is a function of a mean payoff that depends on interactions between individuals, which in turn depend on the population state.

In continuous time, the mean payoff of a strategy has been defined as its growth rate. Then, in an infinitely large well-mixed population, the deterministic dynamics is described by the replicator equation (Taylor and Jonker [57], Zeeman [61]). Evolutionary concepts such as evolutionarily stability (Maynard Smith and Price [40]), continuous stability (Eshel [5]) or convergence stability (Christiansen [3]) were first studied in this framework (see, e.g., Taylor [58], Hofbauer and Sigmund [19]).

For evolutionary game dynamics in a finite population, we have to resort to stochastic processes and probability concepts. Consider, for instance, a well-mixed population of fixed size $N\geq 2$ in discrete time, where each individual can use as strategy either $C$ for cooperation or $D$ for defection in a Prisoner’s dilemma. With any update rule from one time step to the next and in the absence of mutation, the population state over time is described by a Markov chain that has two absorbing states corresponding to an all cooperating population and an all defecting population. Let $\rho_{C}$ (respectively, $\rho_{D}$) be the probability that a single individual of type $C$ (respectively, type $D$) among $N-1$ individuals of type $D$ (respectively, type $C$) generates a lineage forward in time that will take over the whole population. Then, selection is said to favor the evolution of $C$ more than the evolution of $D$ if $\rho_{C}>\rho_{D}$ (Nowak et al. [43], Imhof and Nowak [21]). In the presence of symmetric mutation, the Markov chain is irreducible and, as a result, it possesses a stationary state called the mutation-selection equilibrium. If the average frequency of $C$ in this equilibrium is greater that the average frequency of $D$, then selection is said to favor the abundance of $C$ (Antal et al. [1], see also Kroumi and Lessard [29, 30]).

The above game dynamics were considered first under the assumption of random pairwise interactions. Multi-player games were considered later on to take into account interactions within groups of any fixed size $d\geq 2$. Kurokawa and Ihara [33], for instance, studied the probability of ultimate fixation of a strategy given its initial frequency in this framework in the absence of the mutation, and Gokhale and Traulsen [15] its average abundance in the presence of recurrent mutation. In particular, an exact condition for weak selection to favor the abundance of a given strategy in a large population were deduced in the case $d=3$.

In all the above mentioned studies, the payoffs were supposed deterministic. This assumes a constant environment which is not generally the case in biological populations. Unknown risk of predation and variability in available resources, competition capabilities as well as birth and death rates (May [38], Kaplan et al. [24], Lande et al. [34]) are among multiple factors that motivate taking into account random fluctuations in evolutionary models. Indeed, such fluctuations have an effect on evolutionary outcomes. Early studies on the effect of varying the selection coefficients between generations and varying the offspring numbers within generations in haploid as well as diploid population genetic models in the absence of mutation include Gillespie [13, 14], Karlin and Levikson [25], Karlin and Liberman [26] and Frank and Slatkin [11]). Extensions can be found in Starrfelt and Kokko [55], Schreiber [52] or Rychtar and Taylor [49]. Moreover, the fixation probability for a given type in a population whose size fluctuates dynamically was addressed in Lambert [35], Parsons and Quince [46, 47] and Otto and Whitlock [44], among others, while Uecker and Hermisson [59] studied a population with temporal variation not only in its size but also in selection pressure. Competing populations distributed over habitat patches where environmental conditions fluctuate in time and space were considered too (Evans et al. [6], Schreiber [51]).

In evolutionary games, the payoffs may change at random over time. A stochastic version of the continuous-time replicator equation with a random noise added to the growth rate of every strategy was considered in Fudenberg and Harris [12]. More recently, the effect of stochastic changes in payoffs in two-player linear games in discrete time were studied with particular attention to stochastic local stability of fixation states and polymorphic equilibria in an infinite population (Zheng et al. [63, 64]) and fixation probabilities of strategies in a finite population (Li and Lessard [36], Kroumi et al. [32]).

In a previous paper (Kroumi and Lessard [31]), we have studied the effects of randomness in payoffs on the average abundance of strategies under recurrent mutation in two-player games. Assuming a finite population in discrete time updated according to a Moran model. The payoffs for cooperation and defection in Prisoner’s dilemmas, repeated or not, fluctuate over time such that their means, variances and covariances are of the same small order while higher moments are insignificant. In the mutation-selection equilibrium, we have shown that an increase in the variance of any payoff received by defection against cooperation or defection, or in their covariance, or a decrease in the variance of any payoff received by cooperation against cooperation or defection, or in their covariance, increases the average abundance of cooperation. Then, it is easier for weak selection to favor the abundance of cooperation. Moreover, increasing the scaled mutation rate can lessen or enhance the effect.

In this paper, we extend our analysis of average abundance with random payoffs to many-player games. Interactions occur within groups of any fixed size $d\geq 2$ and the payoffs received by cooperators and defectors are random variables that depend on the group composition. We study the average abundance of $C$ in the stationary state of the mutation-selection equilibrium to derive conditions for weak selection to favor the abundance of $C$. These conditions are examined in detail to understand the effect of the scaled mutation rate and the group size in different games and scenarios.

The remainder of this paper is organized as follows. In Section $2$, we formulate the model. In Section 3, we derive the average abundance of $C$ in the stationary state under symmetric mutation and weak selection. In Section $4$, we examine in detail the particular case of payoffs with additive mean scaled cost and benefit for cooperation. In the next sections, we focus on three special social dilemmas, the public goods game in Section $5$, the stag hunt game in Section $6$, and the snowdrift game in Section $7$. We summarize the conclusions and interpretations of our results in Section $8$.

Consider a population of a fixed size $N$. Each individual can be of only one of two types depending on the strategy used: $C$ for cooperation or $D$ for defection. Interactions between individuals occur within random groups of $d$ players. A cooperator (respectively, defector) receives a payoff $a_{k}$ (respectively, $b_{k}$) when it interacts with $k$ cooperators and $d-k-1$ defectors. The payoffs at a given time step are described in the following payoff array:

Here, we suppose that the payoffs are random variables whose first and second moments are given by

	$\displaystyle E\left[a_{i}\right]=\mu_{C,i}\delta+o(\delta),$		(1a)
	$\displaystyle E\left[a_{i}a_{j}\right]=\sigma_{CC,ij}\delta+o(\delta),$		(1b)
	$\displaystyle E\left[b_{i}\right]=\mu_{D,i}\delta+o(\delta),$		(1c)
	$\displaystyle E\left[b_{i}b_{j}\right]=\sigma_{DD,ij}\delta+o(\delta),$		(1d)
	$\displaystyle E\left[a_{i}b_{j}\right]=\sigma_{CD,ij}\delta+o(\delta),$		(1e)

for $i,j=0,1,\ldots,d-1$. The parameter $\delta\geq 0$, which corresponds to an intensity of selection, measures the order of the first and second moments of the payoffs. The parameters $\mu_{S_{1},i}$ and $\sigma_{S_{1}S_{2},ij}$ for $i,j=0,1,\ldots,d-1$ and $S_{1},S_{2}=C$ or $D$ correspond to scaled means and covariances, respectively. In addition, all higher-order moments of the payoffs are assumed to be negligible compared to $\delta$, that is,

$$E\left[\prod_{i=0}^{d-1}a_{i}^{k_{i}}b_{i}^{l_{i}}\right]=o(\delta)$$

(2)

as soon as $\sum_{i=0}^{d-1}(k_{i}+l_{i})\geq 3$ for integers $k_{i},l_{i}\geq 0$ for $i=1,\ldots,d-1$. Moreover, the payoffs at any given time step are assumed to be independent of the payoffs at all other time steps.

Following interactions within groups, each individual accumulates a payoff $P$ that is translated into a reproductive fitness that takes the form $f=1+P$. Let $f_{C}(x)$ and $f_{D}(x)$ be the reproductive fitnesses of a cooperator and a defector, respectively, at a given time step when the frequency of $C$ in the population is $x=i/N$. Since groups of size $d$ are formed at random, we have

	$\displaystyle f_{C}(x)$	$\displaystyle=1+\sum_{k=0}^{d-1}\frac{\binom{i-1}{k}\binom{N-i}{d-k-1}}{\binom{N-1}{d-k-1}}a_{k}=1+P_{C}(x),$		(3a)
	$\displaystyle f_{D}(x)$	$\displaystyle=1+\sum_{k=0}^{d-1}\frac{\binom{i}{k}\binom{N-i-1}{d-k-1}}{\binom{N-1}{d-k-1}}b_{k}=1+P_{D}(x),$		(3b)

where $P_{C}(x)$ and $P_{D}(x)$ are the average payoffs to $C$ and $D$, respectively.

The update of the population from one time step to the next is done through a single birth-death event according to a Moran model allowing for mutation (Moran [41], Ewens [9]). With probability proportional to its reproductive fitness, an individual is chosen to produce an offspring. With probability $1-u<1$, the offspring is an exact copy of its parent and, therefore, uses the same strategy. With the complementary probability $u>0$, the offspring is a mutant, in which case its strategy is chosen at random. More precisely, a mutant offspring adopts strategy $C$ with probability $1/2$ or strategy $D$ with probability $1/2$. In all cases, the offspring produced replaces an individual that is chosen at random in the population to die, possibly the parent of the offspring but not the offspring itself. This leads to the population state at the next time step.

The state space for the frequency of $C$ in the population at a given time step, represented by $X$, is $S=\{0,1/N,\ldots,(N-1)/N,1\}$. The frequency of $C$ over all time steps is an aperiodic irreducible Markov chain on a finite state space. Owing to the ergodic theorem (see, e.g., Karlin and Taylor [27]), the chain tends to an equilibrium state, called the selection-mutation equilibrium, given by a unique stationary probability distribution, represented by $\{\Pi^{\delta}(x)\}_{x\in S}$ where $\Pi^{\delta}(x)=\mathbb{P}^{\delta}(X=x)>0$ and $\sum\limits_{x\in S}\Pi^{\delta}(x)=1$.

Let $\mathbb{E}^{\delta}$ denote the expectation with respect to the stationary probability distribution if the intensity of selection is $\delta\geq 0$. The average abundance of $C$ is defined as

$$\mathbb{E}^{\delta}[X]=\sum_{x\in\mathbf{S}}x\Pi^{\delta}(x).$$

(4)

We say that weak selection favors the abundance of $C$ if its average abundance under weak enough selection exceeds what it would be under neutrality, that is,

$$\mathbb{E}^{\delta}[X]>\mathbb{E}^{0}[X]=\frac{1}{2}$$

(5)

for $\delta>0$ small enough. Here, $\mathbb{E}^{0}[X]$ represents the average abundance of $C$ under neutrality when $\delta=0$.

From one time step to the next, the frequency of $C$ can increase by $1/N$, decrease by $1/N$, or remain the same. Let us denote this change by $\Delta X=X^{\prime}-X$. We have $\Delta X=-1/N$ when a $D$ offspring replaces a $C$ individual. Let $T^{-}(x)$ be the conditional probability that $\Delta X=-1/N$ given that $X=x$. Then, we have

	$\displaystyle T^{-}(x)$	$\displaystyle=\mathbb{P}^{\delta}\left[\Delta X=-\frac{1}{N}\,\Big{|}\,X=x\right]$
		$\displaystyle=\Bigg{[}(1-u)E\left[\frac{(1-x)f_{D}(x)}{xf_{C}(x)+(1-x)f_{D}(x)}\right]+\frac{u}{2}\Bigg{]}x,$		(6)

where $E$ denotes an expectation with respect to the probability distribution of the payoffs. Similarly, $\Delta X=1/N$ when a $C$ offspring replaces a $D$ individual, and the conditional probability of this event given that $X=x$ is

	$\displaystyle T^{+}(x)$	$\displaystyle=\mathbb{P}^{\delta}\left[\Delta X=\frac{1}{N}\,\Big{|}\,X=x\right]$
		$\displaystyle=\Bigg{[}(1-u)E\left[\frac{xf_{C}(x)}{xf_{C}(x)+(1-x)f_{D}(x)}\right]+\frac{u}{2}\Bigg{]}(1-x).$		(7)

Accordingly, the conditional expected change in the frequency of $C$ is

	$\displaystyle\mathbb{E}_{x}^{\delta}\left[\Delta X\right]$	$\displaystyle=\mathbb{E}^{\delta}\left[\Delta X\,|\,X=x\right]$
		$\displaystyle=\frac{1}{N}\left(T^{+}(x)-T^{-}(x)\right)$
		$\displaystyle=\frac{u(1-2x)}{2N}+\frac{\delta}{N}(1-u)x(1-x)m(x)+o(\delta),$		(8)

where

$$m(x)=E\left[\frac{f_{C}(x)-f_{D}(x)}{xf_{C}(x)+(1-x)f_{D}(x)}\right].$$

(9)

Multiplying both sides in (3) by $\Pi^{\delta}(x)$ and summing up over all states $x$ in $S$ lead to

$$\mathbb{E}^{\delta}\left[\Delta X\right]=\frac{u}{2N}\left[1-2\mathbb{E}^{\delta}[X]\right]+\frac{\delta}{N}(1-u)\mathbb{E}^{\delta}\left[X(1-X)m(X)\right]+o(\delta).$$

(10)

In the stationary state, the frequency of $C$ in the population keeps a constant expected value, that is,

$$\mathbb{E}^{\delta}\left[\Delta X\right]=0.$$

(11)

Therefore, (10) yields

$$\mathbb{E}^{\delta}[X]=\frac{1}{2}+\frac{\delta(1-u)}{u}\mathbb{E}^{\delta}\left[X(1-X)m(X)\right]+o(\delta).$$

(12)

Then, using

$$\mathbb{E}^{\delta}\left[X(1-X)m(X)\right]=\mathbb{E}^{0}\left[X(1-X)m(X)\right]+O(\delta),$$

(13)

we obtain the first-order approximation

$$\mathbb{E}^{\delta}[X]\approx\frac{1}{2}+\frac{\delta(1-u)}{u}\mathbb{E}^{0}\left[X(1-X)m(X)\right]$$

(14)

for the expected frequency of $C$ in the selection-mutation equilibrium. Returning to (5), we have that weak selection favors the abundance of $C$ if

$$\mathbb{E}^{0}\left[X(1-X)m(X)\right]>0.$$

(15)

In the remainder of this paper, we consider a large population size, in which case we have

	$\displaystyle m(x)$	$\displaystyle=\sum_{k=0}^{d-1}\binom{d-1}{k}x^{k}(1-x)^{d-k-1}\left(\mu_{C,k}-\mu_{D,k}\right)$
		$\displaystyle\quad+\sum_{k,l=0}^{d-1}\binom{d-1}{k}\binom{d-1}{l}x^{k+l}(1-x)^{2d-k-l-2}$
		$\displaystyle\quad\quad\times\Big{[}x(\sigma_{CD,kl}-\sigma_{CC,kl})+(1-x)(\sigma_{DD,kl}-\sigma_{CD,kl})\Big{]}.$		(16)

See Appendix $A$.

Now, let us define

$$\psi_{n}^{k}=\mathbb{E}^{0}\left[X^{k}(1-X)^{n-k}\right]$$

(17)

for $k=0,1,\ldots,n$. This is the probability that among $n$ individuals drawn at random with replacement in a neutral population at equilibrium, $k$ of them in particular are of type $C$, while the other $n-k$ individuals are of type $D$. Note that $\psi^{k}_{n}=\psi^{n-k}_{n}$ since mutation is symmetric.

In a neutral population in the limit of a large population size $N$ with $N^{2}/2$ birth-death events as unit of time, each pair of lineages coalesces backward in time at rate $1$ independently of all others according to Kingman’s coalescent (Kingman[28]). Besides, each lineage mutates at the scaled rate $\theta=\lim_{N\rightarrow\infty}Nu/2$ independently of all others and independently of the coalescent process (see, e.g., Ewens[9], p. 340). This implies that the expected value in (17) in the limit of a large population size corresponds to a moment of a Dirichlet distribution (Ewens[9], p. 195). This leads to the following key lemma (see Appendix B for a proof).

Note that we have the approximation

$\displaystyle\psi_{n}^{k}=\left\{\begin{array}[]{ll}\frac{\theta(k-1)!(n-k-1)!}{2(n-1)!}&\mbox{if }k=1,2,\ldots,n-1,\\ \frac{1}{2}-\frac{\theta}{2}H_{n-1}&\mbox{if }k=0,n,\end{array}\right.$

(21)

when $\theta$ is small and terms of order $o(\theta)$ are neglected, and the approximation

$$\psi_{n}^{k}\approx\frac{1}{2^{n}}$$

(22)

for $k=0,1,\ldots,n$ when $\theta$ is large and terms of order $O(\theta^{-1})$ are neglected. Here, $H_{k}=\sum_{i=1}^{k}1/i$ denotes the $k$-th harmonic number (Conway and Guy [4]).

Using (3) and (17), the expected value on the left-hand side of (15) can be expressed as

	$\displaystyle\mathbb{E}^{0}\left[X(1-X)m(X)\right]$	$\displaystyle=\sum_{k=0}^{d-1}\binom{d-1}{k}\psi_{d+1}^{k+1}\left(\mu_{C,k}-\mu_{D,k}\right)+\sum_{k,l=0}^{d-1}\binom{d-1}{k}\binom{d-1}{l}$
		$\displaystyle\quad\times\left[-\psi_{2d+1}^{k+l+2}\left(\sigma_{CC,kl}-\sigma_{CD,kl}\right)+\psi_{2d+1}^{k+l+1}\left(\sigma_{DD,kl}-\sigma_{CD,kl}\right)\right].$		(23)

In this expression, the coefficient of $\sigma_{CC,kl}$ is always negative, while the coefficient of $\sigma_{DD,kl}$ is always positive, for $k,l=0,1,\ldots,d-1$. As for the effect of a change in any covariance between a payoff to $C$ and a payoff to $D$, $\sigma_{CD,kl}$, on average abundance, it depends on the sign of $\psi_{2d+1}^{k+l+2}-\psi_{2d+1}^{k+l+1}$, for $k,l=0,1,\ldots,d-1$. Owing to the above lemma, we have the approximation

$$\psi_{2d+1}^{k+l+2}-\psi_{2d+1}^{k+l+1}\approx\frac{(k+l+1-d)\prod_{i=1}^{k+l+2}(\theta+i-1)\prod_{j=1}^{2d-k-l-1}(\theta+j-1)}{\prod_{l=1}^{2d+1}(2\theta+l-1)}$$

(24)

if the population size is large enough. This approximation is positive if $k+l>d-1$, negative if $k+l<d-1$, and null if $k+l=d-1$.

Therefore, the following conclusion ensues.

In this section, we focus on a particular case where the scaled mean payoffs to $C$ and $D$ are given by

	$\displaystyle\mu_{C,k}$	$\displaystyle=\frac{k}{d-1}\mu_{b}-\mu_{c},$		(26a)
	$\displaystyle\mu_{D,k}$	$\displaystyle=\frac{k}{d-1}\mu_{b},$		(26b)

respectively, for $k=0,1,\ldots,d-1$. This corresponds to a public goods game where a cooperator pays a mean scaled cost $\mu_{c}$, while the scaled mean benefit of cooperation $\mu_{b}$ is distributed equally among the other $d-1$ members of the group.

Under constant payoffs, that is, $\sigma_{S_{1}S_{2},kl}=0$ for all $S_{1},S_{2}=C,D$ and $k=0,1,\ldots,d-1$, weak selection never favors the abundance of $C$ since then the espression (3) reduces to the first summation given by

$$\sum_{k=0}^{d-1}\binom{d-1}{k}\psi_{d+1}^{k+1}\left(\mu_{C,k}-\mu_{D,k}\right)=-\mu_{c}\sum_{k=0}^{d-1}\binom{d-1}{k}\psi_{d+1}^{k+1}=-\mu_{c}\psi_{2}^{1},$$

(27)

which cannot be positive. Here, we have used the identity

	$\displaystyle\sum_{k=0}^{n}\binom{n}{k}\psi_{n+i}^{k+j}$	$\displaystyle=\sum_{k=0}^{n}\binom{n}{k}\mathbb{E}^{0}\left[X^{k+j}(1-X)^{n+i-k-j}\right]$
		$\displaystyle=\mathbb{E}^{0}\left[X^{j}\left(1-X\right)^{i-j}\sum_{k=0}^{n}\binom{n}{k}X^{k}\left(1-X\right)^{n-k}\right]$
		$\displaystyle=\mathbb{E}^{0}\left[X^{j}\left(1-X\right)^{i-j}\right]=\psi^{j}_{i},$		(28)

for $0\leq j\leq i$.

In the remainder of this section, we will show that introducing uncertainty in the payoffs to $D$ can make it possible for weak selection to favor the abundance of $C$.

In this section and the next two ones, we are interested in classical social dilemmas where cooperation incurs a random cost $c>0$ but provides a random benefit $b>c$ in groups of size $d$. Moreover, we assume that

	$\displaystyle E[b]$	$\displaystyle=\mu_{b}\delta+o(\delta),$		(37a)
	$\displaystyle E[c]$	$\displaystyle=\mu_{c}\delta+o(\delta),$		(37b)
	$\displaystyle E[b^{2}]$	$\displaystyle=\sigma^{2}_{b}\delta+o(\delta),$		(37c)
	$\displaystyle E[c^{2}]$	$\displaystyle=\sigma^{2}_{c}\delta+o(\delta),$		(37d)
	$\displaystyle E[bc]$	$\displaystyle=\sigma_{bc}\delta+o(\delta)$		(37e)

and

$$E\left[b^{i}c^{j}\right]=o(\delta),$$

(38)

for any non-negative integers $i$ and $j$ such that $i+j\geq 3$.

We consider first a linear public goods game in which the benefit of cooperation $b$ by an individual at a cost $c$ is distributed equally among the $d-1$ other individuals in the same group and all effects of cooperation are additive (Hamburger [18], Fox and Guyer [10], Nowak and Sigmund [42], Wild and Traulsen [60]). In this case, the payoffs to $C$ and $D$ for an individual in interaction with $k$ cooperators and $d-k-1$ defectors are

	$\displaystyle a_{k}$	$\displaystyle=\frac{k}{d-1}b-c,$		(39a)
	$\displaystyle b_{k}$	$\displaystyle=\frac{k}{d-1}b,$		(39b)

whose scaled means are given by

	$\displaystyle\mu_{C,k}$	$\displaystyle=\frac{k}{d-1}\mu_{b}-\mu_{c},$		(40a)
	$\displaystyle\mu_{D,k}$	$\displaystyle=\frac{k}{d-1}\mu_{b},$		(40b)

and scaled variances and covariances by

	$\displaystyle\sigma_{CC,kl}$	$\displaystyle=\frac{kl}{(d-1)^{2}}\sigma^{2}_{b}+\sigma^{2}_{c}-\frac{k+l}{d-1}\sigma_{bc},$		(41a)
	$\displaystyle\sigma_{CD,kl}$	$\displaystyle=\frac{kl}{(d-1)^{2}}\sigma^{2}_{b}-\frac{k}{d-1}\sigma_{bc},$		(41b)
	$\displaystyle\sigma_{DD,kl}$	$\displaystyle=\frac{kl}{(d-1)^{2}}\sigma^{2}_{b},$		(41c)

for $k,l=0,1,\ldots,d-1$. Then, the first summation in (3) is the same as in (27). On the other hand, using the identity $\frac{l}{d-1}\binom{d-1}{l}=\binom{d-2}{l-1}$, we obtain

	$\displaystyle\sum_{k,l=0}^{d-1}\binom{d-1}{k}\binom{d-1}{l}\frac{l}{d-1}\psi_{2d+1}^{k+l+i}$
	$\displaystyle=\sum_{k=0}^{d-1}\sum_{l=1}^{d-1}\binom{d-1}{k}\binom{d-2}{l-1}\mathbb{E}^{0}\left[X^{k+l+i}\left(1-X\right)^{2d+1-k-l-i}\right]$
	$\displaystyle=\mathbb{E}^{0}\left[X^{i+1}(1-X)^{3-i}\sum_{k=0}^{d-1}\binom{d-1}{k}X^{k}\left(1-X\right)^{d-1-k}\sum_{l=0}^{d-2}\binom{d-2}{l}X^{l}\left(1-X\right)^{d-2-l}\right]$
	$\displaystyle=\mathbb{E}^{0}\left[X^{i+1}(1-X)^{3-i}\right]=\psi_{4}^{i+1},$		(42)

for $i=0,1,2,3$. Then, the second summation in (3) can be written as

$$\begin{split}&\sum_{k,l=0}^{d-1}\binom{d-1}{k}\binom{d-1}{l}\left[\left(\frac{l}{d-1}\psi_{2d+1}^{k+l+2}+\frac{k}{d-1}\psi_{2d+1}^{k+l+1}\right)\sigma_{bc}-\psi_{2d+1}^{k+l+2}\sigma^{2}_{c}\right]\\ &=\left(\psi^{2}_{4}+\psi^{3}_{4}\right)\sigma_{bc}-\psi^{2}_{3}\sigma^{2}_{c}\\ &=\psi^{2}_{3}\left(\sigma_{bc}-\sigma^{2}_{c}\right).\end{split}$$

(43)

In the last passage, we have used the identity

$$\psi^{i}_{n}+\psi^{i+1}_{n}=\psi^{i}_{n-1}.$$

(44)

Inserting (27) and (43) in (14) and using Lemma $1$, the first-order approximation of the average abundance of $C$ is given by

$$\mathbb{E}^{\delta}[X]\approx\frac{1}{2}+\frac{\delta N}{8(2\theta+1)}\left(\sigma_{bc}-\sigma^{2}_{c}-2\mu_{c}\right).$$

(45)

In a stag hunt game, an individual in a group of size $d$ receives a benefit $b$ only if all the individuals in the group cooperate, each one at a cost $c$ (Skyrms [53], Pacheco et al. [45]). Then, the payoffs to $C$ are given by $a_{k}=-c$ if $k=0,1,\cdots,d-2$ and $a_{d-1}=b-c$, while the payoffs to $D$ are $b_{l}=0$ for $l=0,1,\ldots,d-1$. In this case, the scaled means of the payoffs to $C$ or $D$ according to the numbers of cooperating partners in the same group are given by

$$\mu_{C,k}=\left\{\begin{array}[]{ll}-\mu_{c}&\mbox{if }k<d-1,\\ \mu_{b}-\mu_{c}&\mbox{if }k=d-1,\end{array}\right.$$

(47)

and $\mu_{D,l}=0$, and the scaled variances and covariances by

$$\sigma_{CC,kl}=\left\{\begin{array}[]{ll}\sigma^{2}_{c}&\mbox{if }k<d-1\mbox{ and }l<d-1,\\ \sigma^{2}_{c}-\sigma_{bc}&\mbox{if }k=d-1\mbox{ and }l<d-1\mbox{ or }k<d-1\mbox{ and }l=d-1,\\ \sigma^{2}_{c}-2\sigma_{bc}+\sigma^{2}_{b}&\mbox{if }k=l=d-1,\end{array}\right.$$

(48)

and $\sigma_{CD,kl}=\sigma_{DD,kl}=0$, for $k,l=0,1,\ldots,d-1$.