Large-time asymptotics of a public goods game model with diffusion

We consider a spatially inhomogeneous public goods game model with diffusion. By utilising a generalised Hamiltonian structure of the model we study the existence of global classical solutions as well as the large time behaviour: first, the asymptotic convergence of the PDE to the corresponding ODE system is proven. This result entails also the periodic behaviour of PDE solutions in the large time limit. Secondly, a shadow system approximation is considered and the convergence of the PDE to the shadow system in the associated fast-diffusion limit is shown. Finally, the asymptotic convergence of the shadow to the ODE system is proven.


Introduction
In this work, we are interested in a PDE version of an optional public good game [3] where f and z are relative fractions of populations and we assume Here, is a bounded domain of R d with smooth boundary and outer unit normal ν.
Moreover, d f , d z > 0 are positive diffusion coefficients.
For the remaining parameters, we assume and the function G(z) is given by Note that in the parameter range (3), the function G(z) has exactly one sign change in z ∈ (0, 1) and looks qualitatively like Fig. 1, see Sect. 2 for the details. Moreover, N ≥ 3 in (3) and (4) denotes the number of players, see Sect. 2 for more details on the considered public good game [3]. Adding diffusion in these kind of models has already been considered in [1,2,8], mostly to model microbial interactions mediated by diffusible molecules, since standard game theory cannot describe such behaviour. Moreover, in contrast to human behaviour or animal colonies, microbial communities rarely rely on direct contact since microbes primarily communicate though diffusible molecules. This diffusive behaviour is the reason why such molecules are often termed public goods, see [1,2,8,12] and the references therein.
This work focuses on the study of the dynamics of PDE-problems of the type (1). Herein, the considered optional public good game [3] should be viewed as an interesting example case that leads us to study the general question of links between PDE and ODE model. In fact, we expect our mathematical analysis to similarly apply to related models to (1), which shares the below considered key properties. While global existence of classical solutions of (1) is straightforward, we are in particular interested in the asymptotic large-time behaviour of the solutions and their qualitative properties.
More precisely, a main question of this paper asks if PDE model (1), despite having sign changing terms at the right hand side of both equations, exhibits the same largetime behaviour as the corresponding ODE-model, which was originally studied in [3].
The key structural property, which will allow to characterise the large-time behaviour of the PDE model (1) is that the original ODE-model [3] features in the parameter range (3) a generalised Hamiltonian structure of the form where H 1 and H 2 are defined below in (17) and (18). The first theorem shows that PDE solutions become spatially homogeneous as t ↑ +∞ subject to (5). The proof requires the technical Lemma 9, which provides sufficient conditions to the positive definiteness of the (Hessian of the) Hamiltonian H ( f , z). Theorem 1 (Global existence and convergence to the ODE) Given f 0 , z 0 ∈ C 2 ( ) with finite Hamiltonian H ( f 0 , z 0 ) < +∞ and ∂ ∈ C 2 . Assume the Hessian of the Hamiltonian H ( f , z) to be positive definite, which holds, for instance, under the assumptions of Lemma 9.
Then, a unique global-in-time classical solution to (1) exists. Moreover, given a PDE solution ( f (·, t), z(·, t)) to system (1) Here, Remark 2 A global existence result could also be proved by the method of invariant sets (see e.g. [11]) since 0 < f 0 , z 0 < 1 implies 0 < f , z < 1 for all times. Yet, by using the Hamiltonian structure of the system, we can show (7) and get more information about the global dynamics of system (1) as stated by the following results. Moreover, the Hamiltonian approach can be extended to systems without invariant sets.
Motivated by [5], we notice that any solution to the ODE model (5) is periodic and that any PDE orbit is absorbed into one of the periodic ODE orbits O = {(f (t),z(t))} t≥0 . Thus, we derive the following consequence of Theorem 1: ,z(t))} t≥0 as defined in Theorem 1 be composed of more than one point, i.e. be a non-trivial orbit. Then, the PDE solution ( f (·, t), z(·, t)) is periodic in the large-time limit and there exists a "phase shift" λ > 0 such that Next, we consider the shadow system (see e.g. [9]) corresponding to system (1), which is (formally) obtained in the limit d z ↑ +∞: where now Z = Z (t). Shadow system (9) has homogeneous Neumann boundary conditions for F and considered subject to the initial data Shadow systems are used to approximate the parabolic problem by the equilibrium problem obtained in the limit d z ↑ +∞, see e.g [6]. Accordingly, F = F(x, t) is a space and time dependent function while Z = Z (t) depends only on t. The following theorem justifies rigorously the shadow system approximation scheme. Note that we can equally consider and prove the following results for the shadow system obtained in the limit d f ↑ +∞.
Theorem 4 (Convergence to the shadow system) Suppose the assumptions of Theorem 1 and assume u 0 , v 0 ∈ W 3,s ( ), s > d with smooth boundary ∂ . Let ( f , z) and (F, Z ) be the solutions to (1) and (9), respectively. Then, for any T > 0, holds Remark 5 Note the additional initial regularity u 0 , v 0 ∈ W 3,s ( ) constitutes the minimal regularity, which is required to prove Theorem 4. However, standard parabolic regularity implies arbitrarily regularity of solutions for arbitrarily smooth boundaries ∂ . Parabolic smoothing also implies that the initial regularity could be relaxed to u 0 , v 0 ∈ L ∞ ( ) as in (2) if statement (11) is relaxed to t ∈ [τ, T ] for τ > 0.
The asymptotics of the shadow system are also given by the ODE system as t ↑ ∞.
Theorem 6 (Convergence from the shadow to the ODE system) Under the assumptions of Theorem 1, let (F, Z ) = (F(x, t), Z (t)) be the solution to the shadow system (9).
be the solution to the ODE system (5) subject to initial datâ Discussion of the results Our paper deals with the asymptotic behaviour of solutions to the PDE model (1). We prove global existence of classical solutions to (1) and convergence to the corresponding ODE model (5). Next, we derive the interesting Corollary 3, which implies that if we have two PDE solutions of (1) that start at different times, asymptotically they will be close in the C 2 -norm modulo a suitable phase shift. After that, we prove that solutions to (1) converge to those of the corresponding shadow system obtained in the limit d z ↑ ∞ and that solutions of the shadow system converge to those of the corresponding ODE orbit. These results can also be seen as follows: If we start from the PDE system (1) but with different initial data, we will get two different solutions that have nevertheless two properties in common. First, they are both attracted from the corresponding shadow system and either they will pass close by or through the solutions of this shadow system. Second is the fact that even though we started with different initial data, asymptotically we will have convergence to the corresponding ODE in both cases as t ↑ ∞ but these two limits may be different, i.e. there could be a phase shift between them. This behaviour is similar to the Lotka-Volterra systems which was noticed in [7]. Outline This paper is organised as follows: In Sect. 2, we recall the modelling background and establish some basic properties of system (1). Theorem 1, Corollary 3, Theorems 4 and 6 are proven in Sects. 3, 4, 5 and 6, respectively.

Preliminaries: modelling and formal properties
Public goods games are generalisations of the prisoner's dilemma to an arbitrary number of players, see e.g. [4] and the reference therein. In the model presented in [3], N players are chosen randomly from a large population. Every round, these players may either contribute an amount c or nothing at all to a common pool. η c denotes the number of the players who cooperate and N − η c is the number of the players that defect. At every round the common pool is increased by an interest rate r and then used to pay back to the players. The payoffs for cooperators P c and defectors P d are given by for the model to be a public goods game, [4]. However, in this game it turns out that defecting is the dominating strategy. Hence, the authors of [3] proposed an extended model allowing players to decide whether to participate or not. Those who are unwilling to do so are called "loners" and they will receive a fixed payoff P l = σ c with 0 < σ < r − 1. The payoff P l ensures that an entirely cooperating group will profit more than loners while loners will profit than a group solely formed of defectors. The model of [3] thus considers three types of persons: the loners (refusing to join the group), the cooperators (who join and contribute) and the defectors (who just join). These groups correspond to payoffs P l , P c , P d and the relative frequencies of these strategies shall be denoted by x, y, z and satisfy the condition x + y + z = 1. More precisely, it was derived in [3] that The sign of P d − P c , i.e. the sign of the function G(z) plays a key factor in determining whether or not it is better to switch strategy, that is to change from deflection to cooperation or vice versa. It is straightforward to check that the function G can also written as a polynomial with real coefficients: (14) Note that for 2 < r < N those coefficients change sign exactly twice and that Descartes' rule of signs implies that G(z) has either two or zero positive roots. In fact, Lemma 7 below shows that lim z→1− G(z) = 0 from negative values. Hence, since clearly G(0) > 0, the function G(z) undergoes exactly one sign change on z ∈ (0, 1) as in Fig. 1 above.
Note that it can be easily verified that when r ≤ 2 then G(z) does not have any root in (0, 1) and G(z) = 0 ⇒ z = 1, which means that defecting is the dominate strategy.
By using the constraint x + y + z = 1, the average payoff can be written as, By introducing f = x x+y as a new variable and considering the replicator dynamicṡ z = z(σ − P), the authors of [3] obtained the following system to be considered for Note that as G(z) changes the sign once, also the factor (σ − f (r − 1)) changes its sign once according to the value of f ∈ (0, 1) and (3).
For the ODE-system (15), the authors proved in [3] that for r ≤ 2 there are no fixed points for the system in (0, 1) 2 , while when r > 2 and 0 < σ < r − 1 then there exists a unique fixed point in the interior of (0, 1) 2 , which is stable and surrounded by closed orbits. Moreover, they proved that all interior orbits are closed.
In fact, system (15) can be written as the generalised Hamiltonian system (5). We remark that in [3], the authors preformed one further transformation of the system by dividing the right hand side terms by the variable φ( f , z) as defined in (5) and then considering the resulting standard Hamiltonian system with a well-known form of prey-predator systems. However, this transformation is not necessary for our arguments but would introduce singular right hand side terms, for which already the existence of weak solutions to a corresponding PDE model are unclear. (It might be possible to define renormalised solutions).
The Hamiltonian, which transforms (15) into the Hamiltonian system (5) is given by where R(z) is defined as a primitive of ∂ R ∂z , which in return is introduced by the following definition It can be shown that ∂ R ∂z is a bounded function on z ∈ [0, 1] (see Lemma 8 below) and that the non-negativity H 2 ≥ 0 follows from choosing a sufficiently large positive integration constant in the definition of R(z).
Before we state further properties of R, we note first that and system (15) can be indeed written as a Hamiltonian system of the form (5) (5).
where the coefficients a j := r − 1 − j r N change sign exactly once between a 1 = r − 1 − r /N > 0 and a N −1 = −1 + r /N < 0, which reflects the single sign change of G(z).
Alternatively, we set N − 1 − j = k and define b k = −1 + r /N (k + 1) for for r > 2 and thus G(z) < 0 for z sufficiently close to 1.

Proof We calculate from (19)
where P(z) is a polynomial in z of order N − 3. Thus, we have which is a bounded rational function on z ∈ [0, 1].

Positive definiteness of the Hamiltonian
In the following, we need that the Hessian of the Hamiltonian is a positive definite matrix. This is obviously true if and only if ∂ 2 H 2 ∂z 2 > 0. As example, for N = 3, we calculate especially (1 + z) 2 = 0 (23) and hence positive definiteness for all z ∈ [0, 1]. For N = 4, we obtain also positive definiteness since However, the following Lemma 9 proves not only sufficient conditions for the positive definiteness of the Hessian D 2 H , but also that ∂ 2 H 2 ∂z 2 < 0 is possible. Lemma 9 (Positive definiteness of the Hessian of the Hamiltonian) Let r satisfy max{ N 3 , 2} < r < N for N ≥ 3. Then, On the other hand, for r close to 2 and N large, we find ∂ 2 H 2 ∂z 2 < 0.
Proof The proof applies different estimates on two intervals for r , first N 2 ≤ r < N and secondly max{ N 3 , 2} < r < N 2 . The presentation of the proof will be divided accordingly.
We will begin our analysis for N 2 ≤ r < N . By differentiating (19) with respect to z and using the representation (21), we derive the following formula: and as defined in (21). We notice in the range N 2 ≤ r < N that Hence, in (27) the only negative term is the middle one, i.e. b 0 (N z N −1 − 1), while the two sums in (27) contain only non-negative terms for N 2 ≤ r < N . Therefore, in order to prove ∂ 2 H 2 ∂z 2 > 0 in (26), it is sufficient to prove where the first term on the above relation is just the last term of the first sum of (27). Relation (28) with recalling b N −2 from (21) can be rewritten as Since the last term is non-negative, it is sufficient to show By observing that the above bracket is monotone increasing in r for N ≥ 3, we can estimate further below by setting r = N 2 and z = 1 to obtain The above estimates proves that the left hand side of (30) is O(1) as z → 1 for N ≥ 6 while it is O(1−z) for N = 5. Consequently the term (29) is also O(1−z). Hence, this is sufficient to prove (28) in the considered range N 2 ≤ r < N for N ≥ 5. Moreover, since positive definiteness for N = 3 and N = 4 was already shown in (23) and (24), respectively, this completes the proof in the range N 2 ≤ r < N . Now, we will treat the second range max{2, N 3 } < r < N 2 , which is more technical. Since positive definiteness for N = 3 and N = 4 was already shown in (23) and (24), we shall furthermore assume N ≥ 5. In fact, for N ≥ 5, we are able to treat the range N 3 < r < N 2 which implies max{2, N 3 } < r < N 2 . We consider .
(31) Next, we observe that in the range N 3 < r < N 2 holds Therefore The second term S 2 (z) changes sign on [0, 1] from a positive constant (at z = 0) to a negative value (at z = 1) with a single root. In order to control S 2 and S 3 and show (26), we only need to consider the last two terms of both S 1 and S 4 . Hence, we shall prove the non-negativity of the following partial sum, denoted by S p (z), as a sufficient condition From these six expressions only S 2 and S 3 can be negative. More precisely, Accordingly, in order to prove (32), we split the interval [0, 1] into the intervals I = [0, z * ], where instead for proving (32), it will be sufficient to show S 1a + S 1b + S 3 > 0 on (0, z * ] and the interval II = [z * , 1], where we require all six terms of S p to prove the sufficient condition (32).
We begin with the first interval I = [0, z * ]. Since S 2 ≥ 0 on I, it is sufficient to show S 1a + S 1b + S 3 > 0 on z ∈ (0, z * ], i.e. (33) First, we observe that these three terms are all monotone increasing in r for N ≥ 4. With r > max{ N 3 , 2} ≥ 2, it is thus sufficient to set r = 2 in (33) in order to prove (32) on the interval I. Hence, after multiplying with N , we obtain the sufficient condition which proves (33). Finally, since S p (0, N , r ) = −b 0 > 0, we have (32) in the interval I.
We continue with the second interval II = [z * , 1]. Here, since N z N −1 − 1 ≥ 0 and S 2 ≤ 0, the fact that the coefficients b k in (31) are monotone increasing in r implies that all six terms in S p (z, N , r ) are monotone increasing in r . Thus, it is sufficient to prove S p (z, N , 2) > 0: By collecting the terms proportional to 1 − 2 N and 1 − 4 N , respectively, we obtain Hence, for N ≥ 5, we estimate first which implies the strict inequality R 1 > 0 for all z ∈ [z * , 1) and lim z→1− R 1 (1−z) 2 > 0. Secondly, we calculate which implies the strict inequality R 2 > 0 for all z ∈ [z * , 1) and lim z→1− This proves (32) on the Interval II. Finally, Fig. 2 illustrates in the limiting case r = 2 the sign of ∂ 2 H 2 ∂z 2 (by plotting 2 in order to avoid plotting singularities at z = 0, 1). We observe that for 16 N , there are values near z ∼ 0.7 where ∂ 2 H 2 ∂z 2 < 0.

Proof of Theorem 1
We begin by a brief outline of the proof of Theorem 1: first, we reformulate (1) in terms of the Hamiltonian structure of ODE model (5). Then, we notice that the spatially integrated ODE Hamiltonian constitutes also a Lyapunov functional for the PDE system (1), i.e.
Next, we prove that the ω-limit set is non-empty, compact and connected. Afterwards, for trajectoriesw of the ω-limit set, we show that H (w) dx is well defined. and we get that the ω-limit set is invariant. Finally, we notice that asymptotically, we have a spatially homogeneous and periodic in time orbit.
Hamiltonian structure and global classical solutions Recalling the ODE model (5), we rewrite the PDE system (1) as where we denote Straightforward calculation yields where (36) holds provided that H zz > 0 (see e.g. Lemma 9 ) because Hence, H( f , z) constitutes a Lyapunov functional for the PDE system (1). Moreover, the Lyapunov functional H( f , z) ≥ 0 is non-negative since H 1 ≥ 0 and H 2 ≥ 0. Next, we remark that global-in-time solutions to the PDE system (1) follow from standard parabolic theory (e.g. invariant regions, see [11]) or weak comparison principle arguments ensuring 0 ≤ f (x, t), z(x, t) ≤ 1 a.e. x ∈ for all t ≥ 0 provided that 0 ≤ f (x, 0), z(x, 0) ≤ 1. Moreover, standard parabolic regularity implies classical C 2 solutions due to the assumed regularity of ∂ .
Orbits and ω-limit set We define the solutions orbits of the PDE system (1): which is compact and connected in (C 2 ( )) 2 .
At this point we can see that H(w) = H (w)dx is not necessarily well-defined for w = w ∞ ∈ ω(w 0 ) since we only know 0 ≤ w ∞ ≤ 1. In order to exclude that possibility, we consider for any w ∞ ∈ ω(w 0 ) a solution w = w(·, t) to the PDE system (1) for t ≥ 0 with lim t→∞ w = w ∞ . Since we assumed initial data with H(w(·, 0)) ≤ C, by using Fatou's Lemma, the following arguments shows w ∞ ≡ 0, 1: With H 1 and H 2 being non-negative, see (17), (18), this implies The first bound writes as Similarly we derive that z ∞ ≡ 0, 1 and thus conclude that w ∞ ≡ 0, 1. Returning to the above solutionw =w(·, t) subject tow| t=0 =w 0 ∈ ω(w 0 ), it follows that at some points x ∈ , the initial dataw 0 Thus, the strong maximum principle implies that for all t > 0 holds 0 <w(·, t) < 1 and thus H(w) is a well defined for t > 0.

Proof of Corollary 3
The proof of Corollary 3 is based on the ideas of the second part of Theorem 1.1 in [5], which we outline for the convenience of the reader. First, we prove the following for any T > 0. We begin by observing that the previous Theorem 1 and parabolic regularity implies for any solution ( f (·, t), z(·, t)) to (1) the existance of a constant C > such that for positive times, for instance, for t ≥ 1 holds Then, by the theorem of Ascoli-Arzelá, the sequence {t k } ↑ ∞ admits a subsequence {t k } ⊂ {t k } and a solution (f (·, t),ẑ(·, t)) to (1) such that for any T > 0 (cf. [5, Lemma 3.6]) satisfies From relation (37), we derive that Hence, the solution (f ,ẑ) must be spatially homogeneous and, thus, it is a solution to the ODE model (15) and we denote it by (f (t),z(t)) in the following. Moreover, (38) follows from (39).
Next, we apply semi-group estimates for the Laplace operator subject to homogeneous Neumann boundary conditions, see e.g. [10,13] ∇e t φ s ≤ C(q, s)e −μt max{1, t where μ 2 denotes the second eigenvalue of − (with the Neumann boundary conditions). Actually, the first equation in (1) implies where Hence, by taking the gradient of (43), it follows from (42) that for q = ∞ = s, Note that again the constant C > 0 is independent of d z > 1. Similar, we consider where and hence ∇z(·, t) ∞ ≤ C, t ≥ 0.

Proof of Theorem 6
Proof of Theorem 6 The shadow system (9) takes the form First, we verify that the Hamiltonian (16)-(18) also applies to the shadow system (47): for all (∇ f , ∇z) = 0, since Z = Z (t) and H z (Z ) are spatially homogeneous. Given the same Hamiltonian as for the PDE model (1), we can follow all the steps of the proof of Theorem 1 to prove Theorem 6. In particular, global existence of solutions to the shadow system (47) and the characterisation of the ω-limit set can be performed in the same way.