Quantifying the driving factors for language shift in a bilingual region

Significance Languages are an important part of our culturally diverse world, yet many of today’s languages are in danger of dying out. To save endangered languages, one must first understand the dynamics behind language shift: what are the driving factors of people giving up one language for another? Here, we model language dynamics in time and space starting from empirical data. We show that it is the interaction with speakers of the same language that fundamentally determines spread and retreat of a language. This means that a minimum-sized neighborhood of speakers interacting with each other is essential to preserve the language. Many of the world’s around 6,000 languages are in danger of disappearing as people give up use of a minority language in favor of the majority language in a process called language shift. Language shift can be monitored on a large scale through the use of mathematical models by way of differential equations, for example, reaction–diffusion equations. Here, we use a different approach: we propose a model for language dynamics based on the principles of cellular automata/agent-based modeling and combine it with very detailed empirical data. Our model makes it possible to follow language dynamics over space and time, whereas existing models based on differential equations average over space and consequently provide no information on local changes in language use. Additionally, cellular automata models can be used even in cases where models based on differential equations are not applicable, for example, in situations where one language has become dispersed and retreated to language islands. Using data from a bilingual region in Austria, we show that the most important factor in determining the spread and retreat of a language is the interaction with speakers of the same language. External factors like bilingual schools or parish language have only a minor influence.

Many of the world's around 6,000 languages are in danger of disappearing as people give up use of a minority language in favor of the majority language in a process called language shift. Language shift can be monitored on a large scale through the use of mathematical models by way of differential equations, for example, reaction-diffusion equations. Here, we use a different approach: we propose a model for language dynamics based on the principles of cellular automata/agent-based modeling and combine it with very detailed empirical data. Our model makes it possible to follow language dynamics over space and time, whereas existing models based on differential equations average over space and consequently provide no information on local changes in language use. Additionally, cellular automata models can be used even in cases where models based on differential equations are not applicable, for example, in situations where one language has become dispersed and retreated to language islands. Using data from a bilingual region in Austria, we show that the most important factor in determining the spread and retreat of a language is the interaction with speakers of the same language. External factors like bilingual schools or parish language have only a minor influence.
language shift | diffusion | language dynamics | quantitative linguistics | cellular automata I t is estimated that around 90% of the world's 6,000 languages will be replaced by a few dominant languages by the end of the 21st century (1). This replacement, which is called "language shift" (2), leads to a loss of cultural diversity. To prevent this loss and preserve endangered languages, researchers have been trying to find and quantify the factors behind language shift. Language shift (speakers giving up use of one language in favor of another) is driven by a variety of influences, for instance, demographic and social factors (3)(4)(5). To quantify the influence of each of these factors and to study language shift on a large scale, mathematical models and computer simulations have been proposed (6,7). These models generally fall into two categories: (i) macroscopic reaction-diffusion equations that describe the concentration (fraction) of speakers in the population; (ii) microscopic agent-based models that simulate the actions of individual speakers ("agents") changing their language with a certain probability at each interaction. For evaluating both types of model, parameters are required that can be empirically measured so that they can be fitted to data (8). This means that data covering language use over time and space are needed, but such data are often not available in sufficient resolution. Therefore, mathematical models have so far only rarely been checked against data on actual language use.
In this work, we combine mathematical modeling with very detailed empirical data. Applying diffusion theory from physics, we propose a simple model to describe the dynamics of language shift on a microscopic scale based on the principles of cellular automata/ agent-based modeling (9,10). The historical data come from southern Carinthia, Austria, which provides an extremely welldocumented linguistic ecosystem with the interaction of two languages on one and the same territory. Carinthia was a federal state of the Austro-Hungarian Empire until 1918 and of the Federal Republic of Austria afterward. It is geographically separated by a high mountain range, the Karawanks, from the neighbor country Slovenia where Slovenian is the national language. In southern Carinthia, which comprises the districts Klagenfurt and Völkermarkt and parts of the districts Hermagor and Villach (Fig. 1A), the population spoke and speaks partly German and partly Slovenian, the territories being intermixed (11). However, the number of Slovenian speakers in Carinthia has drastically decreased between 1880 and 2001 ( Fig. 1 B and C), and language shift is taking place. We use the data from this case to evaluate our proposed model and its assumptions. Checking against empirical data also allows us to explicitly identify the factors influencing language shift and quantify their impact.

Limits of the Classic Macroscopic Reaction-Diffusion Approach
In the past, language spread and retreat were mostly investigated on a macroscale using differential equations. Macroscopic approaches gained popularity after Abrams and Strogatz (12) published a short seminal paper in 2003 describing the retreat of languages with what they called lower status by differential equations. Their differential equation system considered only temporal and no spatial development, but the paper has drawn a tail of publications in its wake, many of them including spatial development. Spatial and temporal development of languages are usually combined in reaction-diffusion equations (13)(14)(15)(16)(17)(18)(19) of the form ∂u/∂t = D·∂ 2 u/∂x 2 + f(u). These types of equation are also used in other fields, for example, biology or chemistry, to describe all kinds of spread phenomena (20).
Considering a language with higher status, for example, German vs. Slovenian in Carinthia, the development of the fraction u G (x,t) of German speakers in the total population can be written as a

Significance
Languages are an important part of our culturally diverse world, yet many of today's languages are in danger of dying out. To save endangered languages, one must first understand the dynamics behind language shift: what are the driving factors of people giving up one language for another? Here, we model language dynamics in time and space starting from empirical data. We show that it is the interaction with speakers of the same language that fundamentally determines spread and retreat of a language. This means that a minimum-sized neighborhood of speakers interacting with each other is essential to preserve the language. reaction-diffusion equation following Fisher's equation for advantageous genes (21): Here, D G is the diffusivity of German language and k is the conversion rate from Slovenian to German. The fraction of Slovenian speakers is given by u S = 1 − u G . In Carinthia, the language front between German and Slovenian essentially advances only southward motivating a one-dimensional treatment. Eq. 1 then results in a traveling front of the higher status language (in this case German) with velocity v: We defined the front as the line bordering all cells with more than 50% Slovenian speakers each without including outlier cells detached from the contiguous language area ( Fig. 2A). From the data, the velocity of the language front can be derived (Table S1) and the product of diffusivity and conversion rate D G ·k can be determined (Supporting Information). If data on the conversion rate are available, the diffusivity of the majority language can be estimated.
In the period 1971-2001, no contiguous language area exists due to the large decrease in the number of Slovenian speakers (Fig. 2B). The cells with significant fractions of minority language have become dispersed and the minority language has retreated to language islands. Hence there is no continuous language front that clearly separates the two language areas of Slovenian and German. From this, it is evident that the concept of a language front fails when the minority language no longer covers a contiguous area, and reaction-diffusion equations and the resulting language front are not applicable in this case.
Noting the limits of treatment by reaction-diffusion equations, we take a different approach. We simulate the microdevelopment on the basis of the smallest registered population units, hence providing a spatially much more detailed description of language spread and retreat than the macroscale description by reactiondiffusion equations. As a result, we can follow not only global development over time and space but also local processes such as the deviating dynamics in urban areas.

Materials and Methods
We start by taking the data of the censuses in the former Austro-Hungarian Empire that were held in 10-y intervals beginning from 1880 until 1910. Such censuses were continued in the Federal Republic of Austria from 1971 until 2001. In between, the two world wars and the after-war turbulences prevented consistent censuses and data with the same level of detail is not available. In the censuses together with several other items, the vernacular language of each person was asked for and registered. In the censuses until 1910, only one language could be recorded in the questionnaire, whereas in the censuses after 1971 also bilingualism could be indicated. For the sake of simplicity and consistency and data with the same level of detail is not available, we included bilingual Slovenian-German speakers in the minority group, that is, Slovenian speakers. In this paper, we do not consider bilingualism as a separate speaker state. We neither consider the different Slovenian and German dialects that are not encoded in the census data. From the census results, we read the number of German and Slovenian speakers in each of the ∼1,500 population units (hamlets, villages, and towns) in southern Carinthia.
Mapping the Data. For our study of language dynamics in Carinthia, the area investigated is subdivided into a quadratic grid with cells sized 1 × 1 km. Carinthia spans ∼2.4°of longitude or 184 km (east border to west border) and 0.8°of latitude or 84 km (north border to south border). We thus cover the total area of Carinthia with a regular grid of 84 × 184 = 15,456 cells of 1 × 1 km each. Of these 15,456 cells, 9,549 comprise the actual area of Carinthia. A total of 1,170 of the cells inside the Carinthian borders is populated by data for speakers in southern Carinthia. Data for the number of speakers was extracted from census records for the period from 1880 to 1910 (22-25) and the raw dataset for 1971-2001 (special evaluation commissioned from Statistics Austria). Digitized data can be accessed on figshare (https://figshare.com/articles/Language_use_in_Carinthia/ 4535399). Speaker data are assigned to cells based on the geographic coordinates of the population unit. Geographic coordinates (WGS 84) were obtained from KAGIS Kärnten Atlas (https://gis.ktn.gv.at/). The 10-y periods between the censuses are divided into 1-y time steps for the simulation.
Language Shift Model. As a first basic model, we assume that in each cell after 1 y's time the probability p α (r, t + 1) to speak a language α will be proportional to the local number of speakers of that language in the preceding  Diffusivities D α according to simulations for both periods. Changing both diffusivities in the same direction and by the same amount changes little in the quality of the fit, whereas even a small change to the ratio of diffusivities has a big impact on fit quality. A 10% change of both diffusivities in the same direction leads to a 1% change in the sum of absolute errors. On the other hand, if only one diffusivity is changed by 10% (i.e., the ratio between diffusivities changes as well), this leads to a 10% change in the sum of absolute errors.
year n α (r, t) plus an increase through interaction F α (r, t) with speakers of the same language in the neighborhood cells. p α (r, t + 1) is normalized to the total number of speakers and the total interaction in that cell. We obtain Eq. 3 for the probability p α (r, t + 1) to speak language α (where α = G or S, German or Slovenian) in the cell located at position r at time t + 1 p α ðr, t + 1Þ = n α ðr, tÞ + F α ðr, tÞ n S ðr, tÞ + F S ðr, tÞ + n G ðr, tÞ + F G ðr, tÞ . [3] To calculate p α for the first year of the simulation, F α and n α are calculated from the initial census data. Afterward, F α and n α are calculated for each year from the result of the preceding year as follows.
The number of speakers of a language α at position r at time t, n α (r, t), is given by Eq. 4: the probability p α (r, t) to speak the language α at time t multiplied by the total number of people in the cell n total (r, t), which for each time step and cell is given by linear interpolation between censuses: n α ðr, tÞ = n total ðr, tÞ · p α ðr, tÞ. [4] Each interaction term F α (r, t) is a sum over the contributions of all other cells surrounding the initial cell at position r. The interaction F α with speakers of the same language α in the neighboring cells at r j is as follows: F α ðr, tÞ : = F α ðr, n α , tÞ = X rj ≠r c α À r, r j , n α , t Á . [5] The contributions c α (r, r j , n α , t) of all other cells positioned at r j surrounding the initial cell at position r are modeled by Gaussian functions identical to distributions describing the diffusion of particles in physics or chemistry: where D α are the diffusivities of each language, that is, measures for their spread. The diffusivities can also be seen as a measure for the region of influence of a language. We set Δt = 1 y because c α is calculated individually for each year from the result of the preceding year. The Gaussian function is a simple choice to model the interaction with neighboring cells and provides a good fit with the census data. In an extension of our model, this interaction could be modeled by other functions such as leptokurtic (long-tail) distributions or combinations of functions to describe more complex interaction patterns, for example, both long-range and short-range interaction.
Evaluation Procedure. Simulations were performed using GNU Octave 4.0.0. The data from the first census in each period (1880 and 1971) were set as the initial state from which the number of speakers in each cell changes according to Eqs. 3 and 4, assuming a linear population development between censuses. To evaluate the goodness of fit between simulated data and census data, we use ordinary least squares (OLS) to minimize the squared sum of errors: where O i is an observed data point (census data) and E i is an estimated data point (simulated data). t is the number of times the observed data can be compared with the estimated data. t runs from 1 to m − 1, where m is the number of censuses within the period. The data from the first census in each period are excluded as they are equivalent to the initial state of the system; hence there is no error for the initial state and we sum only over the remaining censuses. Optimization was done using the Nelder-Mead method (26). Additionally, we used least absolute errors (LAE) as follows: that is, minimizing the sum of the absolute errors to check the reliability of the fit. General model performance is evaluated by comparison with a baseline. Comparison values can be found in Table S2.

Results
Language Shift in the Period 1880-1910. The widths of both Gaussians, and hence the diffusivities for German and Slovenian, are fitted to the number of speakers in each cell as given by the census data. Fits to the census data were performed for the period from 1880 to 1910. The best solution was achieved with the values given in Table 1. Fig. 3 shows the increase (red) and decrease (blue) of the number of Slovenian speakers in southern Carinthia for census data and simulated data. We obtain satisfactory agreement between the empirical data and the predicted data on a microscopic scale. In detail, this can be seen in Fig. S1 where the model's errors are shown for cells with different numbers of Slovenian speakers. The total number of Slovenian speakers as predicted by the simulation also agrees with the census data (Fig. S2). Fig. S3 shows the residuals (census data minus simulated data). Thus, our model is able to follow how either language has spread and retreated in the time period 1880-1910.
Extension of the Model Through Habitat Parameters. In a second step, for the period from 1880 to 1910, the influence of habitat conditions, such as the influence of urban areas, that is, major towns, the language of schools, and language in parishes, was investigated. To this end, we introduced a habitat parameter h i into Eq. 3, which modifies the effect of local speakers by an exponential function with the argument (±Hh i ). The multiplicative factor H indicates the presence (H = 1) or absence (H = 0) of a local habitat condition, that is, H = 1 for the two largest towns Klagenfurt and Villach or if a bilingual school or Slovenian parish existed. Otherwise, H is set to zero. The exponential function was chosen as a modifier because it is a simple function, which for small h i adds nh i to the speaker effect if H = 1, while recovering the basic model (Eq. 3) if H = 0. We obtain the following equation for the probability p α (r, t + 1) of speaking a language α at position r and time t + 1: n α ðr, tÞ · expð±Hh i Þ + F α ðr, tÞ n S ðr, tÞ · expð+Hh i Þ + F S ðr, tÞ + n G ðr, tÞ · expð−Hh i Þ + F G ðr, tÞ . [9] Optimization was performed as before. Of the three investigated parameters (urban areas, bilingual schools, and parish language), only that of bilingual schools showed a small influence (Supporting Information). However, the influence is so small that Fig. 3B  does not visibly change with the introduction of the bilingualschools habitat parameter.
Language Shift in the Period 1971-2001. After simulating the language dynamics during the Austro-Hungarian Empire, we now turn to language development in the second period after the two world wars. The development from 1971 to 2001 was first pursued using the same basic model (Eqs. 3 and 4). Table 1 shows the numerical results and Fig. 4 A and B shows the increase (red) and decrease (blue) of the number of Slovenian speakers in southern Carinthia for census data and simulation. We also investigated the influence of two habitat parameters for this period: urban areas and parish language. Only the urban habitat parameter resulted in a noticeable difference in goodness of fit ( Fig. 4C and Supporting Information). Errors depending on the number of Slovenian speakers present are shown as before in Fig. S1.
Local Differences in the Language Diffusivities. We have cut out three regions in the districts of Völkermarkt, Klagenfurt, and Villach to search for local differences in the diffusion behavior: is there a difference between rural and urban regions? In all three regions in the districts of Völkermarkt, Klagenfurt, and Villach, the diffusivity of the Slovenian language D S is between 25 and 50% lower than the diffusivity of the German language D G , with the largest difference (factor of 2) for the urban region of Klagenfurt. We suppose that the discrepancies between the different regions are due to local differences in language spread and retreat because of differences in geography and population distribution: faster diffusion in urban areas, German diffusing particularly faster than Slovenian in urban region of Klagenfurt.

Conclusion and Discussion
Macroscopic vs. Microscopic Models. In the past, language dynamics have been commonly described on a macroscopic scale by reaction-diffusion equations that model the fraction of speakers of a language in the population. However, this treatment breaks down when the spread of one language and the retreat of the other one no longer follows a traveling front because one language has become dispersed and has retreated to language islands (Fig. 2B). Additionally, reaction-diffusion equations are not applicable at all in areas without any language front movement such as behind mountain chains (green arrows in Fig. 2A). In contrast, the development can still be followed and predicted in both cases with our microscopic model. The microscopic model also takes into account the interaction with all neighboring cells, whereas in the case of a macroscopic language front the interaction only comes from one direction. Thus, microscopic models yield a more detailed and complete description of language spread and retreat than macroscopic treatment by reactiondiffusion equations.
A challenge for microscopic models on a realistic basis is obviously the need for empirical detailed data (as were at hand for this work) from which to determine the diffusivities. For this reason, language censuses have to be conducted in regular intervals and with fine-grained spatial resolution.
What Drives Language Shift? We have shown that the data predicted by our basic model ( Eqs. 3 and 4) show satisfactory agreement with the historical data for the period between 1880 and 1910. Even in different socioeconomic conditions (the second period between 1971 and 2001), the predicted data still match the empirical data. This means that the basic model can reliably reproduce language dynamics of the studied language competition between Slovenian and German.
The model is also able to reveal similarities between physical phenomena like atomic diffusion and social phenomena like language shift: by modeling linguistic interaction as a Gaussian function as in models of physical diffusion, we obtain good agreement between the predicted and the empirical data. Thus, we have illustrated that it is possible to use physical models to simulate social dynamics on a large scale over time and space.
The basic model uses only two parameters to calculate the probability of speaking a language: the number of speakers in the preceding year and interaction between speakers. Both of these can be directly calculated from census data, ensuring our model is applicable even in situations where data on other factors influencing language use (e.g., perceived status of a language) is not available or even possible to obtain. Without interaction (i.e., using only the number of speakers), the probability of speaking a language (Eq. 3) remains constant. Consequently, interaction with other speakers is an essential drive for the linguistic change in each cell. This point has been argued by linguists (27) and is validated by our simulation. The number of speakers of a language in the population units (hamlets, villages, towns) neighboring the given cell is therefore an important influence on language dynamics. This means that a minimum-sized neighborhood of speakers of the minority language interacting with each other is necessary to preserve the language.
In addition, the simulation shows that other habitat conditions (the language of schools, and in parishes) are of minor influence. There is, however, a noticeable effect of urban areas, which have their own dynamics: between 1880 and 1910, Slovenian decays slightly faster in the larger towns than predicted by the basic model; between 1971 and 2001, the development is reversed, that is, the number of Slovenian speakers increases at a higher rate in large towns than predicted by the basic model (Supporting Information). This reverse in development might be attributed to language playing a larger role in people's identity in an increasingly mobile society (after 1971) compared with a largely rural society (as between 1880 and 1910). When language makes up a larger part of one's identity, there might be a higher tendency to preserve or revive it. This preservation happens, for example, through language associations and cultural clubs, which commonly originate in large towns and consequently have their largest impact there (3). With our model, it is possible to follow these different local developments and quantify the strength of their influence.
As interaction is the driving force for linguistic change in our model, it also offers a tool for possible future work on how interaction shapes language use: what happens when the interaction with speakers of the same language is considerably higher than the interaction with speakers of a different language? How much interaction with the same language (vs. interaction with a different language) is needed for the preservation of the minority language?