Nonlinear Dynamics in 2-Population Mean Field Games Competing for a Resource: Influence on the Tragedy of the Commons

Abstract Mean field Game (MFG) Theory derived by combining stochastic optimal control with a statistical description of multiple populations of competing agents is increasingly used in many applications, especially those involving acquisition of a resource. Using an ergodic, 2-population MFG model it is demonstrated that nonlinear dynamics caused by the model’s underlying bifurcation structure can dominate predictions as the influence of the optimal control is enhanced by lowering the stochastic diffusivity \(\sigma\). Multiple ergodic states are predicted, and rapid segregation of the distributions is observed over small ranges of \(\sigma\) close to bifurcation points. Introducing bias toward agents with more expertise breaks the bifurcation structure, however the effects of the bifurcations persist in the resulting multiple disconnected states. Including a common resource that is self-regenerating models the Tragedy of the Commons (TOTC) where the resource is depleted by over-acquisition by the competing populations. The dependence of the TOTC limit on \(\sigma\) and the regeneration rate \((R_{N} )\) exhibits the residual influence of the bifurcation structure and controls the sensitivity of the limiting value of \(R_{N} .\) Large differences in the maximum resource acquisition between the two populations are predicted. Similar content being viewed by others 1 Introduction Understanding the competition between groups to acquire a resource is of fundamental importance in a variety of situations in society. Whether it is the competition for a natural resource, customers competing in a market-place or for resources distributed by a government, the results of these competitions effect both the competitors and, in many instances, the availability of the resource. The concept of the Tragedy of the Commons (TOTC) popularized by Hardin [1] (also see [2, 3]) is of special importance. Here the self-interested actions of populations working to acquire a resource lead to its overuse, potentially to the point of its total depletion. The purpose of this work is to investigate the influence of nonlinear interactions on the results of Mean Field Game (MFG) Theory applied to the competitions between two groups of agents or populations competing for a given resource. The case of the TOTC is described by a self-regenerating resource where excessive acquisition by the populations leads to its depletion. There is a rich literature on modeling the TOTC using agent-based models, systems dynamics descriptions, and reaction–diffusion models. For example, reaction–diffusion models have been extensively studied for modeling competitions [4] using simple linear kinetics and more complex descriptions. Mean Field Game theory was originated by [5,6,7] and Huang et al. [8] to extend the concept of a game with a finite number of competing players or agents to a large number or continuum of agents each following the same game strategy. MFG Theory models are derived by coupling a stochastic differential equation describing the actions by a population of agents with their trajectory determined by stochastic optimal control. The result is two coupled nonlinear partial differential equations describing the evolution of the population: the Hamilton–Jacobi-Bellman (HJB) equation for the value function for control and the Fokker–Planck-Komogorov (FPK) equation for the distribution function of the population. Both equations are written in terms of time and continuous independent variables (denoted here by the single variable x) with meaning that depends on the model for the contest or game. For games with agents competing for a resource, x is interpreted as the variable effort or expertise of an agent at acquiring the resource, with \({ }0 \le x \le 1\). Mean field game theory formulations are typically time dependent having the complex structure that the FPK equation is expressed forward in time while the HJB equation is expressed backward in time over the interval \(0 \le t \le T\), where \(T\) is the time horizon for control. The analysis here focuses on stationary ergodic states corresponding to \(T \to \infty\) [9, 10]. The framework of MFG is well suited to addressing questions of population competing for a resource as it has been shown that MFG’s yield a continuum approximation to a Nash equilibrium for a finite population where no agent would gain by changing her strategy if the strategies of all the other agents stay constant. This connection to a Nash equilibrium gives special meaning to the results described here. Since its introduction there has been an explosion in their development and application of MFG theory; see the volumes and references in [11,12,13]. Mean Field Game Theory also has been applied recently to modelling TOTC [14,15,16,17]. The work of Dayaniki and Laurière [14] is especially relevant here as these authors presented a framework for single and two-population games for a common resource that includes both competition and cooperation between the populations. Their analysis is based on solution in time of forward and backward stochastic differential equations (FBSDEs); however, the authors confined there analysis to cases where uniqueness of solutions is guaranteed. The calculations here demonstrate dramatic effects of nonlinear interactions on the predictions for the TOTC. Two developments in MFG theory are of particular importance to this work. The first was the extension of MFG Theory to multiple populations with distinct strategies and, in our case abilities, for acquiring the resource [18,19,20]. An important consequence of considering two-populations in MFG theory was discovered by Achdou et al. [21] who demonstrated conditions where the two populations segregate in the independent variable, i.e., the nonzero values of the two distributions separate in x (this concept will be made mathematically rigorous below). The second important development was made by Cirant and Verzini [22] who recognized that even simple, two-population mean field games can have multiple ergodic (quasi-time-independent) solutions and that the existence of these solutions is connected to segregation when the control strategy dominates the stochastic dynamics of the populations. The question of temporal stability is key to classifying the multiple ergodic stationary states but is complicated by the model’s origins in stochastic optimal control. Guéant et al. [23] introduced two concepts of stability. They defined “physical stability” as the property when perturbations to a stationary state introduced at the temporal boundaries decay over time. The forward/backward temporal structure of the HJB/FPK equations sets these conditions at \(t = T \) for the value function and \(t = 0 \) or the distribution function. Guéant [24] proved physical stability for a one population game with a unique solution. Eductive stability, the second concept introduced in Guéant et al. [23], refers to the dynamics of a modified dynamical system in which the backward time dependence of the HJB equation is reversed. Then a stationary solution of the original MFG equations is said to be eductively stable if perturbations that evolve forward in time that are “near” the stationary state converge to it. This concept is especially useful for numerical calculations of stationary states by time-integration in this fictious time variable [21]. Both stability concepts are generalized here to ergodic states computed in the limit \(T \to \infty\), where the indefinitely large time horizon makes it reasonable to compute the long-time stability of a small perturbation to the stationary ergodic state using the classical methodology applied to stability analysis of dynamical systems [25]. This approach is adapted here to account for the forward/backward time evolution of the MFG equations. The purpose of this paper is to use a combination of bifurcation analysis and numerical calculations to demonstrate that the interactions of the multiple, ergodic solutions of the two-population MFG qualitatively impact predictions of population segregation and resource consumption. The organization of the paper is as follows. The MFG theory model is presented next for a two-population game with competition framed for three different models of resource availability: (i) the amount of resource available is constant, independent of the game, (ii) the resource is externally replenished, and (iii) the resource is replenished naturally at a rate proportional to its availability. Case (i) results in the analysis being independent of the resource and allows connection to previous results [22],case (ii) models two populations competing for a finite resource Brown [26] discussed a one population version of this model),and case (iii) corresponds to a two-population model appropriate for studying the TOTC. Asymptotic results are summarized as a basis for understanding the numerical results computed by finite element approximations and Newton’s method for solution of nonlinear boundary value problems, coupled with numerical bifurcation analysis for computing ergodic states, linear stability analysis of these states, and nonlinear time integration. 2 Mean Field Models A MFG mathematical model is developed for two populations of agents who are competing to maximize their acquisition of a common resource. Each population is represented by a continuous distribution where the actions of individual agents are part of the continuum. The MFG model is composed of two coupled partial differential equations for each population, consisting of a Hamilton–Jacobi-Bellman equation (HJBE) derived to maximize a functional of a stochastic optimal control problem and a Fokker–Planck-Kolmogorov equation (FPKE) that defines the distribution derived using Ito calculus from the stochastic differential equation (SDE) for the optimal control. We write the SDE in terms of a stochastic variable \(x_{t}\), \(0 \le x_{t} \le 1,\) that measures the effectiveness of an agent at procuring the resource. For the i-th population, the SDE is where the control trajectory \(\alpha_{i} \left( {x,t} \right)\) is determined by the optimal strategy to maximize the utility functional, and \(dW\) is additive Brownian noise with intensity or diffusivity \(\sigma\) that is assumed the same for both populations. Increasing \(\sigma\) increases the randomness of the agent’s actions and lessens the effectiveness of the control. The control trajectories and the drift velocity in the FPKE are determined from the Hamiltonian defined for the i-th population as: where \(p_{i}\) is a generalized momentum variable and \(S\left( {\alpha_{i} } \right)\) is the cost of control. Defining \(U_{i} \left( {x,t} \right)\) as the value function that maximizes the utility function for the i-th population gives \(p_{i} \equiv \frac{{\partial U_{i} }}{\partial x}\). The analysis is based on the much-used quadratic expression for the cost of control [11, 27] \(S\left( {\alpha_{i} } \right) \equiv \frac{{\alpha_{i}^{2} }}{2\gamma }\), where \(\gamma\) is a constant taken here as \(\gamma = 1\). Then from Eq. (2) the optimal control strategy is \(\alpha_{i}^{*} \left( {x,t} \right) = b\left( {x,t,\alpha_{i}^{*} \left( {x,t} \right)} \right) = \gamma^{ - 1} \left( {\partial U_{i} /\partial x} \right)\) and also is the drift velocity in the FPKE. The time-dependent HJBE’s for the value functions are where \(L_{i} \left( {x,t,N,W_{i} } \right)\) models the net effort of an agent, \(W_{i}\)(x,t) is the distribution function, and the variable \(N\left( t \right)\) measures the amount the resource available for acquisition by both populations. Equations (3a) are solved backward in time from an initial condition defined at the time horizon \(t = T\), \(0 \le t \le T\). The value function \(U_{i} \left( {x,t} \right)\) is determined only up to an additive constant which is found by setting The corresponding FPKE’s are where \(\lambda_{Wi} \) is a Lagrange multiplier included to satisfy the normalizing constraint on the distribution The formulation is completed with Neumann boundary conditions by setting \(\partial U_{i} /\partial x \) and \(\left( {\left( {W_{i} /\gamma } \right)\partial U_{i} /\partial x - \left( {\sigma^{2} /2} \right)\partial W_{i} /\partial x} \right)\) to zero at \(x = 0\) and \(x = 1\). Integrating (3) over \(0 \le x \le 1\) and applying (3d) and the Neumann conditions leads to the result \(\lambda_{Wi} = 0\). These variable are retained in the numerical analysis to enforce the constraint (3d) but dropped from the analysis that follows. The analysis is based on computing stationary solutions of Eq. (3) in the ergodic limit \(\left( {T \to \infty } \right)\) where \(W_{i} \left( {x,t} \right) = W_{si} \left( x \right)\) and \(U_{i} \left( {x,t} \right) = U_{si} \left( x \right) - \lambda_{Ui} t\) so that the stationary solution to Eq. (3) is defined by the six variables \(\left( {W_{si} \left( x \right),U_{si} \left( x \right),\lambda_{Ui} } \right), i = 1,2\). Going forward, the subscript “s” denoting the stationary solution is deleted. The stationary problem is: with the boundary conditions described above. 3 Modeling Competitive and Cooperative 2-Populations The utility functions \(\left\{ {L_{i} } \right\}\) are written as where the functions for the rate of acquisition, \(A_{i} \left( {x,W_{i} \left( x \right),N} \right),\) and its cost, \(C_{i} \left( {x,W_{i} \left( x \right),N} \right)\) are selected to model several distinguished cases. The general forms of these functions are: where \(\left( {q_{oi} ,c_{oi} } \right)\) are constants that scale the functions. The dependence on the available resource N is only included in the acquisition to model an increased acquisition rate as the resource becomes more plentiful. No dependence on N is included in the cost function which only includes the cost attributed to competition from the competing distribution; however, the competing distribution implicitly depends on N through its acquisition function. Also, the form of Eq. (3b) allows direct comparison with previous work. The constants \(\left\{ {\gamma_{ij} } \right\}\) set the relative influence of the distribution functions in \(A_{i} \) and \(C_{i}\) and \(\left\{ {\beta_{ai} ,\beta_{ci} } \right\}\) set their exponential dependences. Cooperation between members of the i-th population is model by \(\gamma_{ii} > 0\) which introduces \(W_{i} \left( x \right)\) into \(A_{i}\). Conversely, competition between the two populations corresponds to \(\gamma_{ij} > 0, i \ne j\). Altruistic cooperation between the two populations would be modelled by including \(W_{j} \left( x \right)\) in \(A_{i}\) \(\left( {i \ne j} \right) \) and is not studied here; Dayaniki and Laurière [14] include this effect. Bias toward agents with more expertise is included in the acquisition and cost functions through the functions \(h_{ai} \left( x \right)\) and \(h_{ci} \left( x \right)\); for example, \(dh_{ai} /dx > 0\) corresponds to biasing acquisition toward higher effectiveness x and is equivalent to introducing nonhomogeneous dependence on the independent variable x. These functions serve as imperfections to the bifurcation structure for a perfect (\(\frac{{dh_{ai} }}{dx} = \frac{{dh_{ci} }}{dx} = 0)\) 2-population game. The model is completed by adding a balance equation for the resource. Motivated by McPike [28] N(t) is modeled as where \(R_{N}\) is the rate of production of the resource. The expected values \(E_{i} \left[ \cdot \right] \equiv \mathop \smallint \limits_{0}^{1} \cdot W_{i} \left( x \right)dx\) measure the acquisition of the resource by the two populations and introduce a non-local coupling in \(W_{i} \left( x \right)\). The first two terms in Eq. (7) represent the logistic or Verhulst rate equation, which, without acquisition by the populations, predicts a sigmoidal shaped evolution of N between zero and its equilibrium value \(N_{eq}\). In Eq. (7), \(f_{o} \) models a constant rate of addition of the resource and the term \(\left\{ \cdot \right\}\) accounts for the depletion of the resource due to acquisition by the two populations. Equation (7) assumes that the sizes of the populations are equal; different population sizes can be introduced by including another factor or by adjusting the relative magnitudes of \(\left( {q_{o1} ,q_{02} ,c_{1} ,c_{2} } \right)\). The ergodic version of Eq. (7) assumes \(\partial N/\partial t = 0\). Without external input \(f_{o}\), the TOTC occurs when \(N \to 0\) as the total resource acquisition \(A_{tot} \equiv \left\{ {E_{1} \left[ {q_{1} } \right] + E_{2} \left[ {q_{2} } \right]} \right\}\) reaches this limit as \(\mathop {\lim }\limits_{N \to 0} \left( {A_{tot} /N} \right) \to R_{N}\). Including \(f_{o}\) and setting the rate of production to zero, \(R_{N} = 0,\) gives \(N = \frac{{f_{o} }}{{\left( {A/N} \right)}}\), the form used in [26] to model a resource that is replenished at a constant rate. Simply setting \(N = N_{o}\) models a constant resource level that is invariant to acquisition. The analysis here focuses on these three cases. Calculations with \(N = N_{o}\) connect to the limit studied in [22] where the existence of multiple solutions and segregation for decreasing diffusivity \(\sigma\) are demonstrated. Computations demonstrate the sensitivity of these solution to relatively small changes in \(\sigma\). Introducing the nonhomogeneous factor \(h_{ai} \left( x \right)\) in the acquisition function breaks the bifurcation branches and demonstrates their lingering effects on the solution structure with varying \(\sigma\). Calculations with the constant replenishment rate \(\left( {f_{o} > 0} \right)\) demonstrate differential resource acquisition driven by differences in acquisition and cost functions and by the effects of the underlying nonlinear structure. Finally, the effect of the bifurcation on the TOTC is revealed by calculations with \(f_{o} = 0 \) and \(R_{N} > 0\). The focus is on parameter values where the solutions to (4) are not unique. Following Cirant [19], uniqueness is guaranteed when \(L_{i} \left( {x,W_{i} \left( x \right),W_{j} \left( x \right),N} \right), j \ne i,\) satisfies where \(\left\{ {W_{i} \left( x \right)} \right\}\) and \(\left\{ {\overline{W}_{i} \left( x \right)} \right\}\) are any continuous functions on \(x \in \left[ {0,1} \right]\) that integrate to unity. The derivation of this condition follows from postulating two unique solutions \(\left( {\overline{W}_{i} ,\overline{U}_{i} ,\overline{\lambda}_{i} } \right) and \) \(\left( {\hat{W}_{i} ,\hat{U}_{i} ,\hat{\lambda}_{i} } \right) \) of Eq. (4), multiplying Eq. (4a) by \(\left( {\overline{W}_{i} \left( x \right) - \hat{W}_{i} \left( x \right)} \right)\) and summing over the two populations. Using integration by parts and Eq. (4c) reduces the sum to two terms; as an always positive term proportional to \(\left( {d\overline{U}_{i} /dx - d\hat{U}_{i} /dx } \right)^{2}\) and the expression (8). It follows that if (8) is negative-indefinite then \(d\overline{U}_{i} /dx = d\hat{U}_{i} /dx \) and applying (4b) that \(\overline{U}_{i} \left( x \right) = \hat{U}_{i} \left( x \right) \). Substituting this result into the FPPEs it is easily shown that \(\overline{W}_{i} \left( x \right) = \hat{W}_{i} \left( x \right)\) . So that the solution must be unique. Equation (8) is satisfied when the symmetric matrix with elements \(\left\{ {\frac{{\partial L_{i} }}{{\partial W_{j} }} + \frac{{\partial L_{j} }}{{\partial W_{i} }}} \right\}\) is negative definite. Using the definitions (6a-6b) this condition yields the criteria: The first condition decreases the resource acquisition (6a) with increasing \(W_{i} \left( x \right){\text{ and is uninteresting in the context of this work}}.\) The second condition also is satisfied by setting \(\gamma_{21} = \gamma_{12} = 0\) which reduces the 2-population system to two independent (neglecting potential coupling through N) populations. This case is not considered further. 3.1 Critical Points for Bifurcation When nonhomogeneous effects are excluded (\(h_{ai} \left( x \right) = h_{ci} \left( x \right) = 1)\), the MFG model, Eqs. (4–7) admits the trivial solution: where \(N_{0}\) is computed from Eq. (7). The three scenarios for resource replenishment correspond to. - (i) N = \(N_{0} ;\) - (ii) \(f_{o} > 0\) and \(R_{N} = 0; \quad{N_{0} = \frac{{f_{o} }}{{\left[ {q_{o1} \left( {1 + \gamma_{11} } \right) + q_{o2} \left( {1 + \gamma_{22} } \right)} \right]}}}\) - (iii) \(f_{o} = 0\) and \(R_{N} > 0\); \(\quad{N_{0} = N_{eq} \left\{ {1 - \frac{{\left[ {q_{o1} \left( {1 + \gamma_{11} } \right) + q_{o2} \left( {1 + \gamma_{22} } \right)} \right]}}{{R_{N} }}} \right\}}\) Weakly nonlinear solutions of the MFG model bifurcating from the trivial state are computed as a function of \(\sigma\) as a Taylor series in amplitude \(\varepsilon\) emanating from the critical values \(\left\{ {\sigma_{on} } \right\}\) [25] written as Taylor series: Substituting (9) into the MFG model and collecting terms of \(O\left( \varepsilon \right)\) gives with Neumann boundary conditions on \(U_{i1} \left( x \right) and W_{i1} \left( x \right). \) The constants in Eq. (12a) are defined as \(a_{i} \equiv q_{oi} \gamma_{ii} N_{0}\) and \(c_{i} \equiv c_{oi} \gamma_{ij} \beta_{ci}\). The equation for the \(O\left( \varepsilon \right)\) correction \(N_{1}\) is derived by expanding the stationary version of Eq. (7) in ϵ. Simplifying these expressions using Eqs. (12b) and (12d) yields \(N_{1} = 0\) in all three cases. Following an early application of bifurcation analysis to a problem including Lagrange multipliers [29], it is useful to express Eqs. (12) in operator notation using \(D^{2} \equiv d^{2} /dx^{2} , \Sigma \equiv \frac{{\sigma_{0}^{2} }}{2}, \Gamma \equiv \gamma^{ - 1} , and {{\varvec{\Phi}}}^{T} \left( {U_{1} \left( x \right),U_{2} \left( x \right),W_{1} \left( x \right),W_{21} \left( x \right),\lambda_{U1} ,\lambda_{U2} } \right)\): with \(\left( {{\Phi }\left( x \right),{\Psi }\left( x \right)} \right) \equiv \mathop \smallint \limits_{0}^{1} {\Phi }\left( x \right){\Psi }\left( x \right)dx\). Defining the inner product as \({{\varvec{\Phi}}},{{\varvec{\Psi}}} = \mathop \sum \limits_{i = 1}^{6} \left( {{\Phi}_{i} ,{\Psi}_{i} } \right)\), where \(\left\{ {{\Phi}_{i} } \right\}\) are components of \({{\varvec{\Phi}}}\), gives the adjoint problem \({\varvec{L}}^{\varvec{*}} \left( {{\varvec{\Phi}}} \right) = 0\) and the adjoint operator \({\varvec{L}}^{\varvec{*}} \left( {{\varvec{\Phi}}} \right):\) Equations (12a) and (12c) are satisfied by solutions of the form which also satisfy the integral constraints (12b) and (12d), so that \(\lambda_{Ui}^{1} = 0, i = 1,2\). Substituting eqs. (15) into eqs. (12a) and (12c) gives homogeneous equations for (\(\alpha_{n} , \beta_{n} )\) that have a nontrivial solution only if yielding the bifurcation points or, if the level of self-cooperation is the same in each population \(\left( {a_{1} = a_{2} = a} \right)\), Taking the positive root with a = 0 and \(\beta_{ci} = 2\) reduces Eq. (17b) to the expression derived by Cirant and Verzini [22] once differences in the definition of \(\sigma\) are considered. The multiplicity of eigenvalues (17) suggest a complex eigen-structure of the operator L which is seen by computing the eigenvalues and eigenvector of \({\varvec{L}}\left( {{\varvec{\Phi}}} \right) - \mu {{\varvec{\Phi}}}\) with \({{\varvec{\Phi}}}\) given by (15). The eigenvalues \(\left\{ {\mu_{n} } \right\} \) are computed as For \(\sigma > \sigma_{01}\), the first critical point, two real eigenvalues of opposite sign pass through zero at \(\sigma = \sigma_{01}\) and become purely imaginary for \(\sigma < \sigma_{01}\). The bifurcation point \(\sigma_{01}\) is a Riesz index two double point with two zero eigenvalues, one eigenfunction and one generalized eigenfunction (p. 104 [25]) The single eigenfunction is described by the coefficients in Eq. (15) as \(\alpha_{2n} = {\Lambda }\alpha_{1n} , \beta_{in} = \left( {2/\sigma_{on}^{2} \gamma } \right)\alpha_{in} , \) and \(\lambda_{Uni} = 0, i = 1,2., {\text{where}} {\Lambda } \equiv \left( {c_{o1} \gamma_{21} /c_{o2} \gamma_{12} } \right)^{1/2}\). The generalized eigenfunction and the eigenfunction and generalized eigenfunction for the adjoint problem have forms similar to Eq. (15), reducing each problem to a 4 × 4 algebraic set for the coefficients \(\left\{ {\alpha_{n} ,\beta_{n} } \right\}\). Cirant and Verzini [22] applied bifurcation analysis appropriate for a simple eigenvalue to the case \(a_{1} = a_{2} = 0 \) and \(\beta_{ci} = 2 {\text{and }}\) predicted supercritically bifurcating branches—evolving to lower values of \(\sigma\). Their analysis is recovered by extending the series solution Eq. (9) to higher orders in \(\varepsilon\). It is not reproduced here; instead, the bifurcation structure is computed numerically. The bifurcation structure implied by their analysis is shown schematically in Fig. 1. Understanding the implications of the bifurcations on the stability of the ergodic states is complicated by the forward/backward structure of the underlying time-dependent equations. The linearized versions are written as: where \({\varvec{D}}\) is a diagonal matrix with the elements \(\left( { - 1, - 1,1,1,0,0} \right)\). Equation (19) is converted to an eigenvalue problem by introducing small amplitude perturbations of the form \(i = 1,2, \) where the forward/backward time integration is accounted for by reversing the exponential time dependences of \(U_{i1} \left( {x,t} \right)\) and \(W_{i1} \left( {x,t} \right)\). Substituting (19) into (18) and computing \(\left\{ {\mu_{n} } \right\}\) so that the homogeneous system has a nontrivial solution yields Expression (21) describes the “physical stability” of the trivial state to infinitesimal disturbances distributed across \(0 \le x \le 1\). For \(\sigma > \sigma_{01}\), \(\widehat{\mu}_{n}^{ \pm } > 0,\) so that the disturbances to both \(U_{i1} \left( {x,t} \right)\) and \( W_{i1} \left( {x,t} \right)\), traveling in opposite directions in time, decay. At \(\sigma = \sigma_{01} , \widehat{\mu}_{n}^{ - } = 0,\) and for \(\sigma < \sigma_{01} t\) he trivial state is unstable. These are the characteristics of a simple bifurcation [25, 30]. The eductive stability of the trivial branch is computed by changing the signs of time derivatives in the HJB Eq. (3a) to give the diagonal in D \(\left( {1,1,1,1,0,00,0} \right), \) yielding differential structures similar to the heat equation. Performing linear stability analysis for these equations with disturbances with exponential time dependence of the form \(exp\left( { - \mu_{n} t} \right), t > 0, \) and for all variables gives the same result Eq. (21), predicting no difference between the two concepts of stability in the ergodic limit. Although the complete asymptotic bifurcation analysis is not reproduced here there are several analytical results that inform the numerical analysis. The eigenfunction and generalized eigenfunctions for the original and adjoint problems computed at \(\sigma = \sigma_{o1}\) have forms similar to Eq. (15) with coefficients proportional to \({\text{cos}}\left( {\pi x} \right)\). The effect on the bifurcation structure of adding x-dependency in the acquisition function is computed with \(h_{ai} \left( x \right) = \left( {1 + \delta \left( {x - 1} \right)} \right), \delta \ll 1\). Expanding \(\delta = \delta \left( \varepsilon \right)\), in addition to the variables in Eq. (11), leads to the inhomogeneous equation of the form \({\varvec{L}}\left( {{\varvec{\Psi}}} \right) = {\varvec{d}}_{{\varvec{\delta}}}\), where \({\varvec{d}}_{{\varvec{\delta}}}^{{\varvec{T}}} = \delta_{1} \left( {q_{01} f_{o} \left( {x - 1} \right),0,q_{02} f_{o} \left( {x - 1} \right),0,0,0} \right)\) and \(\delta_{1}\) is the \(O\left( \varepsilon \right)\) coefficient in the expansion for \(\delta \left( \varepsilon \right)\). The solvability condition at \(\sigma = \sigma_{01}\) is \({\varvec{d}}_{{\varvec{\delta}}} ,{\Phi}_{11}^{*} \left( x \right) = 0\), where \({\Phi}_{11}^{*} \left( x \right)\) is the generalized adjoint eigenfunction corresponding to \(n = 1\). Because all nonzero terms in \({\Phi}_{11}^{*} \left( x \right)\) are proportional to \(\cos \left( {n\pi x} \right)\), the nonzero integrals in the inner product are proportional to \(\mathop \smallint \limits_{0}^{1} \left( {x - 1} \right)\cos \left( {\pi x} \right)dx = - 2/\pi^{2}\) so that the solvability condition gives \(\delta_{1} = 0\) and the bifurcation structure is singularly perturbed (often referred to as broken) by introducing \(\delta\). Introducing \(h_{ai} \left( x \right) = \left( {1 + \delta \left( {x - 1} \right)} \right), \delta \ll 1,\) does not affect the structure of the branches bifurcating from \(\sigma = \sigma_{02}\). A similar calculation of the solvability condition \({\varvec{d}}_{{\varvec{\delta}}} ,{\Phi}_{12}^{*} \left( x \right) = 0,\) where \({\Phi}_{12}^{*} \left( x \right)\) is the generalized adjoint eigenfunction corresponding to \(n = 2\), yields integrals of the form \(\mathop \smallint \limits_{0}^{1} \left( {x - 1} \right)\cos \left( {2\pi x} \right)dx = 0\), leaving \(\delta_{1}\) undefined and suggesting that the structure is unperturbed by \(h_{ai} \left( x \right)\). The role of introducing bias on the bifurcation structure also is shown qualitatively in Fig. 1 for \(\delta \ll 1.\) 3.2 Numerical Analysis The coupled boundary-value problems, integral constraints, and resource equations, Eq. (4) and (7), are solved numerically by combining a finite element-Newton method with numerical bifurcation analysis. Equations (4) are written in Galerkin weak form using a linear finite element discretization in \(n_{d} - 1\) elements (\(n_{d} = 100)\) with the Neumann boundary conditions incorporated naturally. The resulting nonlinear algebraic equations are represented as where \({\varvec{u}}\left( {{\varvec{p}},\sigma } \right)\) is the solution vector in terms of parameters \({\varvec{p}}\) and diffusivity \(\sigma .\) From a first guess \(\left( {{\varvec{u}}^{\left( 0 \right)} } \right)\) Newton’s method was used to solve Eq. (22) according to the iterations where \(R_{u} \left( {{\varvec{u}}^{\left( k \right)} ;{\varvec{p}},\sigma } \right) \) is the Jacobian matrix; the linear equations were solved by LU factorization using MatLab [31]. Critical points were detected by computing the eigenvalues of the matrix \({\varvec{R}}_{{\varvec{u}}} \left( {{\varvec{u}};{\varvec{p}},\sigma } \right)\) evaluated along a solution branch; bifurcation points along the trivial state were located with an accuracy better that 0.1 percent. Pseudo-arc-length continuation [32] was used to compute solution branches and navigate turning points with respect to \(\sigma\). Here the nonlinear system is inflated by parametrizing \(\sigma = \sigma \left( s \right)\), where s is the pseudo-arc-length. An approximate unit tangent vector \(\left( {u\left( {\dot{s}} \right),\sigma \left( {\dot{s}} \right)} \right)\) to the solution curve is defined as which is solved as described by Chan [32] using the LU factorization \({\varvec{R}}_{{\varvec{u}}} \left( {{\varvec{u}}_{o} ;{\varvec{p}},\sigma_{o} } \right)\) computed at \(s = s_{o} .\) The arc-length s is related to \(\sigma \left( s \right) \) using the expression suggested by Keller [33]: where \({\Delta }s\) is varied to control the step size in \(\sigma ; \) Eq. (25) is solved simultaneously with Eqs (20). A prediction of an adjacent solution along the branch is \(u\left( {s_{o} + \Delta s} \right) = u_{o} + \Delta s\dot{u}_{o} and \sigma \left( {s_{o} + \Delta s} \right) = \sigma_{o} + \Delta s\dot{\sigma}_{o} .\) Newton’s method implemented for the coupled equations typically converged in 3–5 iterations. The eigenvalues and eigenvectors of both the linearized equations and the eigenvalue problem defining physical stability are computed for computed stationary states. The eigen-structure at the critical points matches Eq. (21); for example, for \(\sigma > \sigma_{01}\) and decreasing there are two real eigenvalues \(\left( {\mu_{1} ,\mu_{2} } \right)\) satisfying \(\mu_{1} + \mu_{2} = 0\) that pass through zero at \(\sigma = \sigma_{01}\) to become purely imaginary for \(\sigma < \sigma_{01}\). A single null vector \({\varvec{z}}_{1}\) of \({\varvec{R}}_{u} \left( {{\varvec{u}};{\varvec{p}},\sigma } \right) \) corresponds to both zero eigenvalues at \(\sigma = \sigma_{01}\) giving the expected index 2 degeneracy. The corresponding adjoint eigenvector \(and \) the generalized eigenvectors, are computed as (p. 52 [25]): Solutions on bifurcating branches near \(\sigma = \sigma_{0n}\) can be computed by perturbing the first approximation for the Newton iteration by \({\varvec{z}}_{1}\) and varying its amplitude to drive convergence to a solution on the bifurcating branch. Alternatively, the imperfection \(h_{ai} \left( x \right) = x\) can be introduced in Eq. (6a) to move the calculation directly to the perturbed bifurcating branch. Removing the imperfection (reducing \(h_{a} \left( x \right) \) to \(h_{a} \left( x \right) = 1)\) results in the desired unperturbed state. Numerical calculations are performed for both physical and eductive stability of the ergodic states. The finite element discretization of the time dependent version of the MFG model is expressed as a system of differential algebraic equations (DAEs) as where M is the singular mass matrix with components corresponding only to the value functions, distributions, and the resource quantity N(t). \({\varvec{M}}\) was varied to correspond to both forms of the stability analysis. Linear stability analysis was performed by converting Eq. (25) to a generalized eigenvalue problem for the eigenpair \(\left( {\mu_{n} ,{\varvec{z}}_{{\varvec{n}}} } \right) \) written about a stationary state \({\varvec{u}}\left( {{\varvec{p}},\sigma } \right) \) as where \({\varvec{M}}\) is adjusted to perform either physical or eductive stability calculations. Equation (27) also was integrated forward in time by the Backward-Euler method [34] solving the nonlinear algebraic equations at each time step by Newton’s method. These results are equivalent to the forward/forward time integration used by others. 4 Results Solutions are represented either by their distribution functions \(W_{i} \left( x \right)\), cumulative distributions \(F_{i} \left( x \right) \equiv \mathop \smallint \limits_{0}^{x} W_{i} \left( {x^{\prime}} \right)dx^{\prime}\), the medians \(x_{med} \) of the distributions, \(F_{i} \left( {x_{med} } \right) = \frac{1}{2}\), or simply the differences \(\left( {W_{i} \left( 1 \right) - W_{i} \left( 0 \right)} \right)\). The relationship between the distribution functions and resource acquisition \(\left\{ {E_{i} } \right\}\) is or for \(h_{ai} \left( x \right) = 1\), 4.1 Case I: Competitive 2-Population MFG with Constant Resource A purely competitive 2-population MFG is modeled with constant resource level for both populations. The parameter values are \(q_{o1} = q_{o2} = 1, N = 2 \left( {{\text{fixed}}} \right),{\text{ and}} c_{o1} = c_{o2} = 1\), with \(\gamma_{21} = 4\), \(\gamma_{12} = 1\) and \(\gamma_{ii} = 0, i = 1,2\) so that the cost of acquiring the resource for population-2 is negatively impacted by population-1. For these parameters \(\sigma_{01} = 0.9488 \) and \(\sigma_{02} = 0.6709. \) A classic pitchfork bifurcation diagram is computed emanating from \(\sigma = \sigma_{01} , \) as shown in Fig. 2 plotted for \(\left( {W_{1} \left( 1 \right) - W_{1} \left( 0 \right)} \right)\) versus \(\sigma^{ - 1}\) . A typical solution is shown in Fig. 3 for \(\sigma = 0.851\). Samples distributions, cumulative distributions, and medians are shown in Fig. 3 corresponding to points on the bifurcation diagram shown in Fig. 1 for the \(n = 1\) and \(n = 2\) branches. For states on the first branch (n = 1), the relationship with the form \(\cos \left( {n\pi x} \right)\) is clear. Note that different measures are used to plot each branch. The separation of the medians is indicative of the segregation of the distributions. A more precise measure of segregation used by others [21, 22] is Separation of the medians does not require the distributions to satisfy Eq. (30). Sample solutions along the bifurcating families are shown in Fig. 4. The solutions 4A and 4B on the \(n = 1\) branch (B) demonstrate the segregation of the two populations which occurs over a narrow range of \(\sigma^{ - 1}\) above \(\sigma_{01}^{ - 1}\). Figure 4A shows that the segregation slows at higher \(\sigma^{ - 1}\). The inverted structure of the distributions on branch A is shown in Fig. 4C and is caused by the invariance of \(L_{i} \left( {x,W_{i} \left( x \right),N} \right)\) in \(x\). The solutions 4 K and 4L demonstrate the same inverted structure on the branches emanating from \(\sigma_{02}\) and the influence of the eigenfunction’s dependence on \(\cos \left( {2\pi x} \right)\). On each branch bifurcating from the trivial state, the constants \(\left\{ {\lambda_{Ui} } \right\} \) begin at the values given by Eq. (10) and increase with increasing \(\sigma^{ - 1}\). Because \(\lambda_{Ui} \ne 0\), the ergodic solutions do not correspond to stationary solutions of Eq. (3). The impact of introducing bias by setting \(h_{ai} \left( x \right) = x \) is shown by the perturbed solution branches in Fig. 3. One continuous branch extends over a large range of \(\sigma^{ - 1}\) and is the merger of the trivial solution with branch B for \(n = 1\); representative solutions are shown as Figs. 3G-J. With increasing \(\sigma^{ - 1}\) the distributions rapidly segregate with population-1 clearly dominating, driven by the bias to larger values of x and its advantage in the competition caused by \(\gamma_{21} > \gamma_{12}\). However, a value of \(\sigma^{ - 1}\) is reached where the distribution \(W_{2} \left( x \right)\) begins to shift to higher x (see Fig. 4I), forming a peak separated from the boundaries (Fig. 4J); the median of the distribution shifts with the peak. According to the definition Eq. (17), the distributions are segregated in this state as \(W_{1} \left( x \right)\) remains localized near \(x = 1\) so that \(W_{1} \left( x \right)W_{2} \left( x \right) \to 0\) everywhere with increasing \(\sigma^{ - 1}\). The rapid segregation with increasing \(\sigma^{ - 1} \) also is represented in Fig. 5 by the difference of the distribution medians. Without bias, this difference approaches an almost constant value for large \(\sigma^{ - 1} . \) With bias, a maximum in the separation of the medians appears and then decreases with increasing \(\sigma^{ - 1}\). 4.2 Case II: Competitive 2-Population MFG with Fixed Rate of Resource Replenishment Using Eq. (7) with \( R_{N} = 0\) to model the available resource accounts for decreasing available resource with a constant rate of replenishment \( f_{o}\). If the acquisition functions are independent of the distributions (\(\gamma_{ii} = 0)\) and without bias (\(h_{ai} = 1),\) \(E_{i} \equiv Nq_{oi}\) and is unchanged by varying \(\sigma\). Then the solutions correspond to those for \(N = f_{o}\) only with variations in the parameters \(\left\{ {\lambda_{Ui} } \right\}.\) Including self-enhancement in the acquisition function \(\left( {\gamma_{ii} = 1} \right)\) leads to variation of \(E_{i}\) with varying \(\sigma^{ - 1}\). The supercritical pitchfork bifurcation computed for this model is shown in Fig. 6 by the difference in the medians computed for \(h_{ai} = 1\) and \(h_{ai} = x;\) the parameters are the same as in Case I with \(f_{o}\) = 2 and \(\gamma_{ii} = 1\), resulting in \(N_{o} = 1/2\) and \(\sigma_{01} = 1.003.\) The result is similar to the calculation for \(\gamma_{ii} = 0\) (Fig. 5) with rapid changes in the medians for small variation in \(\sigma^{ - 1}\) leading to almost total segregation of the populations by \(\sigma^{ - 1} = 2\).5. Including bias to agents with more expertise (\(h_{ai} = x) \) again perturbed the bifurcation leading to a continuous family of states evolving to higher values of \(\sigma^{ - 1}\). The populations segregate over a small range of \(\sigma^{ - 1}\), influenced by the remnants of the bifurcation at \(\sigma_{o1}^{ - 1}\). Population-1 is favored with its majority near \(x = 1, \) whereas population-2 is segregated toward \(x = 0\), degrading resource acquisition. As in Case I, when \(\sigma^{ - 1}\) is decreased further, population-2 begins to shift toward higher x lowering the difference between the medians; see Fig. 6. The differences in the acquisition functions \(\left\{ {E_{i} } \right\}\) and \(N\) for Branch B with and without bias are shown in Fig. 7. Without bias, there is a relatively small difference between \(E_{1}\) and \(E_{2}\) as the dependence of distribution on x results only from the bifurcation; along Branch B, \(E_{2}\) increases from one at \(\sigma_{o1}^{ - 1}\) to approximately 1.123 at \(\sigma^{ - 1} = 1.135\) before decreasing toward one for larger \(\sigma^{ - 1}\). The decrease for increasing \(\sigma^{ - 1}\) is driven by the increasing segregation of the distributions which decreases \(N\) (see Fig. 7). With the constant replenishment rate \(f_{o} , \) \(E_{1} = f_{o} - E_{2}\). Without the introduction of bias, there is no distinction between the segregation of a population toward either \(x = 0\) or \(x = 1\); for branch B \(E_{2} > E_{1}\) and for branch A the values are reversed. The only differences between the acquisition functions \(\left\{ {E_{i} } \right\} \) are due to the difference between \(\gamma_{12}\) and \(\gamma_{21}\) which is confined to \(\sigma^{ - 1}\) near \(\sigma_{o1}^{ - 1}\). Adding bias (\(h_{ai} = x) \) breaks the invariance to \(x and \) skews the distributions to higher x. The bias interacts with the advantage to population-1 \((\gamma_{21} > \gamma_{12} )\) resulting in \(E_{1}\) becoming significantly larger than \(E_{2}\). The difference remains for higher values of \(\sigma^{ - 1}\) even as the peak in the distribution function \(W_{2} \left( x \right)\) drifts to larger x as shown in the insert in Fig. 5. 4.3 Case III: Cooperative-Competitive 2-Population Model and TOTC The finiteness of the resource is modeled by the logistic Eq. (7) without replenishment (\(f_{o} = 0)\). Calculations are presented for the parameters used above with \(R_{s} = 6\) and \(N_{eq} = 1\), giving \(N_{o} = 1/3\), \(\sigma_{01} = 0.9861\), and \(\sigma_{02} = 0.6973. \) The bifurcation diagram for branches A and B emanating from \(\sigma = \sigma_{01} \) is shown in Fig. 8 and is similar to Fig. 2 with constant N. However, the branches do not continue indefinitely with increasing \(\sigma^{ - 1}\). As the segregation increases with \(\sigma^{ - 1}\) both acquisition functions \(\left\{ {E_{i} } \right\}\) [Eq. (16)] and the resource N decrease to zero at \(\sigma^{ - 1} = \sigma_{limit}^{ - 1}\), setting the limit to each branch. This behavior is shown for Branch B in Fig. 9A. Calculations with \(h_{ia} \left( x \right) = x\) also are shown in Figs. 8 and 9B. Again, the populations rapidly segregate near the bifurcation point. Adding \(h_{ia} \left( x \right)\) enhances the acquisition of population-1 compared to population-2 and, for \(\sigma^{ - 1}\) near \(\sigma_{o1}^{ - 1} ,\) \(E_{1} \gg E_{2}\); \(E_{1}\) reaches a maximum and then decreases to zero at \(\sigma_{limit} = 0.0526\). At the peak in \(E_{1}\) (\(\sigma_{o1}^{ - 1} \cong 1.5) \) there is almost a six-fold difference between \(E_{1}\) and \(E_{2}\). The same behavior is shown for the evolution of \(\left\{ {E_{i} } \right\}\) and \(N\) for \(R_{N} = 3\) displayed in Fig. 10. Here there is a maximum four-fold difference between \( E_{1}\) and \( E_{2} \) before \(\sigma_{limit} = 0.3999 {\text{is reached}}.\) The limit \(\sigma_{limit}\) defines the TOTC where the resource N has been exhausted. The dependence of \(\sigma_{limit}\) on resource production \(R_{N}\) is displayed in Fig. 11 and demonstrates two regimes. First, near the lower limit \(R_{N} = R_{N}^{crit}\) for \(\left( {\sigma_{o1} - \sigma } \right) \ll 1\), where \(R_{N}^{crit}\) was set by the distributions evolving back to the trivial state \(W_{1} \left( x \right) = W_{2} \left( x \right) = 1\) as N decreased. For higher values of \(R_{N} ,\) \(R_{N} = R_{N}^{crit} \) is determined by the increasing segregation of the distribution functions. The curve \(R_{N} = R_{N}^{crit} \left( {\sigma^{ - 1} } \right) \) define the minimum value of \(R_{N}\) where \(N = 0\). Increasing \(\sigma^{ - 1}\) increases the effectiveness of each population’s acquisition, leading to more resource consumption that can only be accommodated by increasing resource production. The rate of increase of \(R_{N}\) was greatest for \(\sigma^{ - 1}\) near \(\sigma_{o1}^{ - 1}\) as the two populations segregate, driven by the residual of the bifurcating branches. 4.4 Effect of Nonlinear Cost Function The calculations presented above are based on the utility functions for cost and acquisition depending linearly on the distributions \(W_{i} \left( x \right)\). Following [22] calculations were performed with quadratic dependence of the cost function, \(\beta_{ci} = 2\). The bifurcating families corresponding to \(n = 1\) and \(\delta_{ai} = 0\) are shown in Fig. 8 for comparison to the results for \(\beta_{ci} = 1\); the bifurcation point increases to \(\sigma_{o1}^{ - 1} = 0.8991 \) according to Eq. (17b). The evolution of the medians of the distributions are similar for both values of \(\beta_{ci}\). For \(\beta_{ci} > 1\), the shift in \(\sigma_{o1}^{ - 1}\) to higher values leads to the TOTC \(\sigma_{limit} \) occurring at higher values of \(\sigma^{ - 1} .\) Calculations with \(R_{N} = 3\) and \(1 \le \beta_{ci} \le 3\) predict \(\sigma_{limit}\) varying as \(0.3990 \le \sigma_{limit} \le 0.58642\), with \(0.9861 \le \sigma_{o1} { } \le 1.2658\). Not surprisingly, the distributions at \(\sigma^{ - 1} = \sigma_{limit}^{ - 1}\) are similar, independent of \(\beta_{ci}\). Interestingly, as implied by the analysis by Cirant and Verzini [22] including cubic dependences in the cost function \((\beta_{ci} = 3)\) yields a subcritical bifurcation at \(\sigma^{ - 1} = \sigma_{o1}^{ - 1}\) where the branches begin by evolving to lower values \(\sigma^{ - 1} . \) Numerically tracking the A & B branches shows a very small range, \(\sigma^{ - 1} \le \sigma_{o1}^{ - 1}\), before the branches turn back to increasing \(\sigma^{ - 1}\). This small range of \(\sigma^{ - 1}\) defines a hysteresis loop making the segregation of the distributions appear almost “instantaneous” as \(\sigma^{ - 1}\) is increased from \(\sigma_{o1}^{ - 1}\). 5 On the Stability of the Ergodic States The multiple ergodic states raise questions about their relative stability. Eigenvalue calculations of the physical stability for states along each branch reproduce the pattern predicted for the trivial state; the supercritical bifurcating branch emanating from \(\sigma_{o1}^{ - 1} \) has all positive eigenvalues (stable), whereas the branch emanating from \(\sigma_{o2}\) has a single negative eigenvalue (unstable), corresponding to the eigen-structure of the trivial state for \(\sigma^{ - 1} > \sigma_{o1}^{ - 1}\). The physical stability of separated bifurcation branches caused by introducing bias follows the classic theory for a simple bifurcation with the states along both the continuous family emanating from small values of \(\sigma^{ - 1}\) and the branch formed from branch A being stable. As a result, there is a range of \(\sigma^{ - 1} \) where two physically stable states exist. Repeating the calculations with \(M \) modified for eductive stability gives the same stability predictions. 6 Discussion The popularity of Mean Field Game Theory has been growing rapidly, and two-population games are being deployed in many applications. A simple model for two-population model for competition for a constant or potentially depleting resource is analyzed to explore the role of nonlinear interactions on the behavior of the distributions in the ergodic limit. Even in this simple formulation—linear coupling between populations and quadratic control—the computations demonstrate the potential for nonlinear dynamics to dramatically influence the response of the system. For “perfect systems” without explicit dependence on the independent variable x and with three different models for consumption of the common resource, results demonstrate supercritical bifurcations in \(\sigma^{ - 1}\). In all cases there is rapid segregation of the two distributions over small ranges of \(\sigma^{ - 1}\) driven by the competition between the populations, is demonstrated by the distribution medians seen in Figs. 5, 6, And 8. The effects of the multiple states are visible in MFG models perturbed away from the perfect states as demonstrated by introducing bias in the resource acquisition that favors agents at large x, modelling agents with more expert or effort. These states exhibit the rapid segregation with increasing \(\sigma^{ - 1}\). In addition, for still higher values of \(\sigma^{ - 1}\), nonlinear couplings result in drift of the population segregated from lower x (lower expertise) to higher x, decreasing the difference between the medians of the distributions, although the distributions remain segregated. The predictions for resource acquisition by the two populations is affected by the nonlinear dynamics. When the resource is regenerated at a fixed rate but without bias, the acquisition differences are small. Without explicit dependence on x of the acquisition function, the segregation of the populations caused only small differences in resource acquisition \(\left\{ {E_{i} } \right\}\). Introducing bias with \(h_{ai} \left( x \right) = x\) leads to strong differentiation between the populations and large differences in \(\left\{ {E_{i} } \right\}\). These differences become pronounced when the resource is self-regenerating so that over acquisition leads to its depletion. With the bias introduced the rapid segregation caused by the remnants of the bifurcation leads to maxima in resource acquisition \(\left\{ {E_{i} } \right\}\) followed by continuous decline until \(N = 0\) is reached. Increasing the control (increasing \(\sigma^{ - 1}\)) requires increasing rates of production by the resource to avoid the TOTC. For a natural resource (for example fish populations) that was in equilibrium before agents began to acquire it, the limiting value \(\sigma_{limit}^{ - 1}\) changes slowly with increasing reproduction rate \(R_{N}\), making the TOTC highly sensitive to changes in effectiveness measured by \(\sigma\). Mitigation strategies are required to avoid the TOTC limit in this region of parameter space. Models for such strategies can be introduced in either the acquisition (6a) or cost (6b) functions as dependencies on either the size of the resource population N or the acquisition rates \(\left\{ {E_{i} } \right\}\). Analysis of these effects is left to more detailed MFG models developed to explore specific resources such as those in [17, 26, 28]. Data Availability No datasets were generated or analysed during the current study. References Hardin G (1968) The tragedy of the commons. Science 162:1243–1248 Ostrom E (1990) Governing the commons: the evolution of institutions for collective action. Cambridge University Press Ostrom E (1999) Coping with tragedies of the common. Annu Rev Polit Sci 2(1):493–535 Lam KY, Liu S, Lou Y (2020) Selected topics on reaction-diffusion-advection models from spatial ecology. Math Appl Sciences Engng 1(2):150–180 Lasry J-M, Lions P-L (2006) Jeux à champ moyen. II–Horizon fini et contrôle optimal. Comptes Rendus Mathématique 343(10):679–684 Lasry J-M, Lions P-L (2006) Jeux à champ moyen. I–le cas stationnaire. Comptes Rendus Mathématique. 343(9):619–625 Lasry J-M, Lions P-L (2007) Mean field games. Japan J Math 2(1):229–260 Huang M, Malhamé RP, Caines PE (2006) Large population stochastic dynamic games: closed-loop McKean-Vlasov systems and the Nash certainty principle. Commun InF Syst 6(3):221–252 Cardaliaguet P, Lasry J-M, Lions P-L, Porretta A (2012) Long time average of mean field games. Netw Heterog Media 7:279–301 Cardaliaguet P, Lasry J-M, Lions P-L, Porretta A (2013) Long time average of mean field games with a nonlocal coupling. SIAM J Control Optim 51:3558–3591 Bensoussan A, Frehse J, Yam P (2013) Mean field games and mean field type control theory. Springer Cardaliaguet P, Porretta A (2020) An introduction to mean field game theory. In: Cardaliaguet P, Porretta A (eds) Mean field games. Springer Carmona R, Delarue F (2018) Probabilistic theory of mean field games with applications I. Springer Dayaniki G, Laurière M (2025) Cooperation, competition, and common pooled resources in mean field games. arXiv:2504.09043v1 Kobeissi Z, Mazari-Fouquer I, Ruiz-Balet D (2024) Mean-field games for harvesting problems: Uniqueness, long-time behaviour and weak kam theory. arXiv:2406.06057 Kobeissi Z, Mazari-Fouquer I, Ruiz-Balet D (2024) The tragedy of the commons: a mean-field game approach to the reversal of travelling waves. Nonlinearity 37(11):115010 Yoshioka H, Tsujimura M, Yoshioka Y (2024) Numerical analysis of an extended mean field game for harvesting common fishery resource. Comput Math Appl 165:88–105 Bensoussan A, Huang T, Lauriére M (2018) Mean field control and mean field type game models with several populations. Minimax Theory Appl 3:173209 Cirant M (2015) Multi-population mean field game systems with Neumann boundary conditions. J Math Pures Appl 103:1294–1315 Feleqi E (2013) The derivation of ergodic mean field game equations for several populations of players. Dyn Games Appl 3(4):523–536 Achdou Y, Bardi M, Cirant M (2017) Mean filed games models of segregation. Math Model Methods Appl Sci 27(1):75–113 Cirant M, Verzini G (2017) Bifurcation and Segregation in quadratic two-population mean field game systems. ESAIM Control Optim Calc Var 23(3):1145–1177 Guéant O, Lasry J-M, Lions P-L (2010) Mean field games and applications. In: Caromona RA, Cinlar E, Ekeland I, Jouini E, Scheinkma J, Tousi N (eds) Paris-princeton lectures on mathematical finance 2010. Springer, pp 205–266 Guéant O (2009) A reference case for mean field games. J Math Pure Appl 92:276–294 Iooss G, Joseph DD (1980) Elementary stability and bifurcation theory. Springer Brown RA (2026) Mean field game model of the impact of reductions in support on faculty research activity. Proc Nat Acad Sci 123(23):e2538029123 Ullmo D, Swiecicki I, Gobron T (2019) Quadratic mean field games. Phys Rep 799:1–35 McPike R (2021) Integrating ecology and economics in the mathematical modelling of marine ecosystems. University to Strathclyde, Glasgow Ungar LH, Brown RA (1982) The dependence of the shape and stability of captive rotating drops on multiple parameters. Phil Trans R Soc Lond A 306:347–370 Kuznetsov YA (1998) Elements of applied bifurcation theory, 2nd edn. Springer MathWorks Inc (2023) MATLAB version: R2023b Update 6. The MathWorks Inc Chan T (1984) Newton-like pseudo-arclength methods for computing simple turning points. SIAM J Sci Stat Comput 5(1):135–148 Keller HB (1977) Numerical solution of bifurcation and nonlinear eigenvalue problems. In: Rabinowitz P (ed) Applications of bifurcation theory. Academic Press, New York, pp 359–384 Ascher UM, Petzold LR (1998) Computer methods for ordinary differential equations and differential-algebraic equations. SIAM Author information Authors and Affiliations Contributions RAB is solely responsible for the conseptualization, execution, and writing of this manuascript. Corresponding author Ethics declarations Conflict of interest The authors declare no competing interests. Additional information Publisher's Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations. Rights and permissions Open Access This article is licensed under a Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International License, which permits any non-commercial use, sharing, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if you modified the licensed material. You do not have permission under this licence to share adapted material derived from this article or parts of it. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by-nc-nd/4.0/. About this article Cite this article Brown, R.A. Nonlinear Dynamics in 2-Population Mean Field Games Competing for a Resource: Influence on the Tragedy of the Commons. Dyn Games Appl (2026). https://doi.org/10.1007/s13235-026-00711-4 Received: Accepted: Published: Version of record: DOI: https://doi.org/10.1007/s13235-026-00711-4