Why multilingual, and how to keep it—An evolutionary dynamics perspective

While many languages are in danger of extinction worldwide, multilingualism is being adopted for communication among different language groups, and is playing a unique role in preserving language and cultural diversities. How multilingualism is developed and maintained therefore becomes an important interdisciplinary research subject for understanding complex social changes of modern-day societies. In this paper, a mixed population of multilingual speakers and bilingual speakers in particular is considered, with multilingual defined broadly as zero, limited, or full uses of multiple languages or dialects, and an evolutionary dynamic model for its development and evolution is proposed. The model consists of two different parts, formulated as two different evolutionary games, respectively. The first part accounts for the selection of languages based on the competition for population and social or economic preferences. The second part relates to circumstances when the selection of languages is altered, for better or worse, by forces other than competition such as public policies, education, or family influences. By combining competition with intervention, the paper shows how multilingualism may evolve under these two different sources of influences. It shows in particular that by choosing appropriate interventional strategies, the stable co-existence of languages, especially in multilingual forms, is possible, and extinction can be prevented. This is in contrast with major predictions from previous studies that the co-existence of languages is unstable in general, and one language will eventually dominate while all others will become extinct.


Introduction
As the world becomes increasingly globalized, more people become multilingual across the continents. It is reported that now more than half of the world's population speak at least two languages [1]. While many languages are in danger of extinction, multilingualism is being adopted as a common way of communication among different language groups, and is playing a unique role in preserving language and cultural diversities [2]. Therefore, how multilingualism is developed and maintained becomes an important research subject in linguistics, social studies, and beyond [3]. The promotion and protection of multilingualism has also been a hot topic long discussed in public [4]. populations under different conditions are analyzed. The stability conditions of equilibrium states are also justified for different interventional strategies. In Dynamic Simulation, a twodimensional dynamic simulation scheme based on the proposed evolutionary model is described. Two sets of simulation results are presented showing potential geographical impacts on the evolution of multilingual populations. To keep it simple and easy to follow, the paper is focused more on bilingual populations, with models for general multilingual populations provided as an appendix in the end. To make the paper more accessible, all formal definitions, theorems, and especially proofs are also organized into the appendix sections. Nonetheless, these appendices are considered to be necessary contents as well for supporting the results and conclusions presented in the paper.

Evolutionary models
To keep the description simple and easy to follow, consider a bilingual population, the simplest yet most common form of multilingualism. Assume it is large and well mixed, i.e., every individual speaker can interact with all others in the population. If an individual speaker uses two languages A and B with frequencies x A and x B , respectively, this individual is called an (x A , x B )speaker, where 0 � x A , x B � 1, and x A + x B = 1. Thus, a (1, 0)-speaker is an A speaker, a (0, 1)speaker is a B speaker, and an (x A , x B )-speaker with x A , x B 6 ¼ 0 is a speaker of both A and B, with frequency x A for A and x B for B. Likewise, if language A and B are used with frequencies y A and y B in average in the whole population, this population is called a (y A , y B )-population, where 0 � y A , y B � 1, and y A + y B = 1. Thus, a (1, 0)-population is an A population, a (0, 1)-population is a B population, and a (y A , y B )-population with y A , y B 6 ¼ 0 is a bilingual population of A and B, with average frequency y A for A and y B for B. Here, if the population size is n and speaker i uses A and B with frequencies x ðiÞ A and x ðiÞ B , respectively, then y A ¼ S i x ðiÞ A =n and y B ¼ S i x ðiÞ B =n.

Evolution with competition
Let P A (y A ) and P B (y B ) be the payoff functions for A and B speakers in a (y A , y B )-population, respectively, defined in terms of population percentages and social or economic preferences, with P A increasing in y A and P B in y B , meaning that the larger the population percentage of a language, the more benefit the language provides. Then, the payoff function for a general (x A , x B )-speaker in a (y A , y B )-population can be defined in terms of the average use of A and B by this speaker: For example, the payoff of a (0.5, 0.5)-speaker, i.e., a half-A and half-B speaker, in a (0.5, 0.5)population, i.e., a half-A and half-B population, should be 0.5P A (0.5) + 0.5P B (0.5).
Considering (x A , x B ) as the strategy of an individual and (y A , y B ) the strategy of the population, an evolutionary game (Appendix A) can be defined, where every individual tries to maximize his/her payoff. The latter can be achieved when an optimal strategy ðx � A ; x � B Þ is found for every individual. The strategy for the whole population then becomes ðx � A ; x � B Þ as well, and a Nash equilibrium is reached. A Nash equilibrium of this game is thus a strategy ðx � A ; x � B Þ such that Assume for a (y A , y B )-population that y A and y B vary over time as they reach their equilibrium. Then, a system of replicator equations (Appendix A) can also be defined for the population: meaning that the changing rate of average frequency y A (or y B ) of language A (or B) depends on the payoff P A (or P B ) compared with the payoff P B (or P A ): if it is higher, y A (or y B ) increases; otherwise, it decreases. Based on evolutionary game theory, a Nash equilibrium of the evolutionary game as conditioned in (2) is always a fixed point of the system of replicator equations given in (3) [18][19][20] (Appendix A).

Evolution with intervention
Now consider the situation where some societal influences other than simple competition intervene in the use of languages. Assume that the societal influences or interventions are implemented to counter the arbitrary increase or decrease of either language in a bilingual population. Let � P A ðy A Þ and � P B ðy B Þ be the payoff functions for A and B speakers in a (y A , y B )population, respectively, defined in terms of societal influences and interventions, with � P A decreasing in y A and � P B in y B , meaning that the smaller the population percentage of a language, the more incentive or less penalty the speakers of the language receive. Here, incentive means more benefit from the society such as more work opportunities or public supports, etc., and penalty means less advantage in the society such as less access to public education or discouragement from the public, etc. Then, the payoff function for a general (x A , x B )-speaker in a (y A , y B )-population can be defined in terms of the average use of A and B by this speaker: It follows that another evolutionary game can be defined, with the Nash equilibrium being a strategy ðx � and also a corresponding system of replicator equations: Note that the game with competition in (2) has three possible equilibrium strategies. They can be obtained as the fixed points of the system of replicator equations in (3): It is not so hard to prove that the strategies in (a) and (b) are evolutionarily stable while the one in (c) is not (see Appendix C). This means that under competition, one of the languages is expected to die eventually while the other takes over the whole population. The co-existence of the two languages is not sustainable.
On the other hand, the game with intervention in (5) also has three possible equilibrium strategies. They can be obtained as the fixed points of the system of replicator equations in (6): However in this case, it can be proved that the strategy in (c) is evolutionarily stable while the ones in (a) and (b) are not (see Appendix C). This means that with intervention, it is possible to keep the population of either language not too large or too small. A bilingual population can be maintained.

Evolution with competition and intervention
Now the two types of games described above can be combined to obtain a complete model for the evolution of a bilingual population. A simple combination can be obtained by merging the systems in (3) and (6) into the following system: where λ 2 [0, 1] is a scalar. If λ = 1, the system is reduced to the one in (3) with competition only. If λ = 0, it is reduced to the one in (6) with intervention only. If λ 2 (0, 1), the system is in some sense a convex combination of competition and intervention, with more weight on competition when λ > 0.5 while more on intervention when λ < 0.5. For convenience, however, a simpler system, with an equal weight given to competition and intervention, will be considered in the following discussion: In any case, with a combined system, the impacts from language competition and influences of societal intervention are both included for the uses of the languages. The first part of the system, based on language competition, pushes the population to monolingual, while the second part, counted for societal interventions, prevents the population from losing either language. The combination of the two is then expected to provide a more comprehensive account on the changes of the languages and the evolution of the bilingualism. Note that the system in (10) can be reformulated in a more compact form: It is then a system of replicator equations with payoff functionsP A andP B , and corresponds to an evolutionary game, with the Nash equilibrium being a strategy ðx � wherep is the payoff function for the game, and for an (x A , x B )-speaker in a (y A , y B )-populationp In general, P A and P B can be any increasing functions and � P A and � P B be any decreasing functions. However, based on previous studies, they can be defined as the following: where c and � c are scaling constants, α, � a, a determine the order of dependency of the payoffs on the population percentages. Since 1 < α � 2, the payoffs from P A and P B increase with increasing population percentages. On the other hand, since 0 � � a < 1, the payoffs from � P A and � P B decrease with increasing population percentages. The parameters s A , s B are used to define the payoffs from competition. They are indicators of social or economic impacts on the payoffs. The larger these values, the more benefit for the corresponding language groups. The parameters � s A ; � s B are used to define the payoffs from interventions. They are rates for language reversing due to interventions. The larger these values, the faster the reversing rates.
The definitions of the functions in (16) and (17) are actually based on previous work on language competition and, in particular, the work by Abrams and Strogatz 2003 [13]. In fact, with these functions, the systems in (3) and (6) are both equivalent to an Abrams-Strogatz system as shown below, differing only in the ranges of the α values and the parameters c, s A , s B versus � c, � s A , � s B : where P is a function, P(y, s) = cy α s for some constant c. Abrams and Strogatz 2003 surveyed a large number of regions in UK and South America for language competition and estimated the α value around 1.31 ± 0.25 [13]. In general, the system can be adopted for modeling language competition for 1 < α � 2 as well as language reversing for 0 � α < 1, as discussed above.

Dynamic analysis
Continue the assumption of bilingual populations and consider their possible dynamic behaviors based on the models described in previous sections. Assume that the payoffs P A ; P B ; � P A ; � P B are defined by the functions in (16) and (17).

Dynamics with competition
For competition only, an evolutionary game and a corresponding system of replicator equations are given in (2) and (3). With P A and P B given in (16), the game in (2) is actually a so-called potential game [18,20] (Appendix B), i.e., there is a potential function f such that Indeed, f(y A , y B ) = (P(y A , s A ) + P(y B , s B ))/α, where P(y, s) = cy α s, 1 < α � 2. Therefore, an equilibrium strategy for the game must be a KKT point of the following potential maximization problem: Further, if the equilibrium strategy is evolutionarily stable, it must be a strict local maximizer of the potential maximization problem, and vice versa [18,20] (Appendix B). For a given strategy (y A , y B ), the Hessian H of f is given as follows.
According to the theory on constrained optimization [21,22], as a KKT point of the potential maximization problem in (20), if an equilibrium strategy for the game in (2) is a local maximizer, the Hessian projected on the null space of the active constraints at this strategy is necessarily negative semidefinite; conversely, if the projected Hessian is negative definite, the equilibrium strategy must be a strict local maximizer of the potential maximization problem [21,22]. Thus, the evolutionary stability of the equilibrium strategy can be justified by the negative definiteness of the projected Hessian at the strategy (Appendix B). The game in (2) has three possible equilibrium strategies corresponding to three fixed points of (3), as given in (7). It is easy to see that strategies (a) and (b) are evolutionarily stable (Appendix C). In either case, one language dies while the other takes over the whole population. If s A > s B , the chance for (a) will be greater than (b), and if s B > s A , the chance for (b) will be greater than (a). For (c), it is easy to verify that the projected Hessian at this strategy is always positive definite (Appendix C). It follows that the strategy can never be a local maximizer of the potential maximization problem in (20), and can never be evolutionarily stable.
The above behaviors of the game can be further demonstrated in Fig 1, where P A and P B are plotted against y A , and the changing directions of y A are pointed with arrows. The graph in Fig  1(a) shows the behavior of a population with α = 2, when P A and P B are simple linear functions. The curves of P A and P B are displayed over y A 2 [0, 1] with s A = 0.75 and s B = 0.25, when language A dominants B. There is an equilibrium strategy y � A in (0, 1) such that P A ðy � A Þ ¼ P B ðy � B Þ, with y � A ¼ 0:25 and y � B ¼ 0:75. It is unstable as y A decreases to 0 if it starts from below y � A while increases to 1 if it starts from above y � A . Note that y � A is relatively small. Therefore, y A increases to 1 even if y A starts small.
The graph in Fig 1(b) shows the behavior of a population with α = 3/2, when P A and P B are square root functions. The curves of P A and P B are displayed over y A 2 [0, 1] with s A = 0.25 and s B = 0.75, when language B dominants A. There is an equilibrium strategy y � A in (0, 1) such that P A ðy � A Þ ¼ P B ðy � B Þ, with y � A ¼ 0:9 and y � B ¼ 0:1. It is unstable as y A decreases to 0 if it starts from below y � A while increases to 1 if it starts from above y � A . Note that y � A in this case is very large. Therefore, y A decreases to 0 even if y A starts large.

Dynamics with intervention
For intervention only, an evolutionary game and a corresponding system of replicator equations are given in (5) and (6). With � P A and � P B given in (17), the game in (5) is also a potential game, i.e., there is a potential function � f such that Therefore, an equilibrium strategy for the game must be a KKT point of the following potential maximization problem: Further, if the equilibrium strategy is evolutionarily stable, it must be a strict local maximizer of the potential maximization problem, and vice versa [18,20] (Appendix B). For a given strategy (y A , y B ), the Hessian � H of � f is, Again, as a KKT point of the potential maximization problem in (23), if an equilibrium strategy for the game in (5) is a local maximizer, the Hessian projected on the null space of the active constraints at this strategy is necessarily negative semidefinite; conversely, if the projected Hessian is negative definite, the equilibrium strategy must be a strict local maximizer of the potential maximization problem [21,22]. Thus, the evolutionary stability of the equilibrium strategy can be justified by the negative definiteness of the projected Hessian at the strategy (Appendix B). The game in (5) has three possible equilibrium strategies corresponding to three fixed points of (6), as given in (8). It is easy to see that strategies (a) and (b) are unstable (Appendix C). For (c), it is easy to verify that the projected Hessian at this strategy is always negative definite (Appendix C). It follows that this strategy must be a strict local maximizer of the potential maximization problem in (23), and must be evolutionarily stable.
The above behaviors of the game can be further demonstrated in Fig 2, where � P A and � P B are plotted against y A , and the changing directions of y A are pointed with arrows. The graphs in when language A is more promoted than B. There is an equilibrium strategy

Dynamics with competition and intervention
Now consider the dynamic behaviors of a bilingual population under the influences of both inter-language competition and external interventions, as modeled by the system in (10) and the corresponding game in (14). The payoff functionsP A andP B are simply combinations of those for competition and intervention as given in (12) and (13). With these two payoff functions, the game in (14) is again a potential game, and therefore, there must be a functionf such that with f and � f defined in previous sections. It follows that an equilibrium strategy for the game must be a KKT point of the potential maximization problem: Further, if the equilibrium strategy is evolutionarily stable, it must be a strict local maximizer of the potential maximization problem, and vice versa [18,20] (Appendix B). For a given strategy (y A , y B ), the HessianH off is  Again, as discussed in the previous two sections, if an equilibrium strategy for the game in (14) is a local maximizer of the potential maximization problem in (26), the Hessian projected on the null space of the active constraints at this strategy is necessarily negative semidefinite; conversely, if the projected Hessian is negative definite, the equilibrium strategy must be a strict local maximizer of the potential maximization problem [21,22]. Thus, the evolutionary stability of the equilibrium strategy can be justified by the negative definiteness of the projected Hessian at the strategy (Appendix B).
For convenience, set c ¼ � c ¼ 1 for the following analysis. Let ðy � A ; y � B Þ, y � A ; y � B 6 ¼ 0, be an equilibrium strategy for the evolutionary game in (14). Then, This is a necessary and sufficient condition for ðy � A ; y � B Þ, y � A ; y � B 6 ¼ 0, to be an equilibrium strategy for the evolutionary game in (14). For a specific population, for example, for α = 3/2 and � a ¼ 0, it can be simplified to and for α = 3/2 and � a ¼ 1=2, to Note that at ðy � (26) has only one active constraint y A + y B = 1, and the Hessian projected on the null space of this constraint will beĤ This is a sufficient condition for an equilibrium strategy ðy � A ; y � B Þ, y � A ; y � B 6 ¼ 0, to be evolutionarily stable. For a specific population, for example, for α = 3/2 and � a ¼ 0, the condition can be simplified to: and for α = 3/2 and � a ¼ 1=2, to An equilibrium strategy ðy � A ; y � B Þ, y � A ; y � B 6 ¼ 0, corresponds to a bilingual population with both language A and B co-existing. In previous sections, competition-only and interventiononly populations have been discussed. For the case with competition-only, such a strategy is never stable, and therefore, the two languages can never co-exist, even in bilingual forms. For the case with intervention-only, such a strategy is always stable, and the two languages can coexist. By combining the two, it is hoped that such an equilibrium strategy still exists and is also stable. Indeed, such an equilibrium can be achieved by choosing appropriate interventional strategies: , and therefore, the stability condition in (33) With such a choice of � s A and � s B , one can prove that ðy � A ; y � B Þ is also unique (see Appendix D for verification and Fig 3(a) for demonstration).
is an equilibrium strategy for all three type of games, the competition-only game, the intervention-only game, and the combination of the two, i.e., the game with both competition and intervention. Mathematically, In addition, the stability condition in (33) is satisfied for ðy � A ; y � B Þ. One can then prove that ðy � A ; y � B Þ is evolutionarily stable and also unique (see Appendix D for verification and Fig 3(b) for demonstration).
Strategy 3: In general, given a desired equilibrium strategy ðy � A ; y � B Þ, y � A ; y � B 6 ¼ 0, not necessarily optimal for either the competition-only or intervention-only game, it can be made a true equilibrium strategy for the combined game by choosing appropriate interventional strategies: Then, the condition in (33) is satisfied at ðy � A ; y � B Þ, and the strategy is also evolutionarily stable (see Appendix D for verification and Fig 4 for demonstration).
Note that in the last case, when ða À 1Þ=ð1 À � aÞ is close to 1, the functionsP A andP B are likely to have more than one intersections over the interval (0, 1) for y A , as demonstrated in the examples in Fig 4. This implies that there may be more than one equilibrium strategies, among which is the desired one. In this case, the interval (0, 1) is divided by the y � A values of these strategies, and the y � A value of the selected strategy is converged only within a small subinterval around it. When ða À 1Þ=ð1 À � aÞ is much smaller than 1, however, the functions P A andP B are more separated, and usually there is only a single intersection, corresponding to a unique equilibrium strategy.
The dynamic behaviors of the combined game in (14) can be demonstrated through examples shown in Figs 3 and 4. The graph in Fig 3(a) shows the behavior of a population with α = 3/2 and � a ¼ 0.  Fig 3(b) shows the behavior of a population with α = 3/2 and � a ¼ 1=2.   half circles. The desired strategy is one of them. This happens when ða À 1Þ=ð1 À � aÞ is close or equal to 1.

Dynamic simulation
Computer simulation can always be used to validate theories when physical experiments are not accessible, and to reveal additional insights into system behaviors which may otherwise be hard to see. This section describes further dynamic behaviors of bilingual populations through computer simulation with the combined model given in (10) and (14). To carry out the simulation, the population is assumed to be distributed in a two-dimentional space, and the changes of the bilingual level of the population will then be tracked across time and space.
More specifically, to carry out the simulation, a 2D torus-shaped lattice of n × n cells is constructed first, with each cell assumed to be occupied by an individual speaker. A random individual can then be selected repeatedly from the lattice, and a game is played for the individual against the population of the lattice. Let (x A , x B ) be the current strategy for the individual, and (y A , y B ) the strategy for the population. Let p A ¼P A ðy A Þ and p B ¼P B ðy B Þ be the payoffs for A and B speakers, respectively. Then, the payoff for the individual, Initially, all individuals are assigned with a random strategy. The game is played n × n times for the population to complete a generation. The game is repeated for 100 generations to make sure the population reaches its equilibrium. In general, the game can be played in a neighborhood of each selected individual. Let the neighborhood be an m × m sub-lattice, with the selected individual located at the center. Then, the game can be carried out for each selected individual only against the population in its neighborhood of this size, with the population strategy (y A , y B ) computed from the population in the neighborhood. Such a game may in fact be more realistic, as people usually interact only with a small group of others around them.
Recall that the population is assumed to be large and well mixed in the game models discussed in previous sections. This means that every individual should be able to interact with all others in the population. With such an assumption, the dynamic behaviors and especially the equilibrium states of the population can be computed by directly solving a corresponding system of replicator equations. The 2D simulation described here can be considered as a discrete solution of the system of replicator equations if the individual in each cell of the lattice is allowed to interact with the individuals in all other cells. However, the 2D simulation is more flexible, allowing restrictions on interactions. It can therefore be used to see how a population behaves when each individual interacts only with a selected group of other individuals, which a continuous model would not be able to reveal. In the first column of graphs, the one on the top shows the result from the game with the neighborhood size equal to 75 × 75, when each individual interacts with all others in the whole population. The equilibrium frequency y � A of the population in this case is approximately equal to 0.75, which agrees with the direct prediction from the continuous model described in previous sections. In addition, the distribution of the individual frequency x � A in the population is very homogeneous, with x � A � 0:75 across the board, suggesting that language A and B coexist in the population in an evenly distributed bilingual form. The second graph in this column shows the result from the game with the neighborhood size equal to 25 × 25, when the interactions among individual speakers are restricted. The equilibrium frequency y � A of the population remains about the same, approximately equal to 0.75. However, the individual frequency x � A becomes less constant. Some regions have higher individual frequencies than others, and local groups are formed with varying individual frequencies, as shown in the graph. The graph in the bottom of this column shows the result from the game with the neighborhood size further reduced to 5 × 5. While the population frequency y � A is not significantly changed, the individual frequency x � A shows even bigger variations, with even smaller local spots formed with higher or lower individual frequencies than average. The graphs in the second column show similar results. Again, when the neighborhood size is set to 75 × 75, the equilibrium frequency y � A of the population reaches approximately 0.4, which agrees with the direct prediction from the continuous model described in previous sections. The distribution of the individual frequency x � A in the population is homogeneous, with x � A � 0:4 across the board, as shown in the graph on the top. However, when the neighborhood size is reduced to 25 × 25, the individual frequency x � A becomes less constant, although y � A remains about the same. Some regions have higher individual frequencies than others, and local groups are formed with varying individual frequencies, as shown in the second graph. When the neighborhood size is further reduced to 5 × 5, the individual frequency x � A varies in greater degrees in the population. Even smaller local spots are formed with higher or lower individual frequencies than average. The dynamic behaviors shown from these simulations agree with our intuition or experience in language development: Indeed, when communications are restricted to local groups, language variations will remain.

Discussion
The competition-only model as described in the Evolution with Competition section predicts that between two competing languages, if based only on competition, one would eventually die while the other takes over the whole population, and the co-existing state is unstable. This result is not so surprising, for it has already been discussed in previous studies, although from different perspectives. On the other hand, the intervention-only model as described in the Evolution with Intervention section shows that if controlled only by interventions, it is possible to prevent either language population from becoming too large or too small, and keep the population in a stable co-existing state. By combining the two, the model with both competition and intervention as described in the Evolution with Competition and Intervention section gives a more complete description on the evolution of multilingual populations when it is under the influences of both language competition and societal intervention. It predicts that languages may co-exist stably in multilingual forms if appropriate interventional strategies are employed. In addition, the interventional measures may not only be able to prevent language extinction but also direct populations to desired equilibrium states.
These predictions are based on the assumptions that the payoff functions P A and P B are increasing functions and � P A and � P B are decreasing functions, meaning that in some sense, competition favors large populations, while intervention subsidies small populations. For this reason, the functions in (16) and (17) are adopted for defining P A , P B , � P A , and � P B . Besides, P A and P B in (16) can be considered as a general set of functions that includes the one used in the Abrams-Strogatz model with α = 1.31 ± 0.25, which is well tested against real language data [13]. The functions for � P A and � P B in (17) are also related to some applications in genetic selection, where they are used to model genetic mutations with � a ¼ 0 [6,19]. However, the � s A and � s B values for genetic mutation are certainly in a much smaller scale than those for language conversion.
Note that this paper has not discussed in detail how the language is acquired, learned, and used, which may in fact relate directly to what language strategies really mean and how they can be changed. There are some obvious questions related to this issue. For example, what does it mean that an individual uses language A with 40% frequency? Is only reading language A without speaking counted as using language A? What happens if an individual can never learn a second language? How long does it take for an individual to change his/her language use from one frequency to another? These are legitimate concerns, but may require more linguistic characterizations on language learning and communication, and will hopefully be addressed in future work.
The work in this paper is to develop a general framework for modeling and simulating evolutionary dynamics of multilingual populations. Work remains to be done to analyze some real language groups, while several issues are yet to be resolved: First, some terms used in the proposed models need to be further discussed for exact implications such as "social or economic preferences", which determine the parameters s A and s B , and "public policies", "educational influences", and "family influences", which the parameters � s A and � s B depend on. Second, even if exact meanings of these terms are given, the corresponding parameters including α and � a are not easy to estimate. They may vary with varying populations and varying time periods. The language data is also difficult to collect across a meaningful period of time. Third, the factors that affect language evolution can be more than what a few parameters can cover. While the parameters used in the proposed models are key to determining general dynamic behaviors of multilingual populations, more parameters may be introduced, and additional mechanisms such as the reward and punishment schemes discussed in [23,24] may be employed.
Computer simulation is a tool to extend the power of theoretical models. It can be used to not only validate the theory but also implement and test additional hypothesis such as the influences from the geographical or demographical structure of the population on multilingual evolution. The simulation done in this study is motivated by previous work on simulation of evolutionary dynamics such as Nowak and May 1992, Durrett and Levin 1997, and more recently, Chen, et al. 2008, Wang, et al. 2012, and Wang, et al. 2016. However, the simulation algorithms can be different for different models and even for the same model. Central to a simulation algorithm is the rule to update the current strategy. Different updating rules may result in very different outcomes. In this work, language competition is considered to be moderate. Therefore, if p A > π, i.e., speaking A benefits more, then x A will stay unchanged if it is already bigger than the average frequency y A . It will increase to y A only when it is less than y A . The rule could be more aggressive such that x A will increase by a certain amount even if it is bigger than y A . The results would look different from what have been presented.
Research on language evolution has been pursued extensively in recent years including modeling and simulation. Work has focused on two different but related aspects of language evolution. One is on how a given language develops and evolves with respect to the "genetic changes" of its linguistic elements such as lexicon, syntax, semantics, pronunciations, etc [30][31][32][33][34][35]. The other is on how a mixed population of multilingual speakers changes dynamically in terms of the percentage of speakers of each language in the whole population [13][14][15][16][17][36][37][38][39][40]. The work described in this paper belongs to the second category. Work on language death by Abrams and Strogatz 2003 has been followed and extended by several research groups [13][14][15]17]. Other approaches have also been proposed, using Lofka-Volterra equations to characterize the evolutionary nature of language dynamics [36,38], introducing reaction-diffusion models to include geographical influences on language competition [37,39], and directly simulating the language dynamics with an agent-based or cellular automata models [16,40]. The model proposed in this work shares some features of these previous approaches, but is focused more on the evolution of multilingualism and in other words, the co-existence of multiple languages in multilingual forms. It is intended to account for possible outside controls or interventions as well as inter-language competitions. Evolutionary games are defined with languages as competing strategies, and corresponding replicator equations are formed to characterize the dynamic behaviors of multilingual populations. The games are also recognized as a special class of potential games and therefore connected to a special class of potential maximization problems, which makes it possible to analyze and even control the equilibrium and stability conditions of given language populations.
To keep the description simple and easy to follow, the proposed models are introduced with and analyzed on bilingual populations only, but they are not limited to bilingual populations. They can actually be extended straightforwardly to general multilingual populations, as given in Appendix E. In order to make the paper more accessible, efforts have been made to keep the discussions from being too technical, with more formal definitions, theorems, and especially proofs all provided as appendices in the end of the paper. Nonetheless, these appendices are considered to be necessary contents as well for supporting the results and conclusions presented in the paper.

Appendix A: Evolutionary games [18, 19]
Let x, y 2 S be the strategies of an individual and a given population, respectively, where S = {x 2 R n : S i x i = 1, x i � 0, i = 1, . . ., n}. Let π(x, y) = S i x i p i (y) be the payoff function for the individual, where p i is the payoff function for the ith pure strategist.
Definition 1 (Evolutionary Game). A Nash equilibrium for the evolutionary game with the payoff function π is a strategy x � such that Definition 2 (Evolutionary Stability). An equilibrium strategy x � is said to be evolutionarily stable if there is a positive number � � < 1 such that for all x 6 ¼ x � and � � � �, Definition 3 (Replicator Equation). The following system of equations is called a system of replicator equations: where π(x, x) = S i x i p i (x).

Theorem 1. The equilibrium strategy of an evolutionary game is a fixed point of the system of replicator equations. It is an asymptotically stable fixed point if it is an evolutionarily stable equilibrium strategy.
Theorem 2. An equilibrium strategy x � is evolutionarily stable if and only if there is a small neighborhood N of x � such that π(x � , x) > π(x, x) for all x 2 N \ S, x 6 ¼ x � .

Appendix B: Potential games [18, 20]
Definition 4 (Potential Game). A game defined by a payoff function π(x, y) = S i x i p i (y) is called a potential game if there is a function f such that @f(y)/@y i = p i (y) for all i = 1, . . ., n.
Definition 5 (Potential Maximization Problem). The following optimization problem is called a potential maximization problem: where S = {y 2 R n : S i y i = 1, y i � 0, i = 1, . . ., n} and @f(y)/@y i = p i (y) for all i = 1, . . ., n. Theorem 3. Let x � be a strategy for a given potential game. Then, x � is an equilibrium strategy for the game if and only if it is a KKT point of the corresponding potential maximization problem.
Theorem 4. Let x � be an equilibrium strategy for a given potential game. Then, x � is evolutionarily stable if and only if it is a strict local maximizer of the corresponding potential maximization problem.

Cases for competition only
Consider the game in (2) with the payoff functions P A and P B being increasing functions as defined in (16). Let ðy � A ; y � B Þ be an equilibrium strategy in one of the three cases given in (7).  (20) at this strategy is y A + y B = 1. The null space of the constraint at this strategy can be represented by a basis matrix z = (1, −1) T . Then, the projected Hessian on this space at ðy � which is always positive definite. Therefore, by the second order necessary conditions for constraint optimization, ðy � A ; y � B Þ can never be a local maximizer of the problem in (20), and therefore, can never be evolutionarily stable.

Cases for intervention only
Consider the game in (5) with the payoff functions � P A and � P B being decreasing functions as defined in (17). Let ðy � A ; y � B Þ be an equilibrium strategy in one of the three cases given in (8).  For case (c), y � A 6 ¼ 0 and y � B 6 ¼ 0. The only active constraint of the problem in (23) at this strategy is y A + y B = 1. The null space of the constraint at this strategy can be represented by a basis matrix z = (1, −1) T . Then, the projected Hessian on this space at ðy � which is always negative definite. Therefore, by the second order sufficient conditions for constraint optimization, ðy � A ; y � B Þ must be a strict local maximizer of the problem in (23), and must therefore be evolutionarily stable.

Appendix D: Dynamics of interventional strategies
Consider the game in (14) with the payoff functionsP A andP B defined in (15), where P A and P B are given in (16) This equation can be written equivalently to It can then be re-arranged to which is equivalent to Define a function ψ for the difference between the left and right hand sides of the stability condition in (33) without the constant terms ð1 À � aÞ and (α − 1): Let y A ¼ y � A , y B ¼ y � B , and substitute � s A and � s B into the function to obtain: The formula can be further simplified to Note that

Appendix E: General models for multilingualism
Assume that there is a multilingual population of m languages. Let x = (x 1 , . . ., x m ) T be the strategy vector for an individual speaker, where x i is the frequency of the individual to speak language i, S i x i = 1. Let y = (y 1 , . . ., y m ) T be the strategy vector for the population, with y i being the average frequency of the population to speak language i, S i y i = 1. Thus, an x-speaker is an individual speaking language i with frequency x i , and a y-population is a population speaking language i with average frequency y i . Assume that the evolution of the given population is influenced only by competition on population and social or economic preferences. Let P i (y i ) be the payoff function for the individual who speaks only i when the average frequency of speaking i in the population is given by y i . Assume that P i , i = 1, . . ., m, are monotonically increasing functions. Then, the payoff function π for an x-speaker in y-population can be defined by the average payoff, π(x, y) = S i x i P i (y i ), with which an evolutionary game can be formulated. A Nash equilibrium of this game is a strategy x � such that π(x � , x � ) � π(x, x � ) for any strategy x. A corresponding system of replicator equations can be defined as _ y i ¼ S j y i y j ðP i ðy i Þ À P j ðy j ÞÞ; i ¼ 1; . . . ; m Similarly, assume that the evolution of the given population is influenced only by possible societal interventions. Let � P i ðy i Þ be the payoff function for the individual who speaks only i when the average frequency of speaking i in the population is given by y i . Assume that � P i , i = 1, . . ., m, are monotonically decreasing functions. Then, the payoff function � p for an x-speaker in y-population can be defined by the average payoff, � pðx; yÞ ¼ S i x i � P i ðy i Þ, with which an evolutionary game can be formulated. A Nash equilibrium of this game is a strategy x � such that � pðx � ; x � Þ � � pðx; x � Þ for any strategy x. A corresponding system of replicator equations can be defined as _ y i ¼ S j y i y j ð � P i ðy i Þ À � P j ðy j ÞÞ; i ¼ 1; . . . ; m ð57Þ LetP i ðy i Þ ¼ lP i ðy i Þ þ ð1 À lÞ � P i ðy i Þ, i = 1, . . ., m, and correspondingly, pðx; yÞ ¼ lpðx; yÞ þ ð1 À lÞ� pðx; yÞ, where λ 2 [0, 1]. Then, a general evolutionary game can be formulated for the population with the influences of both language competition and societal intervention. A Nash equilibrium of the game is a strategy x � such thatpðx � ; x � Þ �pðx; x � Þ for any strategy x. A corresponding system of replicator equations can be defined as _ y i ¼ S j y i y j ðP i ðy i Þ ÀP j ðy j ÞÞ; i ¼ 1; . . . ; m: The payoff functions P i and � P i can be defined more specifically such as P i ðy i Þ ¼ cy aÀ 1 i s i , 1 < α � 2 and � P i ðy i Þ ¼ � cy � aÀ 1 i � s i , 0 � � a < 1, where c and � c are scaling constants, s i are indicators for social or economic preferences, and � s i are language reversing rates due to interventions, 0 � s i ; � s i � 1.