Coexistence of coordination and anticoordination in nonlinear public goods game

There is a plethora of instances of interactions between players, in both biological and socio-economical context, that can be modeled as the paradigmatic PGG. However, in such interactions, arguably the PGG is often nonlinear in nature. This is because the increment in benefit generated, owing to additional cost contributed by the players, is realistically seldom linear. Furthermore, sometimes a social good is created due to interspecific interactions, e.g. in cooperative hunting by animals of two different species. In this paper, we study the evolutionary dynamics of a heterogenous population of cooperators and defectors playing nonlinear PGG; here we define heterogenous population as the one composed of distinct subpopulations with interactions among them. We employ the replicator equations for this investigation, and present the non-trivial effects of nonlinearity and size of the groups involved in the game. We report the possibility of discoordination, and coexistence of coordination and anti-coordination in such nonlinear PGG.


Introduction
In order to sustain a public good (also referred as common), either all the participating agents coordinate to cooperate or some agents voluntarily cooperate while the rest anticoordinate to defect; of course, the former is a more desirable and efficient scenario, but realistically the latter scenario is frequently encountered as revealed by the presence of free riders. Interestingly, the agents under question need not exclusively be Homo Economicus: they can be rather cognitively inferior biological entities-like microbes, plants, and animals-who play intriguing roles in the tragedy of the commons [1][2][3] encountered in evolutionary biology [4].
The traditional mathematical setup used to analyze the fate of a common in both classical game theory and evolutionary game theory is the paradigmatic public goods game (PGG) [5]. In PGG, a cooperating individual contributes a cost towards creation of a public good, while a defecting individual (free rider) does not contribute anything. In the simplest form of the PGG-the so-called linear PGG [6][7][8][9]-the benefit received by each player is decided by a constant real number termed synergy factor: the total contribution by the cooperators is multiplied by a synergy factor and the resultant is equally distributed among all the players.
The linear public good game and prisoner's dilemma [10]-its N-person version to be precise [11]-are known to be equivalent [12]. Considering two person prisoner's dilemma, for simplicity, one may recall that in prisoner's dilemma defection is the only symmetric Nash equilibrium (and evolutionary stable strategy). It is interesting to note that although other well-studied dilemmas like the anti-coordination dilemmas (such as the snowdrift game [13,14]) or the coordination dilemmas (e.g. the stag-hunt game [15,16]) do not appear in a linear PGG, at least (N-person) stag-hunt dilemma does appear when there exists a requirement of a finite threshold [8] of cooperators to create the public good. A stable anti-coordination strategy profile is not realizable even with such a threshold in the linear PGG.
However, it must be borne in mind that the linear payoff function in the linear PGG is a merely a crude approximation of the actual nonlinear payoff functions often observed in biological populations [17],

Nonlinear PGG in heterogenous population
Let us consider a very large, mathematically infinite, well-mixed heterogeneous population consisting of two subpopulations tagged A and B; each subpopulation, in turn, exclusively consists of cooperators and defectors. The frequencies of the cooperators in subpopulation A and B is denoted by x A and x B respectively, and the remaining frequencies-viz., (1 − x A ) and (1 − x B ) respectively-is obviously that of the defectors in subpopulations A and B. From subpopulation A, N A and from the other, N B number of players are chosen randomly who play a PGG collectively. We denote the (random) number of cooperators in the group of N A players as k A and and in the group of N B players as k B . Each cooperator of type A contributes a cost c A > 0, and similarly each cooperator of type B contributes a cost c B > 0 towards the PGG. The defector are free riders in both the subpopulations.
As motivated in the introduction of this paper, the accumulated contribution from each subpopulation may not increase linearly. Thus, we consider for generality a nonlinear payoff function [27]-with an extra parameter, β ∈ R + is addition to the usual synergy factor, α ∈ R + -that can be concave or convex depending on the parameters. We call β synergy exponent in this paper. Assuming that the final collective benefit is shared equally among all individuals, the payoff of a defector in either population is given by whereas the payoffs of individual cooperators of type A and type B are respectively given by Note that β > 1 means that every additional cost is more prized while β < 1 means just the opposite; and β = 1 leads to the usual linear PGG. Of course, a more general setup demands that the synergy factor (α) and the exponent (β) can be different for different subpopulations. Nevertheless, in this paper, for the sake of simplicity and reducing the number of parameters, we refrain from taking that course. Now in the context of the evolutionary game theory, fitness is conventionally defined as the expected payoff. One can think of x A and x B as the probabilities of a randomly chosen player being a cooperator of respectively type A and B. Hence, assuming that the groups playing the PGG are formed randomly in accordance with binomial sampling, the probability that a focal A type cooperator should find k A − 1 number of A type cooperators in a group of size Moreover, there may be contributions from k B ∈ {0, 1, . . . , N B } cooperators of B type: The probability that the focal A type cooperator receives benefit from k B number of B type cooperators in a randomly formed group of size , of an A type cooperator in a heterogenous population where nonlinear PGG is being played is given by the average payoff: Likewise the fitness, f D A (x A , x B ), of an A type defector, who can get share of the contributions from at most (N A − 1) + (N B ) cooperators in the heterogenous populations, is given by The expressions for the fitness, , of a B type cooperator and defector respectively is obtained by simply replacing the subscripts A by B and vice versa in the terms (except in the arguments of the payoffs) present in the right hand side of equations (4) and (5). We observe that the fitnesses of a particular type individual not only depends on the frequency of the cooperators in its own subpopulation but also on the cooperator frequency in the other subpopulation.
Finally, in order to follow the evolutionary dynamics of the well-mixed heterogenous population, we adopt the well-known replicator dynamics [28,41,42] for two strategy and two subpopulation to write: The above continuous-time equations implicitly assume that we are considering the population to be generation-wise overlapping [43,44], otherwise discrete-time version of the equations would be more appropriate. We expect [45] the dynamic equilibria of the replicator equations to correspond to the game-theoretic concepts of Nash equilibrium and evolutionary stable strategy. This is precisely what we now want to explore: Specifically, we want to investigate which stable states of the population correspond to which kind of frequency distribution in the population, and when the robust collective coordination and anticoordination appear in the nonlinear PGG.

Analyzing the replicator model
Before we analyze equations (6a) and (6b) with equations (4) and (5) (and the analogous expressions for subpopulation B) in mind, let us set N A = N B = N, without too much loss of generality as far as qualitative results are concerned. This makes the presentation of analysis of the model uncluttered and conspicuous.
As it is known that full analytic solution of the replicator equations is not possible, so we resort to the fixed point analysis and numerical investigation of the equations. While technically the dynamics is happening in the Cartesian product space of two 1-simplices, the phase space dynamics under consideration through the replicator equation is effectively a closed unit square, x A -x B space. Given that it is two-dimensional replicator system, we do not expect any limit cycle or chaotic orbits in this system. Consequently, we exclusively focus on the convergent solutions that corresponds to the fixed points.

Existence of fixed points
The existence of the four corner fixed points-(0, 0), (0, 1), (1, 0), and (1, 1)-for all values of the parameters (like α, β, c A , c B , and N) is rather obvious. There in principle may be many more fixed points depending on the simultaneous roots of the two (2N + 1)-order polynomials of two variables: the right hand sides of equations (6a) and (6b). However, out of the many possible fixed points, we are only interested in the ones that lie on the unit square because only those have physical meanings. Following considerations are helpful in finding the interior and the boundary fixed points.

Proposition. There is at most one interior fixed point in the system of equations (6a) and (6b) when
In order to prove this proposition for the case β ̸ = 1, we first need to know the four facts written below: Let us discuss only the proof for fact (a) written above (as the next three facts can be proven analogously). To this end, we take the partial derivatives of Q A (x A , x B ) with respect to its arguments and simplify to arrive at the following: Here for β < 1 (actually for any β). In equation (7), it is apparent that all the terms-except ξ(k)-are always positive: Determination of the sign of ξ(k) requires more consideration which we now undertake. We note that where the integrand g(y) ≡ y β−1 is monotonically decreasing if β < 1 and k ∈ N. Therefore, or in other words, B(k) < B(k − 1) ∀k ∈ N. This implies that ξ(k) is negative ∀k ∈ N. Consequently, the right hand side of equation (7) is negative which means Q A (x A , x B ) is monotonically decreasing in x A direction when β < 1. Hence, fact (a) mentioned above is proven completely and facts (b)-(d) can be proven along the similar line of argument.
Now, coming to the proof of the aforementioned proposition for β ̸ = 1, we note that the internal fixed point of the system of equations (6a) and (6b) are formed by the simultaneous roots of the polynomial equations, Q A (x A , x B ) = 0 and Q B (x A , x B ) = 0. We have already argued that Q A monotonically increases along x A and Q B monotonically increases along x B . It helps to visualize the graphical representation of the functions, Q A and Q B , in a three-dimensional plot where the third axis (perpendicular to both x A and x B axes) represents the value of the functions: The functions, Q A (x A , x B ) and Q B (x A , x B ), are generically two-dimensional surfaces. Owing to the specific monotonicity properties, these two surfaces must intersect along a one dimensional curve. Obviously, if at all this curve intersects x A -x B plane, then it can do so at one and only one point which would correspond to the interior fixed point. Hence the proposition stands proven for the case β ̸ = 1.
The case of the unity synergy exponent is rather trivial. In this case the derivative of Q A (x A , x B ) is zero. Therefore it is neither monotonically increasing nor decreasing, it is a constant. So, it cannot grow from positive to negative, hence there is no internal fixed point. The same argument can be given for Hence the proposition stands proven for the case β = 1 as well.
Whether the interior fixed point (x * A , x * B ) exists or not, depends on the specific values of the system parameters, viz., α, β, c A , c B , and N. For the case β < 1, the existence of this unique interior fixed point, the parameters should be such that Whereas for the case β > 1, the existence condition is given as, The rationale behind the above set of inequalities is simply that, being monotonic, Q A and Q B must grow from a negative value to a positive value or from a positive value to a negative value in order to become zero on the phase space.
Lastly, there can be four more fixed points, viz., the boundary fixed points: Consider the case β < 1: using the above arguments, a nonzero x * A exists when conditions equations (11a) and (11b) are satisfied irrespective of whether conditions equations (11c) and (11d) hold. This means that for β < 1, (x * A , 0) and (x * A , 1) exist whenever conditions equations (11a) and (11b) are satisfied. Similarly, for β < 1, (0, x * B ), and (1, x * B ) exist when conditions equations (11c) and (11d) hold irrespective of whether conditions equations (11a) and (11b) are satisfied. It may be noted that it trivially follows from this argument that interior fixed point always appears along with all the four boundary fixed points.
For completeness, we mention that for the case β > 1, In this case also, the interior fixed point always appears along with all the four boundary fixed points. Finally, it is easy to conclude that when the synergy exponent is unity, the boundary fixed points do not exist.
Summarizing, depending on the parameters' values, our system has a minimum of four (corner) fixed points co-existing and has a maximum of nine fixed points existing simultaneously; only other possibility is coexistence of six fixed points.

Stability of fixed points
By knowing the existence and the stability of the fixed points, we can get the complete flow diagram for the system for a given set of parameter values. Having already discussed the existence of the fixed points, we now proceed to investigate their stability properties. To this end, we calculate the eigenvalues of the Jacobian matrix of the corresponding linearized equations.
We find that all the eigenvalues are real (meaning that a fixed point is either a node or a saddle point) and their signs are determined by Q 0 , However, owing to monotonicity of can equivalently be found from the signs of Q 0 (c A ) and Q 1 (c A ), and Q 0 (c B ) and Q 1 (c B ) respectively. Therefore, it is paramount to know how, for a given α, β, and N, the signs of Q 0 and Q 1 change with c (which refers to c A or c B as the case may be). Form the monotonicity properties of these functions, we know that Q 0 (c) > Q 1 (c)∀c ̸ = 0 when β < 1; and Q 0 (c) > Q 1 (c)∀c ̸ = 0 when β > 1. Additionally, it is easily seen that in the former case the curvatures of Q 0 (c) and Q 1 (c) are negative, while in the latter case the curvatures are positive. Therefore, figures 1(a) and 2(a) schematically showcases the comparative plots of Q 0 (c) and Q 1 (c) for the cases β < 1 and β > 1 respectively. The non-zero roots c (0) th and c (1) th of the two functions are with the property that c th ∀β ∈ R + . Thus, the signs of Q 0 and Q 1 , and hence signs of the eigenvalues, are determined by the choice of the values of the costs relative to c (0) th and c (1) th . Since we can confine ourselves to the case of c B ⩾ c A (without any loss of generality), therefore there can be six distinct ways we can simultaneously pick c A and c B . This is simply because when c A ∈ (0, c     Based on the aforementioned analysis, in figures 1 and 2 we illustrate all possible qualitatively distinct phase portraits that can be realized in the nonlinear PGG played in a heterogeneous population.

Numerics
The analysis based on how c A and c B are placed relative to c (0) th and c (1) th is quite general in the sense that, for any set of values of α, β, and N that completely decides the values of c (0) th and c (1) th , the phase portrait can be predicted using the analysis. However, in the absence of explicit solution of the replicator equations, we must resort to the numerical investigation of the equations in order to see the comparative effects of α and β. For illustrative and discussion purpose, we present nine plots among many generated in the aforementioned fashion. In figure 3, we take three values of N (3, 30, and 300), a fixed value of c A = 10, and three values of c B (10,30,3000). Note that there are actually ten different distinct scenarios (phase portraits in figures 1 and 2) and various combinations of them may appear in α-β parameter space depending on the other parameters' values.

Results
The effect of group size in linear PGG is a topic of active experimental and theoretical research [46][47][48][49][50][51][52][53][54][55][56][57]. Here we can look at the effect of group size in nonlinear PGG in heterogeneous population from a little different perspective. The individuals of the subpopulation may choose to incur costs whose exact value naturally varies from situation to situation. Therefore, given the fact that the costs can realistically range from a small non-zero value to a larger upper value, a pertinent question is: What is the most probable outcome in a nonlinear PGG? Of course, the result depends on the exact value of N-the group size. However, what happens in the limit N → ∞ is insightful.
Consider the case synergy exponent β < 1. Recall (from equations (13) and (14)) that the two threshold parameters, c th ) becomes most probable choices of cost. Ergo, in this case as well, only mutual defection in both the subpopulation turns out to be the most probable state (see figure 2(b)).
The case of the unity synergy exponent, i.e. linear PGG, is relatively rather bland. Here Q 0 (c) = Q 1 (c) = c(α/2N − 1). This equality implies that no fixed point except the four corner ones exist. Straightforward linear stability analysis divulges that while (0, 1) and (1, 0) are always saddle points, (0, 0) and (1, 1) are nodes. The stability property of the fixed points are independent of the cost values: When α/2N > 1, only (1, 1) is stable and when α/2N < 1, only (0, 0) is stable. Thus, for high values (higher than twice the group size) of synergy factor, both the subpopulations' final states are exclusively composed of cooperators. Whereas, for low values (lower than twice the group size) of synergy factor, both the subpopulations become full of defectors.
As soon as we leave the paradigm of linear PGG, a richer set of non-trivial stable population states come into the picture. First, general point to note (see figures 1 and 2) is that while coexisting stable population states are possible only when β > 1; the case of β < 1 leads to unique final states irrespective of what the initial composition of the two subpopulations are. Thus, it appears that the observation of coexisting states in heterogenous populations playing PGG indicates that the situation may be modeled as a nonlinear PGG where the increment in benefit increases owing to contributing additional cost by the individuals in the subpopulations. Of particular interest are discoordination state (which exclusively appears for β < 1) and state with coexisting coordination and anticoordination profiles (which exclusively appear for β > 1). The former is denoted as MM and the latter as {CC, DC, CD, DD} in figures 1-3. Such states are not at all possible in linear PGG (see [8]).
In figure 3, for a given group size, we note that with increase in the cost-difference (c B − c A ), the ranges of α and β-for which the discoordination and the coexisting coordination-anticoordination states exist-progressively shrink. In fact, the regions can disappear altogether for very high cost-difference (see figure 3(b)) vs. figure 3(c)). Whereas the same ranges progressively swell with group-size N at a give cost-difference; interestingly, the increase in N can revive the regions that disappeared at very high values of cost-difference (see figure 3(f)) vs. figure 3(i)). Furthermore, the regions of discoordination and the coexisting coordination-anticoordination are bounded by full cooperation (CC) or full defection (DD) states when the cost-difference is zero, i.e, there is no asymmetry in the contributions by the cooperators of the two subpopulations. As the asymmetry sets in, the parameter-space becomes richer with the appearance of other states (viz. , CM, MD, {DC, CC}, and {DD, DC}) which separate the regions of discoordination and the coexisting coordination-anticoordination from the full cooperation (CC) and the full defection (DD) states.
The order of appearance of the such other states can be understood through the incentive to cooperate per unit cost, defined as: , (e) and (h)) for β < 1. Let us remind ourselves of the obvious fact that in any subpopulation (A or B) state D has less cooperation than that of state M and state M has less cooperation than that of state C. Thus, in the sequence under consideration, where each population state is composed of states of the two subpopulations, it is clear that the level of cooperation increases faster in (first) subpopulation A than in (second) subpopulation B as α increases. Recall that in this case c A < c B . This means that as the synergy factor increases, the subpopulation with higher cost cooperates at higher α as compared to the subpopulation with lower cost; in other words, the cooperators require to forsee more returns (through more synergy factor) in order to contribute more. This is rather obvious when one glances at equations (15a) and (15b): c B > c A implies I cB < I cA for a given α. Therefore, population B has comparatively less incentive to cooperate at a given synergy factor. Of course, for particular any subpopulation, the cooperation level increases simply because the corresponding incentive to cooperate per unit cost increases with the synergy factor.

Discussion
In this paper, we have considered a population composed of two subpopulations of individuals. In each subpopulation, every individual can either cooperate or defect while playing a nonlinear PGG. We have extensively analyzed the nonlinear PGG while considering the simplest possible power-law nonlinearity in the payoff function which is characterzied by two parameters: a synergy factor, α and a synergy exponent, β. The state of the population has been evolved using the appropriate replicator equations. Through elaborate analytical and numerical analyses, we have found all possible asymptotic states of the population; and have explored the effects of group-size (N), synergy factor, and synergy exponent. Although, during the analysis we have kept the synergy factor α same in both the subpopulations, the case of different synergy factors has been implicitly covered because one can always rescale the costs so as to make the synergy factors identical. However, what happens if the synergy exponent β and the group size are different for different subpopulations is a question whose detailed analysis is worth pursuing in future.
We have specifically shown that nonlinearity in PGG in an evolutionary system, having simultaneous presence of inter-and intra-specific interactions, is capable of effecting bistable population states. Moreover, discoordination strategy profiles, and coexistence of coordination and anticoordination strategy profiles are new phenomena that the nonlinear PGG shows-the former one results when the payoff function is concave, and the latter when the payoff function is convex. These states are quite dependent on the group size; in fact, irrespective of the finite values of the synergy factor and the synergy exponent, if very large groups are participating in the PGG, then mutual defection is the only viable state of the population.
It is worth mentioning that our analysis is more general than what is presented herein in the following sense. The continuous-time replicator equation used in this paper implicitly assumes that we are considering the population to be generation-wise overlapping [43,44], otherwise discrete-time version of the equation would be more appropriate. The small time-step limit of the discrete-time replicator equation is actually the adjusted replicator equation [58]:ẋ . (16b) The number of fixed points of equations (16a) and (16b), and their stability properties are identical to the ones for the replicator (equations (6a) and (6b)). This is because the denominators on the righthand sides of the adjusted replicator equation are positive definite and the Jacobian matrices in the linearized equations are diagonal when evaluated at the fixed points.
Before we end, we remark that for the simpler case of homogeneous population, nonlinear PGG is qualitatively different from linear PGG. To see this explicitly, let us assume there are x fraction of cooperators in a well-mixed homogeneous population, i.e. there is only one subpopulation and interspecific interaction is trivially absent. Thus, the payoff of a defector in a group of randomly chosen N players with k cooperators, is Π D (k) = (αkc) β /N. The payoff of cooperator is given by Π C (k) = Π D (k) − c. As a result the fitnesses can be found as and the replicator equation takes the simple form: In the course of analytical and numerical study of this equation, along of the line of what has been presented in the main text, we arrive at figure 4 that summarizes some interesting results succinctly written below.
As is already known in the literature, the linear PGG (β = 1) yields full cooperation in the population for α > N but full defection in the population for α < N. This means that only N-person harmony and N-person prisoner's dilemma are possible scenarios in the linear PGG. A change in the synergy exponents value can change the state of the population. For the case, β > 1, coexistence of full cooperation state and full defection state is brought forth . It may be interpreted as the emergence of coordination strategy where the all the individuals in the population either cooperate or defect together. Hence it is apt to call such a scenario N-person stag-hunt game which may appear in PGG when the payoff function is convex function of total contribution. While it is known that N-person stag-hunt game can be witnessed in linear PGG when a non-zero threshold [8] on the number of cooperators is imposed for building the PGG, the nonlinear PGG does not require any such threshold for the emergence of N-person stag-hunt game. More interestingly, when the synergy exponent is less than unity, N-person hawk-dove game arises in the sense that there is only a mixed equilibrium state of the population is expected: Both the cooperators and the defectors coexist. N-person hawk-dove game is impossible in the linear PGG even with the inclusion of a threshold. Naturally, the phenomenon of threshold when included in the nonlinear PGG in heterogeneous population will lead to far richer dynamical outcomes which we hope to investigate in future. Of course, another natural scope of extension is to relax the ideal conditions of infiniteness and unstructuredness [59] of the population.
Another possible extension would be to incorporate the effects of punishment in the PGG. The framework of binomial distribution, as used in this paper, is not feasible when one deals with three type of players, viz., cooperator, defectors, and punishers (cooperators who punish defectors) [5]. However, if the cost of punishment is allowed to change over time, i.e. the players adapt their sanctioning efforts based on the success of cooperation in the population, then one arrives at the special case of self-organized punishment [5,60]. The case of self-organized punishment can effectively be characterized by using only two strategies (cooperation and defection, although the payoffs are now implicitly time-dependent) and hence we expect that the framework used in this paper can be carried over to such cases. The evolutionary dynamics presented herein can be further extended to encompass eco-evolutionary dynamics to study the effects of punishment therein [61,62] in the presence of density dependent selection [63,64]. In this context, in line with the focus of this paper, it is interesting to point out that alternating rotation of coordinated and anticoordinated behaviors has been reported in the eco-evolutionary dynamics [65].

Data availability statement
No new data were created or analysed in this study.