Fractional Punishment of Free Riders to Improve Cooperation in Optional Public Good Games

Abstract

Improving and maintaining cooperation are fundamental issues for any project to be time-persistent, and sanctioning free riders may be the most applied method to achieve it. However, the application of sanctions differs from one group (project or institution) to another. We propose an optional, public good game model where a randomly selected set of the free riders is punished. To this end, we introduce a parameter that establishes the portion of free riders sanctioned with the purpose to control the population state evolution in the game. This parameter modifies the phase portrait of the system, and we show that, when the parameter surpasses a threshold, the full cooperation equilibrium point becomes a stable global attractor. Hence, we demonstrate that the fractional approach improves cooperation while reducing the sanctioning cost.

1. Introduction

Improving cooperation is a fundamental issue in human activities. The emergence and evolution of cooperation have been studied in life and social science; for instance, biologists have vastly studied altruism [1,2,3], and economists have studied social dilemmas [4,5]. In many cases, this problem has been modeled by using game theory [6,7] and, lately, with the development of evolutionary game theory [8] by using evolutionary dynamics models [9,10,11].

In the literature, two main mechanisms that affect cooperation are found: incentives and abstention. Incentives are known to improve the rate of cooperative behavior in a group [12,13,14,15,16,17], while abstention to participate in a project is useful to maintain cooperation [18,19]. However, incentives are expensive; for this reason, a punishment system is a public good itself [12].

In this paper, we consider improving cooperation by sanctioning only a fraction of the free riders, denominated from here onwards, as fractional punishment. This approach seeks to reduce the number of free riders while minimizing the cost of the sanctioning system. The punishment is built as an extension of the optional public good game model presented in [18].

Therefore, the optional public good game that will be analyzed in this paper consists of a set of randomly selected individuals from the population who shares an institutional service. Each of the individuals decides among three strategies: refuse to play, decide to play and contribute to the common pool (pay for the service), and decide to play but do not contribute (use the service but do not pay for it). In the public good game, the benefit from the common pool is divided equally among all the players. In this work, a fraction of the free riders are randomly selected, and their payoff is reduced to zero.

Sanctioning System

In public good games, the two well-studied forms of sanctioning free riders are peer punishment and pool punishment [20]. In peer punishment, after the game, an individual that contributed (a cooperator) can decide to use part of its payoff to sanction all the free riders in the population. On the contrary, pool punishment is defined before the game; each cooperator must decide to contribute or not to a second pool that will be later on used to sanction the free riders.

A cooperator that is willing to sanction free riders plays a punisher strategy. A punisher has to contribute to the game and to the sanctioning system. A cooperator that decides not to sanction free riders becomes themselves a free rider for the sanctioning system (called second-order free riders). With peer punishment, the cost for a small number of punishers in a group with mostly defectors prevents the punishers from prospering; conversely, in a group without defectors, cooperators and punishers have the same payoff. Therefore, punishers tend to disappear, allowing defectors space to invade the population. With pool punishment, this cannot occur since the punishers must contribute ex ante to the game. Pool punishment, however, is a fixed cost that exists even when there are no defectors in the population [20,21].

Both sanctioning methods were studied and compared in theoretical models and experiments [13,20,22,23]. For instance, it was found that pool punishment can stabilize cooperation if second-order free riders are also punished [21,22]. In addition, the second-order free-riding problem was analyzed in [23,24,25].

Peer and pool punishment are rarely applied in their pure form in practical situations. More recently, some authors have explored different ways to penalize free riders often inspired in real-life situations. In [26], the sanction that a defector has to pay is a fixed amount. Avoiding second-order free riding by introducing an entrance fee was analyzed in [15]. A set of the cooperators was randomly selected to perform as punishers in [27], the cost of punishing was shared between cooperators in [28], and a ceiling payoff for defectors was introduced to promote cooperation in [29,30].

In this work, we assume that the sanctions have to be implemented with the fee collected periodically in a fee pool. This is observed as a common practice in institutions that provide a service. Within this framework, second-order free-riding does not occur since those who pay for the service implicitly are paying for the sanctioning system. However, this also implies that the profit that each player receives is reduced, given that a part of the collected resources must be used to penalize free riders.

The institution decides the fee-pool portion that will be used for each expense: providing the service and sanctioning free riders. The sanctioning system may not be a fixed cost; when free-riding decreases, the fee-pool amount required to sanction free riders will also decrease. Besides, the institution may limit sanctioning expenses by penalizing only a fraction of the free riders.

The aim of this paper consists in modeling mathematically a fractional punishment mechanism commonly applied in real situations to reduce expenses. We prove that fractional punishment improves cooperation and reduces the cost of punishing, and full cooperation can be achieved when the fraction of sanctioned defectors surpasses a minimum. Moreover, depending on the size of the punished fraction, the cooperation can be attained faster. This means that, by adjusting the parameter, the fractional punishment modifies the nature of the equilibrium points of the system and the phase portrait itself. With this analysis, decision-makers of the institutions can produce scenarios to foresee the effect of the fraction sanctioned in the improvement of the cooperation.

This article is organized as follows: in Section 2, the proposed model and dynamic, as well as the definition of the payoff, are presented. In Section 3, the system is analyzed, and the equilibrium points are characterized as a function of the punishment parameter denoted by d. The implications of the model, along with the penalization of a limited set of the free riders, are discussed in Section 4, while some concluding remarks are finally presented in Section 5.

2. Materials and Methods

Consider a large population where from time to time n randomly selected individuals are offered to play in a public good game ($\mathit{n}={\mathit{n}}_{\mathit{x}}+{\mathit{n}}_{\mathit{y}}+{\mathit{n}}_{\mathit{z}}$; $\mathit{n}\ge 2$). Those who refuse to play are known in the game as loners (${\mathit{n}}_{\mathit{z}}$). The loners have a payoff $\sigma$ that is constant and independent of the game. Those who decide to play the public good game will receive their payoff from the game ($s=n_x+n_y$). In the game, each player must decide if either contributes or not with a positive amount c to the common pool. Individuals who contribute are cooperators ($n_x$), and those who do not are defectors ($n_y$). The total contribution to the common pool (obtained from cooperators) is then multiplied by a factor r$(1<r<n)$; in a public good game, the amount in the common pool is distributed equally between all the players, regardless of whether they are cooperators or defectors. Instead, in this work, a randomly selected set of defectors will be sanctioned reducing their payoff to zero.

For modeling the population dynamics, in this paper and following [11], every individual in the population follows a strategy i (for $i=x,y,z$) from the three options available: cooperators x, defectors y, and loners z.

Let $0\le x\left(t\right),y\left(t\right),\mathrm{and}z\left(t\right)\le 1$ be the frequency of each of the corresponding available strategies of the population at a specific time t. To simplify the notation, we drop henceforth the time dependency t and simply write x, y, and z. The frequency distribution of the whole population at a specific time t is defined by the state $[x,y,z]$, which belongs to the unit simplex ${\mathcal{S}}_3$ given by

$${\mathcal{S}}_3=\left\{[x,y,z]\in {\mathbb{R}}^3:x,y,z\ge 0\mathrm{and}x+y+z=1\right\}.$$

(1)

The interior of ${\mathcal{S}}_3$ is defined as the set of points where all strategies are present (i.e., x > 0, y > 0, and z > 0) and the boundary of ${\mathcal{S}}_3$ is defined as the set of points that are not interior points. The boundary of the simplex comprises vertices and borders. The vertices x, y, and z represent a homogeneous population of cooperators, defectors, and loners, respectively. In the borders, $xy$, $yz$, and $zx$, one strategy is absent.

Assuming a large enough population in which generations blend continuously into each other, each strategy evolves with a rate that expresses its evolutionary success following the replicator equation

$$\left\{\begin{array}{c}\dot{x}=x\left(p_x-\overline{p}\right)\hfill \\ \dot{y}=y\left(p_y-\overline{p}\right)\hfill \\ \dot{z}=z\left(p_z-\overline{p}\right)\hfill \end{array}\right.$$

(2)

with $x\left(0\right)$, $y\left(0\right)$, and $z\left(0\right)$ as initial conditions in the simplex ${\mathcal{S}}_3$, with $\overline{p}=x{p}_x+y{p}_y+z{p}_z$, where $p_i$ is the payoff of the ith strategy, and $\overline{p}$ is the average payoff of the population.

It is important to remark that the simplex ${\mathcal{S}}_3$ remains invariant under the flow of Equation (2) (see [11,31]); hence, if the state $[x,y,z]\in {\mathcal{S}}_3$ at $t=0$, then $[x,y,z]\in {\mathcal{S}}_3$$\forall t$, and it is endowed with the standard norm ${\parallel ·\parallel }_2$. Notice that, in Expression (2), an important behavior relies on the payoff $p_i$ of the ith strategy with respect to the average population payoff $\overline{p}$ in the sense that, if $p_i>\overline{p}$, then the proportion of the i strategy increases, otherwise it decreases. Next, the payoff of each strategy is defined.

Considering a group of $n_x$ cooperators, $n_y$ defectors, and $n_z$ loners and assuming that each contribution is equal $(c=1)$, the payoffs of each strategy in an optional public good game without punishment is defined by [18,32]:

$$p_x=r{\displaystyle \frac{n_x}{s}}-1,p_y=r{\displaystyle \frac{n_x}{s}},p_z=\sigma ,$$

(3)

where the parameters r, $\sigma$, and n are considered under the following assumptions:

Assumption 1.

The interest rate on the common pool r satisfies $1<r<n$.

If $1<r$, i.e., if all individuals cooperate, they are better off than if all defects. If $r<n$, each individual is better off defecting than cooperating [18].

Assumption 2.

The payoff of the loner strategy σ satisfies $0<\sigma <r-1$.

This means that the revenue of a cooperator in a group of cooperators (profiting $(r-1)$ each) is better off than loners that receive $\sigma$, but loners are better off than a defector in a group of defectors, where the payoff is equal to 0 [18].

To compute the payoff of a player in a group, consider that this group composition depends on the frequencies of all strategies in the population. In the sample group of size n, where s$(s\le n)$ is the number of those who will participate (cooperators and defectors), and $n-s$ is those who refuse (loners), if $s=1$, the game cannot take place, and the payoff will be equal to a loner payoff regardless of the strategy of the player, and this happens with a probability $z^{n-1}$. If $s\ge 2$, an individual that agrees to play has $s-1$ coplayers from the $n-1$ sample given by the probability [18]

$$\left(\genfrac{}{}{0pt}{}{n-1}{s-1}\right){(1-z)}^{s-1}{z}^{n-s}.$$

(4)

Furthermore, the probability that j of these $s-1$ coplayers will be cooperators and the other $s-1-j$ defectors is [18]

$$\left(\genfrac{}{}{0pt}{}{s-1}{j}\right){\left(\frac{x}{x+y}\right)}^j{\left(\frac{y}{x+y}\right)}^{s-1-j}.$$

(5)

In this work, we consider that only a fraction d ($0\le d\le 1$) of the defectors will be punished, modifying the payoff presented in (3). In particular, we assume that this set of randomly selected defectors will have their corresponding payoff reduced to 0, while the remaining free riders will obtain the normal payoff. Therefore, a defector not knowing if they will be punished or not, in the presence of j cooperators, will have a payoff equal to

$$p_y=(1-d)\left({\displaystyle \frac{rj}{s}}\right)+d0.$$

(6)

The expected payoff for a defector in a group of s$(s=2,...,n)$ players over all possible numbers of cooperators is

$$\begin{array}{cc}\hfill & (1-d)\frac{r}{s}\sum _{j=0}^{s-1}j\left(\genfrac{}{}{0pt}{}{s-1}{j}\right){\left(\frac{x}{x+y}\right)}^j{\left(\frac{y}{x+y}\right)}^{s-1-j}=(1-d)\frac{r}{s}(s-1)\frac{x}{x+y}.\hfill \end{array}$$

Taking into consideration that a defector can be alone $(s=1)$ or in a group $(s\ge 2)$ and that both the number or participants (s) and cooperators $\left(j\right)$ in the group are random variables, the average payoff for a defector over all possible numbers of participants is given by

$$\begin{array}{cc}\hfill p_y(x,z,d)& =\sigma z^{n-1}+\left(1-d\right)\left(r\frac{x}{1-z}\sum _{s=1}^n\left(\genfrac{}{}{0pt}{}{n-1}{s-1}\right){(1-z)}^{s-1}{z}^{n-s}\left(1-\frac{1}{s}\right)\right).\hfill \end{array}$$

Observing that $\left(\genfrac{}{}{0pt}{}{n-1}{s-1}\right)=\left(\genfrac{}{}{0pt}{}{n}{s}\right)\frac{s}{n}$, then the expression for $p_y$ can be rewritten as [18]

$$\begin{array}{cc}\hfill p_y(x,z,d)& =\sigma z^{n-1}+\left(1-d\right)\left(r\frac{x}{1-z}\left(1-\frac{1}{n}\sum _{s=1}^n\left(\genfrac{}{}{0pt}{}{n}{s}\right){(1-z)}^{s-1}{z}^{n-s}\right)\right)\hfill \\ \hfill & =\sigma z^{n-1}+\left(1-d\right)\left(r\frac{x}{1-z}\left(1-\frac{1}{n}\frac{1-z^n}{(1-z)}\right)\right).\hfill \end{array}$$

For the sake of simplicity, we define an auxiliary function $a:=(1-z^n)/\left(n\left(1-z\right)\right)$. The payoff then finally takes the following form:

$$\begin{array}{cc}\hfill p_y(x,z,d)& =\sigma z^{n-1}+\left(1-d\right)r\frac{x}{1-z}\left(1-a\right).\hfill \end{array}$$

(7)

A similar analysis can be performed to obtain the expected payoff of a cooperator in the population. To distinguish the contribution from one cooperator to the contribution of the other cooperators in the group, the payoff can be rewritten:

$$p_x=\left({\displaystyle \frac{r(n_x-1)}{n}}+{\displaystyle \frac{r}{n}}\right)-1.$$

where $n_x-1$ (the remaining cooperators in the game) is defined as j.

$$\begin{array}{c}\hfill p_x=\left({\displaystyle \frac{rj}{n}}+{\displaystyle \frac{r}{n}}\right)-1=\left({\displaystyle \frac{r(j+1)}{n}}\right)-1.\end{array}$$

Therefore, the expected payoff for a cooperator in a group of $s(s=2,...,n)$ players over all possible numbers of cooperators is

$$\begin{array}{cc}\hfill & \sum _{j=0}^{s-1}\left(\genfrac{}{}{0pt}{}{s-1}{j}\right)\left(\frac{r}{s}(j+1)-1\right){\left(\frac{x}{x+y}\right)}^j{\left(\frac{y}{x+y}\right)}^{s-1-j}=\frac{r}{s}(s-1)\frac{x}{1-z}+{\displaystyle \frac{r}{s}}-1.\hfill \end{array}$$

If the cooperator is the only one that wants to play, that is, if $s=1$, the payoff is $\sigma z^{n-1}$. If $(s\ge 2)$, the payoff depends on the number or participants (s) in the group. Hence, the expected payoff for a cooperator is given by

$$\begin{array}{cc}\hfill p_x(x,z)& =\sigma z^{n-1}+\sum _{s=2}^n\left(\genfrac{}{}{0pt}{}{n-1}{s-1}\right){(1-z)}^{s-1}{z}^{n-s}\left(\frac{r}{s}(s-1)\frac{x}{1-z}+{\displaystyle \frac{r}{s}}-1\right)\hfill \\ \hfill p_x(x,z)& =\sigma z^{n-1}+r\frac{x}{1-z}\left(\sum _{s=2}^n\left(\genfrac{}{}{0pt}{}{n-1}{s-1}\right){(1-z)}^{s-1}{z}^{n-s}\right)-r\frac{x}{1-z}\left(\sum _{s=2}^n\left(\genfrac{}{}{0pt}{}{n-1}{s-1}\right){(1-z)}^{s-1}{z}^{n-s}\frac{1}{s}\right)\hfill \\ \hfill & +r\left(\sum _{s=2}^n\left(\genfrac{}{}{0pt}{}{n-1}{s-1}\right){(1-z)}^{s-1}{z}^{n-s}\frac{1}{s}\right)-\sum _{s=2}^n\left(\genfrac{}{}{0pt}{}{n-1}{s-1}\right){(1-z)}^{s-1}{z}^{n-s}.\hfill \end{array}$$

Observing that $\left(\genfrac{}{}{0pt}{}{n-1}{s-1}\right)=\left(\genfrac{}{}{0pt}{}{n}{s}\right)\frac{s}{n}$, we obtain

$$\begin{array}{cc}\hfill p_x(x,z)& =\sigma z^{n-1}+r\frac{x}{1-z}(1-z^{n-1})-r\frac{x}{1-z}\left(\frac{1}{n}\frac{(1-z^n)}{(1-z)}-z^{n-1}\right)\hfill \\ \hfill & +r\left(\frac{1}{n}\frac{(1-z^n)}{(1-z)}-z^{n-1}\right)-1+z^{n-1}.\hfill \end{array}$$

Rearranging and using, as before, the auxiliary function a, we obtain

$$\begin{array}{cc}\hfill p_x(x,z)& =\sigma z^{n-1}+ra+r\frac{x}{1-z}\left(1-a\right)+(1-r)z^{n-1}-1.\hfill \end{array}$$

(8)

The difference in payoff between defector and cooperator strategies shows the relative benefit (or drawback) of defectors over cooperators, and it is essential to characterize the solution orbits. This difference is given by

$$\begin{array}{cc}\hfill g(x,z,d):=& p_y-p_x=1+(r-1)z^{n-1}-\frac{r}{n}\frac{(1-z^n)}{(1-z)}-dr\frac{x}{(1-z)}\left(1-\frac{1}{n}\frac{1-z^n}{(1-z)}\right)\hfill \\ \hfill :=& 1+(r-1)z^{n-1}-ra-dr\frac{x}{(1-z)}\left(1-a\right).\hfill \end{array}$$

(9)

In order to highlight the effect of the parameter d, Expression (9) is rewritten as

$$\begin{array}{cc}\hfill g(x,z,d)& :=m\left(z\right)-dxh\left(z\right),\hfill \end{array}$$

(10)

where $m\left(z\right)=1+(r-1)z^{n-1}-ra$ and $h\left(z\right)=\frac{r}{1-z}\left(1-a\right)$. Observe that, when $d=0$, then $g(x,z,d)=m\left(z\right)$, and the free system introduced in [18] is recovered.

3. Results

In this section, we analyze the system (2). To help the reader, and at the risk of oversimplification, before starting the analysis, we outline the effect of the parameter d on the phase portrait of the system (2). This is presented in Figure 1 for increasing values of d. Note that this system has three obvious equilibrium points, which are the vertices of the simplex: (i) $(\widehat{x},\widehat{y},\widehat{z})=(1,0,0)$ (ii) $(\widehat{x},\widehat{y},\widehat{z})=(0,1,0)$, and (iii) $(\widehat{x},\widehat{y},\widehat{z})=(0,0,1)$. Another equilibrium point exists when all strategies have the same payoff $p_x=p_y=p_z=\overline{p}$, which is an interior equilibrium point (denoted as $q_1$). Later, we show that another equilibrium point arises when incrementing the value of d (this point is denoted by $q_2$).

For $d=0$ (no punishment), the vertices ($x,y,z$) are saddle points and the interior equilibrium ($q_1$) is a center (see Figure 1a), which is proved in [18]. The interesting effect of the punishment parameter d consists in the fact that the interior equilibrium point changes its stability (see Figure 1b). Furthermore, when d is large enough ($d_1$), an interior equilibrium point arises ($q_2$) in the border $xy$ (see Figure 1c). As a consequence, the nature of the equilibrium in the vertex x also changes. As d increases, the equilibrium point $q_1$ shifts to the border $xy$, and for a value $d_2$, both equilibrium points ($q_1$ and $q_2$) merge into a new equilibrium point denoted by $q_3$ (see Figure 1e). The nature of each of the aforementioned fixed points is studied in the next subsections.

For the analysis, we use the following sequence. In Section 3.1, the boundary of the simplex, particularly the border $xy$, where a new equilibrium point arises for a sufficiently large d, is analyzed. In Section 3.2, the interior of the simplex is analyzed. First, in Section 3.2.1 and Section 3.2.2, the effect of d in the position of the interior equilibrium point is considered, while in Section 3.2.3, the stability of the interior equilibrium point is analyzed. The stability of the equilibrium in the vertex x is analyzed in Section 3.3. A summary is presented in Section 3.4.

3.1. Boundary of the Simplex

In the borders, System (2) reduces to two equations. Next, we present the analysis of System (2) restricted to the boundary of ${\mathcal{S}}_3$.

3.1.1. Border $xy$

This border represents a system with cooperators and defectors, where the game is compulsory. Given that, in this border, $x+y=1$, the dynamics can be analyzed with the following equation:

$$\begin{array}{cc}\hfill \dot{x}& =-x(1-x)g(x,d)=-x(1-x)(1-(r/n)-drx(1-1/n\left)\right).\hfill \end{array}$$

(11)

The equilibrium points are $\widehat{x}=1$ ($\widehat{y}=0$), $\widehat{x}=0$ ($\widehat{y}=1$), and $\widehat{x}=(n-r)/\left(r(n-1)d\right)$. The qualitative characterization of the aforementioned equilibrium points were studied in [33] with the following lemma (see the proof in Appendix A.1).

Lemma 1.

Let us consider Equation (11) modeling the border $xy$ where x denotes the cooperators, y the defectors, $x\ge 0$, $y\ge 0$, and $x+y=1$. If the fraction of defectors being sanctioned is denoted by d, with $0\le d\le 1$ and defining $d_1:=(n-r)/\left(r(n-1)\right)$, then (1) for $0\le d\le 1$, the equilibrium point $\widehat{x}=0$ is locally asymptotically stable, (2) $\widehat{x}=1$ is locally asymptotically stable for $d_1<d\le 1$, and (3) the point $\widehat{x}=(n-r)/\left(r(n-1)d\right)$, which exists in the interior of the border for $d_1<d$, is an unstable equilibrium point.

The possible outcomes of the system presented in Lemma 3.1 are sketched in Figure 2. It is important to mention the strong dependence of the characterization of the equilibrium points on the function $g(x,d)$. If 0 $\le d<d_1$, then $g(x,d)\ne 0$ for all values of x, and the two equilibrium points are $\widehat{x}=0$ and $\widehat{x}=1$. In such a case, the $\widehat{x}=1$ is globally unstable, and $\widehat{x}=0$ is globally asymptotically stable (see Figure 2a). For $d_1<d$, an unstable equilibrium point $\tilde{x}$ arises in the interior of the border (see Figure 2b). If $d=(n-r)/r(n-1)x$, then $g(x,d)=0$, the corresponding state x is the equilibrium point $\tilde{x}$.

3.1.2. Border $zx$

This border represents a population of cooperators and loners. When $x+z=1$, it is sufficient to analyze the equation $\dot{x}=((r-1)-\sigma )x(1-x)\left(1-z^{n-1}\right)$ to characterize the system. Since $\left(r-1\right)\ne \sigma$, the equilibrium points are $\widehat{x}=1$ ($\widehat{z}=0$) and $\widehat{x}=0$ ($\widehat{z}=1$). When $\left(r-1\right)=\sigma$, the system remains constant, and every point in the border is an equilibrium point, regardless of the value of z. The Jacobian of $\dot{x}$ has the following form:

$$\begin{array}{cc}\hfill J(x,z)& =(1-2x)\left(((r-1)-\sigma )(1-z^{n-1})\right)\hfill \end{array}$$

(12)

Evaluating (12) shows that the equilibrium point $\widehat{x}=1$ ($\widehat{x}=0$) is stable (unstable) if $(r-1)>\sigma$. The dynamic of the border goes from only loners in the population to full cooperation. On the contrary, if $(r-1)<\sigma$, the equilibrium point $\widehat{x}=1$ ($\widehat{x}=0$) is unstable (stable) and the system reverses direction; it goes from full cooperation to a population of loners.

3.1.3. Border yz

This border represents a population composed of defectors and loners. Since $y+z=1$, the system can be completely represented by $\dot{y}=y(1-y)\sigma (z^{n-1}-1)$. When $\sigma \ne 0$, the equilibrium points are $\widehat{y}=1$ ($\widehat{z}=0$) and $\widehat{y}=0$ ($\widehat{z}=1$). When $\sigma =0$, the system remains constant, and every point in the border is an equilibrium point, regardless of the value of z. The Jacobian of $\dot{y}$ has the following form:

$$J(y,z)=(1-2y)\sigma (z^{n-1}-1)$$

(13)

The equilibrium point $\widehat{y}=1$ ($\widehat{y}=0$) is unstable (stable) if $\sigma >0$. The dynamic in the border goes from full defection to a population of loners. Conversely, if $\sigma <0$, $\widehat{y}=1$ ($\widehat{y}=0$) is stable (unstable) and the system reverses direction, going from only loners to full defection.

Observe that the parameter d does not affect the borders $zx$ and $yz$. Furthermore, since the game setup defines $0<\sigma <r-1$, the dynamic of the border $yz$ goes from full defection to only loners in the population, while the dynamic of the border $xz$ goes from only loners in the population to full cooperation.

3.2. Interior of the Simplex

To determine and analyze the interior equilibrium point, a new variable $f:=x/(x+y)$ is defined [18]. This variable represents the fraction of cooperators among the players in the game and allows one to define a bijection $\mathcal{T}:(x,y,z)\to (f,z)$, provided that the restriction $x+y+z=1$ is satisfied. Hence, the system (2) can be rewritten in a more convenient reduced form [18]:

$$\left\{\begin{array}{ccc}\dot{f}& =& -f(1-f)g(f,z,d)\hfill \\ \dot{z}& =& z(\sigma -\overline{p}(f,z))\hfill \end{array}\right.$$

(14)

where $(f,z)\in [0,1]\times [0,1]$. The function $g(f,z,d)$ is obtained replacing $x=f(1-z)$ in Expression (10), so $g(f,z,d)=1+(r-1)z^{n-1}-ra-drf\left(1-a\right)$. To obtain a more explicit expression than (14), consider that $y=1-x-z$, and $f=x/(1-z)$. The system of Equation (14) can then be rewritten in the following form:

$$\left\{\begin{array}{ccc}\dot{f}& =& -f(1-f)g(f,z,d)\hfill \\ \dot{z}& =& z(1-z)\left(\sigma +fg(f,z,d)-p_y(f,z,d)\right).\hfill \end{array}\right.$$

(15)

For each time t, d parametrizes the solution $(f,z)$ of the system of Equation (15). An equilibrium point of System (15) exists in the interior of the simplex if the following expressions are satisfied:

$$g(f,z,d)=0$$

(16a)

$$\sigma -fg(f,z,d)-\sigma z^{n-1}-(1-d)rf(1-a)=0,$$

(16b)

and since both equations have to be satisfied simultaneously, Equation (16b) is reduced to $\sigma (1-z^{n-1})-(1-d)rf(1-a)=0$, so the value of f can be obtained by

$$f=\frac{\sigma }{(1-d)r}\frac{(1-z^{n-1})}{(1-a)}.$$

(17)

By then introducing f in (16a), an equivalent definition for $g(·)$ that depends only on z and d is obtained:

$$\tilde{g}(z,d)=m\left(z\right)+\tilde{h}(z,d)$$

(18)

where $\tilde{h}(z,d)=-\frac{d}{1-d}\sigma (1-z^{n-1})$ and $m\left(z\right)=1+(r-1)z^{n-1}-ra$ (as introduced in Expression (10)). Now, from (18), the value of z in the equilibrium, defined as $\widehat{z}$, can be numerically obtained as a function of d. Recalling that a depends only on z (as defined in Section 2) and denoting $\widehat{a}$ as the expression of a with $\widehat{z}$, then (16a) and (16b) become

$$\left\{\begin{array}{ccc}1+(r-1){\widehat{z}}^{n-1}-r\widehat{a}-dr\widehat{f}\left(1-\widehat{a}\right)& =& 0\hfill \\ \sigma {\widehat{z}}^{n-1}+\left(1-d\right)r\widehat{f}\left(1-\widehat{a}\right)& =& \sigma \hfill \end{array}\right.$$

(19)

where $\widehat{f}$ corresponds to the value of f at the equilibrium point associated with $\widehat{z}$ for a d. Defining an auxiliary variable as $\tilde{b}:=r\widehat{a}$, the above equations take the form of a system of two equations with two unknowns ($\widehat{f}$ and $\widehat{b}$), as follows:

$$\left\{\begin{array}{ccc}1+(r-1){\widehat{z}}^{n-1}-\tilde{b}-d\widehat{f}\left(r-\tilde{b}\right)& =& 0\hfill \\ \sigma {\widehat{z}}^{n-1}+\left(1-d\right)\widehat{f}\left(r-\tilde{b}\right)& =& \sigma \hfill \end{array}\right.$$

(20)

Extracting $\widehat{f}=(\sigma (1-{\widehat{z}}^{n-1})/\left((1-d)(r-\tilde{b})\right)$ from the second equation and introducing it in the first equation, $\tilde{b}$ is obtained. Replacing $\tilde{b}$ in the second equation, a new form of $\widehat{f}$ in the equilibrium is then obtained:

$$\widehat{f}=\frac{\sigma }{(r-1)+d(\sigma -(r-1\left)\right)}.$$

(21)

Notice that the interior equilibrium point $(\widehat{f},\widehat{z})$ corresponds in the system $(\widehat{x},\widehat{y},\widehat{z})$ to the equilibrium point in the interior of the simplex ${\mathcal{S}}_3$ given by

$$\left(\widehat{x},\widehat{y},\widehat{z}\right)=\left(\frac{\sigma }{\alpha }(1-\widehat{z}),(1-\frac{\sigma }{\alpha })(1-\widehat{z}),\widehat{z}\right)$$

where $\alpha =(r-1)+d(\sigma -(r-1\left)\right)$ and the dependence on the parameter d is given explicitly. To relate both systems, a summary of the equilibrium points with respect to the variables $(x,y,z)$ and its corresponding $(f,z)$ are given in

Table 1. Equilibrium points relationship between $(\widehat{x},\widehat{y},\widehat{z})$ and $(\widehat{f},\widehat{z})$. $\alpha =(r-1)+d(\sigma -(r-1\left)\right)$, $\beta =(N-r)/\left(r\right(N-1\left)d\right)$.

$(\widehat{\mathit{x}},\widehat{\mathit{y}},\widehat{\mathit{z}})$	$(\widehat{\mathit{f}},\widehat{\mathit{z}})$
$(\widehat{x},\widehat{y},\widehat{z})=(1,0,0)$	$(\widehat{f},\widehat{z})=(1,0)$
$(\widehat{x},\widehat{y},\widehat{z})=(0,1,0)$	$(\widehat{f},\widehat{z})=(0,0)$
$(\widehat{x},\widehat{y},\widehat{z})=(0,0,1)$	$(\widehat{f},\widehat{z})=(0,1)$
$(\widehat{x},\widehat{y},\widehat{z})=(\frac{\sigma }{\alpha }(1-\widehat{z}),(1-\frac{\sigma }{\alpha })(1-\widehat{z}),\widehat{z})$	$(\widehat{f},\widehat{z})=(\frac{\sigma }{\alpha },\widehat{z})$
$(\widehat{x}=f,\widehat{y}=1-f,\widehat{z}=0)$	$(\widehat{f},\widehat{z})=(\beta ,0)$

.

The first three equilibrium points in Table 1 are the equilibrium points corresponding to the vertex of ${\mathcal{S}}_3$, the fourth in the interior of the simplex ${\mathcal{S}}_3$, and the fifth in the border $xy$. In the next subsections, the effect of d in the equilibrium point will be analyzed. For that reason, $(\widehat{f},\widehat{z})$ will henceforth be redefined as $({\widehat{f}}_d,{\widehat{z}}_d)$ (for $d\ne 0$) to distinguish from the equilibrium when $d=0$, which will henceforth be defined as $({\widehat{f}}_0,{\widehat{z}}_0)$.

Before further analysis, let us define another threshold for the parameter d that will be used in the following sections. Notice that Expression (18) when $z=0$ is reduced to $\left(1-r/n\right)-\sigma \left(d/1-d\right)$, from which $d_2=(n-r)/(n\sigma +n-r)$ can be obtained.

3.2.1. Effect of d over the Value of z in the Equilibrium

Here, the effect of parameter d on functions $g(f,z,d)$, $\tilde{g}(z,d)$, and in $\widehat{z}$ itself is analyzed. To this end, consider that, for $d_0=0$, there is a value of ${\widehat{z}}_0\in (0,1)$ such that $\tilde{g}({\widehat{z}}_0,0)=m\left({\widehat{z}}_0\right)=0$ (see Expressions (10) and (18)). We claim that the punishment parameter d, when $d>d_0=0$, displaces the value ${\widehat{z}}_0$ to another ${\widehat{z}}_d$ such that ${\widehat{z}}_d\le {\widehat{z}}_0$ with $\tilde{g}({\widehat{z}}_d,d)=0$. To corroborate this affirmation about the effect of d, we introduce the following two lemmas. In all cases, Assumptions 1 and 2 are fulfilled. The prime symbol denotes the derivative of the function with respect to its independent variable, for instance, $m^{\prime }\left(z\right)=\frac{dm\left(z\right)}{dz}$. We are interested in the relationships between $\tilde{g}(z,d)$ and $\tilde{g}(z,0)=m\left(z\right)$ and between ${\tilde{g}}^{\prime }(z,d)$ and $m^{\prime }\left(z\right)$ for a small $d>0$ and for a z in the neighborhood of ${\widehat{z}}_0$ where $z>{\widehat{z}}_0$.

Lemma 2.

Consider ${\widehat{z}}_0\in (0,1)$ associated to $d_0=0$ such that $\tilde{g}({\widehat{z}}_0,0)=m\left({\widehat{z}}_0\right)=0$. Let $0<ϵ<<1$ such that $\mu ={\widehat{z}}_0+ϵ$, and $m\left(\mu \right)$, $\tilde{g}(\mu ,d)$, $m^{\prime }\left(\mu \right)$, and ${\tilde{g}}^{\prime }(\mu ,d)$ are negative. Thus, a) $|\tilde{g}(\mu ,d)|\ge |m\left(\mu \right)|$ and b) $|m^{\prime }\left(\mu \right)|\ge |{\tilde{g}}^{\prime }(\mu ,d)|$.

The effect of the parameter d on the solution ${\widehat{z}}_d$ of $\tilde{g}(z,d)=0$ with respect to the solution ${\widehat{z}}_0$ of $\tilde{g}({\widehat{z}}_0,d_0)=m\left({\widehat{z}}_0\right)=0$ of the free system model is analyzed in the following lemma (see the proof in Appendix A.3).

Lemma 3.

Considering the values of ${\widehat{z}}_0$ and ${\widehat{z}}_d$ associated to $d_0=0$ and $d>d_0$ such that $\tilde{g}({\widehat{z}}_0,0)=m\left({\widehat{z}}_0\right)=0$ and $\tilde{g}({\widehat{z}}_d,d)=0$; hence, ${\widehat{z}}_d\le {\widehat{z}}_0$.

It is important to mention that, due to the value of d considered being small and the smoothness of $m\left(z\right)$ and $\tilde{g}(z,d)$, the high order terms in the expansion in (A4) do not significantly affect the analysis and can be neglected. The results of Lemmas 2 and 3 can be observed in Figure 3, where the functions $m\left(z\right)$ and $\tilde{g}(z,d)$ for a $d=0.01$ with respect to z in the interval $(0,1)$ are presented. The effect of displacing $\widehat{z}$ to the left is observed; i.e., ${\widehat{z}}_d$ with $d\ne 0$ is smaller than $\widehat{z}$ with $d=0$.

Figure 4 shows the evolution of ${\widehat{z}}_d$ with respect to several values of the parameter d (increasing values) remaining constant n, r, and $\sigma$. It can be observed that the main effect of the parameter d is precisely to displace $\widehat{z}$ to the left from its initial position. This property is exploited for the remainder of the article.

3.2.2. Effect of d over the Value of f in the Equilibrium

Here, the effect of the parameter d on the value of f in the equilibrium is analyzed. Recall that, at the interior equilibrium point, $\widehat{f_d}=\frac{\sigma }{(r-1)+d(\sigma -(r-1\left)\right)}$. Furthermore, due to game setup, $\sigma <(r-1)$ (see [18]). As a consequence, we have the expression

$$(r-1)+d(\sigma -(r-1\left)\right)<(r-1)$$

(22)

for all $0<d<1$. This means that, for a $0\le d<\tilde{d}$ and its corresponding $\widehat{f_0}$, $\widehat{f_d}$, and ${\widehat{f}}_{\tilde{d}}$, respectively, the following relationship holds:

$$\widehat{f_0}=\frac{\sigma }{(r-1)}\le \widehat{f_d}=\frac{\sigma }{(r-1)+d(\sigma -(r-1\left)\right)}<{\widehat{f}}_{\tilde{d}}=\frac{\sigma }{(r-1)+\tilde{d}(\sigma -(r-1))}.$$

(23)

Therefore, the parameter d increments the equilibrium value $\widehat{f}$ or, equivalently, improves the level of cooperation. Figure 5 shows the changes in the value of $\widehat{f}$ and $\widehat{z}$ when d increases its value from 0 to $d_2$. It can be seen how $\widehat{f}$ increases while $\widehat{z}$ decreases.

3.2.3. Effect of d in the Interior Equilibrium Point

Using the definition for a (introduced in (7)), Equation (14) (or equivalently Equation (15)) can be rewritten as follows:

$$\left\{\begin{array}{ccc}\dot{f}& =& -f(1-f)g(f,z,d)\hfill \\ \dot{z}& =& z(1-z)\left(\left(1-z^{n-1}\right)\left(\sigma -f(r-1)\right)-dfr\left(1-a\right)\left(f-1\right)\right).\hfill \end{array}\right.$$

(24)

For simplicity, we will drop the arguments in $g(f,z,d)$ and write simply g. To characterize the equilibrium point at the interior of the simplex, we consider the Jacobian of System (24), which has the following form:

$$J=\left(\begin{array}{cc}{J}_{11}& J_{12}\\ J_{21}& J_{22}\end{array}\right)$$

where

$$J_{11}=d\dot{f}/df=-((1-2f)g-f(1-f)dr(1-n)),$$

(25)

$$J_{12}=d\dot{f}/dz=-f(1-f)((n-1)(r-1)z^{n-2}-r{n}^{\prime }+drf{n}^{\prime }),$$

(26)

$$J_{21}=d\dot{z}/df=z(1-z)((1-r)(1-z^{n-1})-(1-a)dr(2f-1))$$

(27)

$$J_{22}=d\dot{z}/dz=(1-2z)((1-z^{n-1})b-e(1-a))+z(1-z)(-(n-1)z^{n-2}b-e\left(-a^{\prime }\right))$$

(28)

with $a^{\prime }=\left(1-z^n-n{z}^{n-1}(1-z)\right)/\left(n{(1-z)}^2\right)$, $b=\sigma -f(r-1)$ and $e=drf(f-1)$.

In the interior equilibrium equilibrium point, $g=0$ and $(1-z^{n-1})b-e(1-a)=0$. Evaluating the Jacobian then yields

$$\begin{array}{c}\hfill J=\left(\begin{array}{cc}f(1-f)dr(1-a)& -f(1-f)((n-1)(r-1)z^{n-2}\\ & -r{a}^{\prime }+drf{a}^{\prime })\\ & \\ z(1-z)((1-r)(1-z^{n-1})& z(1-z)(-(n-1)z^{n-2}b\\ -(1-a)(2drf-dr))& -e\left(-a^{\prime }\right))\end{array}\right).\end{array}$$

(29)

To simplify the analysis, considering the influence of the parameter d, we separate the Jacobian in two matrices $J_R$ and $J_T$ of dimension $2\times 2$ as $J=J_R+d{J}_T$ with

$$\begin{array}{c}\hfill J_R=\left(\begin{array}{cc}0& -f(1-f)\left((n-1)(r-1)z^{n-2}-r{a}^{\prime }\right)\\ & \\ z(1-z)((1-r)(1-z^{n-1})& -z(1-z)(n-1)z^{n-2}(\sigma -f(r-1))\end{array}\right)\end{array}$$

(30)

and

$$\begin{array}{c}\hfill J_T=\left(\begin{array}{cc}f(1-f)r(1-a))& -f(1-f)\left(rf{a}^{\prime }\right)\\ & \\ -z(1-z)\left((1-a)r(2f-1)\right)& z(1-z)rf(1-a)(f-1)a^{\prime }\end{array}\right).\end{array}$$

(31)

It is important to mention that all the entries of the matrices $J_R$ and $J_T$ have bounded derivatives. If $d=0$, the Jacobian matrix is reduced to $J_R$. When (30) is evaluated at ($\widehat{f_0},\widehat{z_0}$), it becomes

$$\begin{array}{c}\hfill J_R=\left(\begin{array}{cc}0& J_{R_{12}}\\ & \\ J_{R_{21}}& 0\end{array}\right).\end{array}$$

(32)

Since $\widehat{f_0}=\sigma /(r-1)$ (see (21)), $J_{R_{22}}=0$. The term $J_{R_{21}}$ is negative because $r>1$. In addition, the term $J_{R_{12}}$ is positive, for ${\widehat{z}}_0\in (0,1)$ and for the parameters n, $\sigma$, and r with values normally used in the literature. Figure 6a shows the value of $J_{R_{12}}$ for $0\le z_0\le 1$ and particularly the value in the equilibrium (${\widehat{f}}_0,{\widehat{z}}_0$).

The eigenvalues of $J_R$ are given by ${\lambda }^2=J_{R_{21}}{J}_{R_{12}}$, observing that $J_{R_{21}}$ is negative and $J_{R_{12}}$ is positive. The eigenvalues of $J_R$ are then pure imaginary. This result is important for the rest of the analysis of the interior equilibrium point when $d\ne 0$. The fact of the eigenvalues of J being pure imaginary is not conclusive when the interior equilibrium point is a center. However, we are assuming here that $d=0$, so the system is reduced to the model in [18], which proves, using a Hamiltonian technique, that the interior equilibrium point is a center.

Next, we analyze the case of the Jacobian matrix J perturbed by the parameter d. To this end, we consider that $d\approx 0$. It is important to remark the implicit dependence of the matrices $J_R$ and $J_T$ on the parameter d.

Theorem 1.

Suppose a Jacobian matrix of System (24) denoted by J and decomposed as $J=J_R+d{J}_T$ with entries $J_{R_{ij}}$ and $J_{T_{ij}}$ for $i,j\in \{1,2\}$ with bounded derivatives. Suppose also an initial equilibrium point $({\widehat{f}}_0,{\widehat{z}}_0)$ of System (24) that corresponds to the parameter $d=0$, and another equilibrium point for $0<d$ denoted by $({\widehat{f}}_d,{\widehat{z}}_d)$. Thus, the Jacobian J, evaluated at $({\widehat{f}}_d,{\widehat{z}}_d)$ with $d\approx 0$, has complex eigenvalues with positive real parts.

Theorem 1 (see the proof in Appendix A.4) shows that, for a small parameter d (considered positive), the eigenvalues of the Jacobian J of System (24) are complex with positive real parts when evaluated at the interior equilibrium point $({\widehat{f}}_d,{\widehat{z}}_d)$. This means that this equilibrium point is an unstable focus; i.e., all solution orbits of System (24) are repelled from the equilibrium point to the boundaries. This can be observed in

Table 2. Effect of the parameter d on the entries of the Jacobian J evaluated at equilibrium point $({\widehat{f}}_d,{\widehat{z}}_d)$ and its corresponding eigenvalues. For $d=0$, $({\widehat{f}}_d,{\widehat{z}}_d)$ = $({\widehat{f}}_0,{\widehat{z}}_0)$. Parameters $N=5$, $r=3$, and $\sigma =1$; $\psi =J_{R_{22}}-{\epsilon }_{22}+{\epsilon }_{11}$ and $\Delta ={(J_{R_{22}}-{\epsilon }_{22}+{\epsilon }_{11})}^2-4(J_{R_{21}}\pm {\epsilon }_{21})(J_{R_{12}}+{\epsilon }_{12})$.

d	$\mathit{J}={\mathit{J}}_{\mathit{R}}+\mathit{d}{\mathit{J}}_{\mathit{T}}$	$({\widehat{\mathit{f}}}_{\mathit{d}},{\widehat{\mathit{z}}}_{\mathit{d}})$	$\mathit{\lambda }=0.5\left(\mathit{\psi }\right)\pm 0.5\sqrt{\mathit{\Delta }}$
$d=0.01$	$\left[\begin{array}{cc}0& 0.24831\\ -0.47471& 0.00046\end{array}\right]+d\left[\begin{array}{cc}0.48162& -0.21732\\ -0.00437& -0.06878\end{array}\right]$	$\begin{array}{c}(0.5025,0.4516)\end{array}$	$\begin{array}{c}{\lambda }_1=0.00210-0.34182i\hfill \\ {\lambda }_2=0.00210+0.34182i\hfill \end{array}$
$d=0.001$	$\left[\begin{array}{cc}0& 0.24688\\ -0.47454& 0.00005\end{array}\right]+d\left[\begin{array}{cc}0.47778& -0.22112\\ -0.00044& -0.06996\end{array}\right]$	$\begin{array}{c}(0.5003,0.4604)\end{array}$	$\begin{array}{c}{\lambda }_1=0.00021-0.34213i\hfill \\ {\lambda }_2=0.00021+0.34213i\hfill \end{array}$
$d=0$	$\left[\begin{array}{cc}0& 0.24671\\ -0.47450& 0\end{array}\right]+d\left[\begin{array}{cc}0.47735& -0.22155\\ 0& -0.07008\end{array}\right]$	$\begin{array}{c}(0.5,0.4613)\end{array}$	$\begin{array}{c}{\lambda }_1=-0.34215i\hfill \\ {\lambda }_2=+0.34215i\hfill \end{array}$

, where some values of d, the corresponding equilibrium points $({\widehat{f}}_d,{\widehat{z}}_d)$, the entries of the Jacobian J, and its eigenvalues are shown. For $d=0$, the eigenvalues of J are pure imaginary corresponding to the equilibrium point being a center. This result is consistent with the literature [18]. With small positive values of the parameter d ($d=0.001$ and $d=0.01$), the eigenvalues have positive real parts, implying that the corresponding equilibrium point is an unstable focus in accordance with Theorem 1 (see also Figure 7a). It is important to mention that, for larger values of the parameter d ( $0\le d<d_2$), the equilibrium point remains an unstable focus. Moreover, it is observed that $J_{R_{22}}$ is positive for $0<d$ and $J_{T_{22}}<J_{T_{11}}$, as mentioned in Theorem 1.

As d increases, the equilibrium point $({\widehat{f}}_d,{\widehat{z}}_d)$ changes its position from the equilibrium point $(\widehat{f_0},\widehat{z_0})$ with $d=0$ in the interior of the simplex towards the border $xy$ by decreasing ${\widehat{z}}_d$ and increasing ${\widehat{f}}_d$ (see Section 3.2.1 and Section 3.2.2).

3.3. The Equilibrium Point $(x,y,z)=(1,0,0)$

This equilibrium represents a homogeneous population composed exclusively of cooperators. When $d=0$, this equilibrium point is a saddle point [19]. Now, we characterize this equilibrium point for $0<d$ and particularly for $d_1<d$. For simplicity, the system $(f,z)$ (see (14)) will be used. As can be seen in Table 1, the equilibrium point $(\widehat{x},\widehat{y},\widehat{z})=(1,0,0)$ corresponds to $(\widehat{f},\widehat{z})=(1,0)$. Evaluating the Jacobian at this equilibrium point is obtained:

$$J=\left[\begin{array}{cc}1-r/n-dr\left(1-1/n\right)& 0\\ 0& \sigma -(r-1)\end{array}\right].$$

(33)

The eigenvalues of (33) are as follows: ${\lambda }_1=1-r/n-dr\left(1-1/n\right)$ and ${\lambda }_2=(\sigma -(r-1))$. Notice that ${\lambda }_2$ is negative given the condition in the game that $0<\sigma <(r-1)$ (Assumption 2) and the choice of d determines the sign of ${\lambda }_1$. Recall that $d_1=(n-r)/r(n-1)$. Thus, for $d<d_1$, ${\lambda }_1$ is positive, and the equilibrium point is unstable (a saddle); however, if $d_1<d$, ${\lambda }_1$ is negative and the equilibrium point is locally asymptotically stable.

To better understand this equilibrium point and its basin of attraction, let us consider the relative-entropy function used in [10,34] as a candidate Lyapunov function.

Theorem 2.

Suppose that $(\widehat{f},\widehat{z})=(1,0)$ and the candidate Lyapunov function $V_{(\widehat{f},0)}=\widehat{f}ln\left(\widehat{f}/f\right)$ such that $V_{(\widehat{f},0)}>0$ for all $(f,z)\ne (\widehat{f},0)$ and $V_{(\widehat{f},0)}=0$ for all $(f,z)=(\widehat{f},0)$. If

$$df>{\displaystyle \frac{n(1-z)(1+(r-1)z^{n-1})-r(1-z^n)}{n(1-z)r-r(1-z^n)}}$$

(34)

then ${\dot{V}}_{(\widehat{f},0)}<0$, and the equilibrium point is locally asymptotically stable (see proof in Appendix A.5).

In Figure 8, the condition (34) is analyzed for increasing values of d, with the restriction $d>d_1$, since this is a requisite for the equilibrium point to be stable. There are two main regions in the simplex. The region colored in orange contains all the states where ${\dot{V}}_{\widehat{f}}<0$, while the region colored in purple contains all the states where ${\dot{V}}_{\widehat{f}}>0$. Both regions are separated by a red line that divides the simplex from a point in the border $xy$ (where $z=0$) defined as $p_1$ to another in the border $yz$ (where $x=0$) defined as $p_2$ (see Figure 8). This division corresponds to the set of points where ${\dot{V}}_{\widehat{f}}$ changes sign.

As d increases the region where ${\dot{V}}_{\widehat{f}}<0$ also increases. This happens because the point $p_1$ moves to the vertex $y=1$. Notice that the value of $p_1$ obtained from ${\dot{V}}_{\widehat{f}}=0$ is $x=(n-r)/r(n-1)d$. This is the position of the unstable equilibrium in the border $xy$ analyzed in Section 3.1.1. The point $p_2$, however, does not depend on d. When $x=0$, the expression is reduced to $1+(r-1)z^{n-1}-ra=0$ (see (10)). Therefore, the value of z in the point $p_2$ is fixed and corresponds to the value of $\widehat{z}$ when $d=0$ as in [18].

Additionally, in Figure 8, some paths with initial value points from both regions—${\dot{V}}_{\widehat{f}}<0$ and ${\dot{V}}_{\widehat{f}}>0$—are shown. Due to the structure of the solution, even trajectories with initial values outside the basin of attraction are drawn inside and eventually end in the equilibrium point. In such a case, the equilibrium point $(1,0,0)$ is globally asymptotically stable.

Notice in Figure 8a that some trajectories oscillate before reaching the equilibrium. This is caused by the unstable equilibrium point in the interior of the simplex analyzed in Section 3.2. Recall that the interior equilibrium shifts to the border $xy$ as d increases, reaching it when $d=d_2$, which, for the parameter in Figure 8, is $d_2=0.2857$.

From a practical point of view, we are interested in trajectories where the frequency of cooperators does not decrease and, as a consequence, in the set of initial values from which the number of cooperators increases along all the trajectory to the equilibrium point. To observe this behavior, two regions are defined in the simplex ${\mathcal{S}}_3$; we denote by $\mathsf{\Omega }$ a set of initial values where the frequency of the cooperator strategy in the solution orbits continuously grows until full cooperation is reached, i.e.,

$$\mathsf{\Omega }:=\{(x\left(t_0\right),y\left(t_0\right),z\left(t_0\right)):x\left(\tilde{t}\right)\le x\left(t\right)\forall t_0\le \tilde{t}\le t.$$

(35)

In Figure 9, the region $\mathsf{\Omega }$ is in green. It can also be observed (and this is an important consideration) that $\mathsf{\Omega }$ is enlarged as the parameter d increases. In addition, the set of initial values in the simplex ${\mathcal{S}}_3$, where the frequency of cooperators decreases before reaching the equilibrium point, i.e., ${\mathcal{S}}_3-\mathsf{\Omega }$, is colored in blue.

3.4. Summary of the Fractional Punishment Effect in the System

In this section, we complete the outline presented in Figure 1 in Section 3. Notice that the whole effect of the parameter d in the dynamic and the equilibrium points of System (2) can be described using the thresholds $d_1$ and $d_2$ as introduced in Section 3.1.1 (Lemma 1) and Section 3.2, respectively.

Recall that, when $d=0$, the interior equilibrium point $q_1$ is a center, and the vertices are saddle points [18]. Figure 1a shows the periodic orbits around the equilibrium.

Positive values of d cause two transformations in the system. First, the equilibrium point $q_1$ becomes unstable; second, the position of the equilibrium changes. It moves towards the border $xy$. Figure 1b shows how, with small values of d, $d<d_1$, the trajectories approaching the border of the simplex oscillate.

As d increases, i.e., $d_1<d<d_2$, a new equilibrium point $q_2$ arises in the border $xy$. This is an unstable equilibrium point that moves to the vertex y as d increments its value. Furthermore, when $q_2$ appears, the vertex x changes its nature from unstable (saddle point) to stable (see Figure 1c,d). Notice that the number of oscillations needed to reach the vertex x decreases as d grows.

When $d=d_2$, the equilibrium point $q_1$ reaches the border $xy$ and merges with the equilibrium point $q_2$ in the border. This new equilibrium $q_3$ remains unstable (see Figure 1e).

From $d_2<d$ to $d=1$, there is no equilibrium in the interior of the simplex. The point $q_3$ continues to move toward the vertex y, but it does not reach the vertex; when $d=1$, the equilibrium point on the $xy$ border takes the value $\widehat{x}=(n-r)/r(n-1)$.

4. Discussion

In this work, we present a mathematical model in which a fractional punishment is incorporated in an optional public good game under the replicator dynamic. To design the model, we consider two common situations regarding the administration of institutions: (1) not all those who did not pay their fees can be punished due to a lack of resources; (2) both the service provided and the cost of the sanctioning system come from the fee paid by the users (fee pool).

Our proposal seeks to handle this situation by limiting the sanction to a random set of defectors, warranting at the same time both the economic feasibility and effective punishment to at least a reduced number of free riders.

To this end, we assume in the game that a randomly selected portion of free riders would be penalized, but we do not include a separate pool to be used to sanction defectors. With these assumptions, we eliminate the possibility of second-order free-riding and the cost of second-order punishment, since cooperators, when paying the fee, are also unintentionally contributing to penalizing free riders. Furthermore, the cost of the sanctioning system is mainly dependent on the size of the portion punished, and it may be adapted in accordance with the resources available to the institution.

Small values of d produce solutions that oscillate approaching the border. The effect of the parameter d on the trajectory towards the equilibrium point is an important issue. A wide excursion could jeopardize the survival of any institution because of the increment of the loners and the possibility that this triggers the disappearance of the active population.

As d increases, the participation in the game is improved by reducing the share of loners within the population. This can be understood by the convergence of the equilibrium point to the $xy$ axis. For an appropriate minimum value of d, a population composed exclusively of cooperators can be achieved. If d surpasses a first threshold value $d_1$, a new unstable equilibrium point arises in the border $xy$ while the equilibrium point in the vertex x becomes stable. Then, when d achieves a second threshold value $d_2$, the interior equilibrium point reaches the border, and both equilibria merge. Notice that, with fractional punishment, the equilibrium reaches the border with $r<n$, whereas in the model introduced in [18], it occurs when $r=n$.

Even with a value of d adequate to achieve full cooperation, the trajectories can oscillate depending on the initial composition of the population. In this context, the set of initial values from which full cooperation is attained through trajectories with an ever-increasing frequency of cooperators is of special interest. This set of points (denoted by $\mathsf{\Omega }$ in the results) increases as d increases.

The fractional punishment is a useful mechanism when resources are scarce; it can improve cooperation at a lower cost. However, with this approach, the pressure on the correct use of the contribution (fee pool) is higher, reaching a balance between the expenditure dedicated to the operation with the expenditure on sanction free-riding is more important. Sanctioning free riders improves cooperation but decreases the profit from the game. We consider that an adequate balance between both expenses is perhaps one of the characteristics of a successful institution.

5. Conclusions

We have analyzed the effect of randomly punishing a percentage of defectors in an optional public good game modeled using replicator dynamics. To this end, the payoff of the strategies is redefined to introduce a punishment parameter. This parametrizes the simplex phase diagram, modifying the trajectories of the strategies in the population as well as the nature of some important equilibrium points. A description of the main characteristics variations of the trajectories in the simplex due to the punishment parameter parametrization was presented, and for the equilibrium point that represents full cooperation, the basin of attraction was described. Future works include modeling how the institution decides to set the value of the parameter d (probably adaptive) and its optimal choice, as well as the redistribution to cooperators of the benefits taken from the punished free riders.